text:
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

% Jan 2, 2019
groundTruthFile = '';
expMapFile = 'C:\NeuroSLAM_Datasets\03_NeuroSLAM_Experiments_Results\QUTCarparData\01_exp_map_ml.txt';
% plot_3d_multilayer_experience_map(groundTruthFile, expMapFile, xExpMapScaling, yExpMapScaling, zExpMapScaling, xExpMapTrans, yExpMapTrans, zExpMapTrans, xGtScaling, yGtScaling, zGtScaling)
plot_3d_multilayer_experience_map_qut(groundTruthFile, expMapFile, 0.8, 0.8, 0.1, 0, 0.5, 0, 20,20, 20);

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function circVar=findCircVar(x,w)
%%FINDCIRCVAR Given a set of samples of a circular distribution, determine
%             the circular variance. As in Chapter 2.3.1 of [1], this is
%             just 1 minus the resultant length of the first trigonometric
%             moment.
%INPUTS: x The 1XN or NX1 vector of possibly weighted samples.
%        w The NX1 weights associated with the samples. If this parameter
%          is omitted or an empty matrix is passed, then the samples are
%          uniformly weighted.
%
%OUTPUTS: circVar The circular variance of the measurements.
%
%The mean resultant length is found using findTrigMomentFromSamp and is
%then subtracted from 1.
%
%REFERENCES
%[1] K. V. Mardia and P. E. Jupp, Directional Statistics. Chichester: John
%    Wiley and Sons, 2000.
%
%April 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2)
    w=[];
end

[~,~,RS]=findTrigMomentFromSamp(1,x,w);
circVar=1-RS;

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [Population,FrontNo,CrowdDis] = EnvironmentalSelection(Population,N,zmin,zmax)
% The environmental selection of NSGA-II

%------------------------------- 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".
%--------------------------------------------------------------------------

    %% Normalization
    PopObj = Population.objs - repmat(zmin,length(Population),1);
    range  = zmax - zmin;
    if 0.05*max(range) < min(range)
        PopObj = PopObj./repmat(range,length(Population),1);
    end
    [~,x]      = unique(round(PopObj*1e6)/1e6,'rows');
    PopObj     = PopObj(x,:);
    Population = Population(x);
    N          = min(N,length(Population));

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort_SDR(PopObj,N);
    Next = FrontNo < MaxFNo;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(PopObj,FrontNo);

    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;

    %% Population for next generation
    Population = Population(Next);
    FrontNo    = FrontNo(Next);
    CrowdDis   = CrowdDis(Next);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function name = p19_name ( )

%*****************************************************************************80
%
%% P19_NAME returns the name of problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string NAME, the name of the problem.
%
  name = 'Sqrt-Prod';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function lambda = invol_eigenvalues ( n )

%*****************************************************************************80
%
%% INVOL_LAMBDA returns the eigenvalues of the INVOL matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  k = floor ( n / 2 );

  lambda(1:k,1) =   +1.0;
  lambda(k+1:n,1) = -1.0;

  return
end


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

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

% Last Modified by GUIDE v2.5 03-Dec-2018 17:16:34

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

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

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

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

    lvdData = varargin{1};
    setappdata(hObject,'lvdData',lvdData);

    populateGUI(handles, lvdData);

    % Update handles structure
    guidata(hObject, handles);

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

function populateGUI(handles, lvdData)
    set(handles.moveStageUpButton,'String','<html>↑');
    set(handles.moveStageDownButton,'String','<html>↓');
    set(handles.stagesListbox,'String',lvdData.launchVehicle.getStagesListBoxStr());

    numStages = length(get(handles.stagesListbox,'String'));
    if(numStages <= 1)
        handles.removeStageButton.Enable = 'off';
    end

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

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = [];
    else
        close(handles.lvd_EditStagesGUI);
    end

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

% --- Executes on selection change in stagesListbox.
function stagesListbox_Callback(hObject, eventdata, handles)
% hObject    handle to stagesListbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns stagesListbox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from stagesListbox
    if(strcmpi(get(handles.lvd_EditStagesGUI,'SelectionType'),'open'))
        lvdData = getappdata(handles.lvd_EditStagesGUI,'lvdData');
        lv = lvdData.launchVehicle;

        selStage = get(handles.stagesListbox,'Value');
        stage = lv.getStageForInd(selStage);

        lvd_EditStageGUI(stage);
        set(handles.stagesListbox,'String',lvdData.launchVehicle.getStagesListBoxStr());
    end

% --- Executes during object creation, after setting all properties.
function stagesListbox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to stagesListbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in addStageButton.
function addStageButton_Callback(hObject, eventdata, handles)
% hObject    handle to addStageButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditStagesGUI,'lvdData');
    lv = lvdData.launchVehicle;

    stage = LaunchVehicleStage(lv);
    useStage = lvd_EditStageGUI(stage);

    if(useStage)
        lv.addStage(stage);

        stageState = LaunchVehicleStageState(stage);
        lvdData.initStateModel.addStageState(stageState);

        set(handles.stagesListbox,'String',lvdData.launchVehicle.getStagesListBoxStr());

        handles.removeStageButton.Enable = 'on';
    end

% --- Executes on button press in removeStageButton.
function removeStageButton_Callback(hObject, eventdata, handles)
% hObject    handle to removeStageButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditStagesGUI,'lvdData');
    lv = lvdData.launchVehicle;

    selStage = get(handles.stagesListbox,'Value');
    stage = lv.getStageForInd(selStage);

    tf = stage.isStageAndChildrenInUse();

    if(tf == false)
        lvdData.initStateModel.removeStageStateForStage(stage);
        lv.removeStage(stage);

        set(handles.stagesListbox,'String',lvdData.launchVehicle.getStagesListBoxStr());

        numStages = lv.getNumStages();
        if(selStage > numStages)
            set(handles.stagesListbox,'Value',numStages);
        end

        if(numStages <= 1)
            handles.removeStageButton.Enable = 'off';
        end
    else
        warndlg(sprintf('Could not delete the stage "%s" because it (or an associated engine or tank) is in use as part of an event termination condition, event action, objective function, or constraint.  Remove all tanks and engines from the stage and remove the stage dependencies before attempting to delete the stage.', stage.name),'Cannot Delete Stage','modal');
    end

% --- Executes on button press in moveStageUpButton.
function moveStageUpButton_Callback(hObject, eventdata, handles)
% hObject    handle to moveStageUpButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditStagesGUI,'lvdData');
    lv = lvdData.launchVehicle;

    selStage = get(handles.stagesListbox,'Value');
    stage = lv.getStageForInd(selStage);

    lv.moveStageUp(stage);

    if(selStage>1)
        set(handles.stagesListbox,'Value',selStage-1);
    end

    set(handles.stagesListbox,'String',lvdData.launchVehicle.getStagesListBoxStr());

% --- Executes on button press in moveStageDownButton.
function moveStageDownButton_Callback(hObject, eventdata, handles)
% hObject    handle to moveStageDownButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditStagesGUI,'lvdData');
    lv = lvdData.launchVehicle;

    selStage = get(handles.stagesListbox,'Value');
    stage = lv.getStageForInd(selStage);

    lv.moveStageDown(stage);

    numStages = lv.getNumStages();
    if(selStage<numStages)
        set(handles.stagesListbox,'Value',selStage+1);
    end

    set(handles.stagesListbox,'String',lvdData.launchVehicle.getStagesListBoxStr());

% --- Executes on key press with focus on lvd_EditStagesGUI or any of its controls.
function lvd_EditStagesGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditStagesGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            uiresume(handles.lvd_EditStagesGUI);
        case 'enter'
            uiresume(handles.lvd_EditStagesGUI);
        case 'escape'
            uiresume(handles.lvd_EditStagesGUI);
    end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%   TTeMPS Toolbox.
%   Michael Steinlechner, 2013-2016
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

function check_precond_laplace(P, grad, P_grad)

d = grad.order;

% first check whether P_grad is correctly gauged
for ii=1:d-1
    orth_err(ii) = norm( unfold(P_grad.dU{ii},'left')'*unfold(grad.U{ii},'left') , 2 ) / norm( unfold(P_grad.dU{ii},'left'), 2 );
end

if max(orth_err) > 1e-12
    orth_err
    warning('precond_laplace.m returned an inaccurate result')
end

end

author: NicolasBoumal

repo: manopt


text:
function matrixData=ivech(stackedData)
% Transform a vector in to a symmetric matrix for use in covariance applications
%
% USAGE:
%   MATRIXDATA = ivech(STACKEDDATA)
%
% INPUTS:
%   STACKEDDATA:   A K(K+1)/2 vector of data to be transformed
%
% OUTPUTS:
%   MATRIXDATA   - a K by K symmetric matrix of the form
%                  [ data(1) data(2)    data(3)     ...               data(K)
%                    data(2) data(K+1)  data(K+2)   ...               ...
%                    data(3) data(K+2)  data(2K)    ...               ...
%                    ...     ....       ...         ...               data(K(K+1)/2-1)
%                    data(K) data(2K-1) ...         data(K(K+1)/2-1)  data(K(K+1)/2) ]
%
% COMMENTS:
%
% See also VECH

% Author: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 2/1/2008

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if size(stackedData,2)>size(stackedData,1)
    stackedData=stackedData';
end

if size(stackedData,2)~=1
    error('STACKED_DATA must be a column vector.')
end

K2=size(stackedData,1);
K=(-1+sqrt(1+8*K2))/2;

if floor(K)~=K
    error(['The number of elemeents in STACKED_DATA must be conformable to' ...
    'the inverse vech operation.'])
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize the output data
matrixData=zeros(K);

% Use a logical trick to inverse the vech
pl=tril(true(K));
matrixData(pl)=stackedData;
diag_matrixData=diag(diag(matrixData));
matrixData=matrixData+matrixData'-diag_matrixData;

author: bashtage

repo: mfe-toolbox


text:
function wtime_test05 ( )

%*****************************************************************************80
%
%% TEST05 times the matrix multiplication problem problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2010
%
%  Author:
%
%    John Burkardt
%
  l_log_min = 1;
  l_log_max = 5;
  l_min = 4^l_log_min;
  l_max = 4^l_log_max;
  rep_num = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST05\n' );
  fprintf ( 1, '  Time the matrix multiplication problem.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute C = A * B\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  where\n' );
  fprintf ( 1, '    A is an L by M matrix,\n' );
  fprintf ( 1, '    B is an M by N matrix,\n' );
  fprintf ( 1, '  and so\n' );
  fprintf ( 1, '    C is an L by N matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Minimum value of L = M = N = %d\n', l_min );
  fprintf ( 1, '  Maximum value of L = M = N = %d\n', l_max );
  fprintf ( 1, '  Number of repetitions of the operation: %d\n', rep_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use triple loop for matrix multiplication.\n' );

  for rep = 1 : rep_num

    for l_log = l_log_min : l_log_max - 1

      l = 4^( l_log );
      m = l;
      n = l;

      a = rand ( l, l );
      b = rand ( l, l );

      seconds = wtime ( );

      for i = 1 : l
        for j = 1 : l
          c(i,j) = 0.0;
          for k = 1 : l
            c(i,j) = c(i,j) + a(i,k) * b(k,j);
          end
        end
      end

      delta(l_log,rep) = wtime ( ) - seconds;

    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Timing results using nested loops:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Vector Size  Rep #1        Rep #2        Rep #3        ' );
  fprintf ( 1, 'Rep #4        Rep #5\n' );
  fprintf ( 1, '\n' );
  for l_log = l_log_min : l_log_max - 1
    l = 4^( l_log );
    fprintf ( 1, '  %8d  %12f  %12f  %12f  %12f  %12f\n', ...
      l, delta(l_log,1:rep_num) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use C = A * B matrix multiplication.\n' );

  for rep = 1 : rep_num

    for l_log = l_log_min : l_log_max - 1

      l = 4^( l_log );
      m = l;
      n = l;

      a = rand ( l, l );
      b = rand ( l, l );

      seconds = wtime ( );

      c = a * b;

      delta(l_log,rep) = wtime ( ) - seconds;

    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Timing results using C = A * B:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Vector Size  Rep #1        Rep #2        Rep #3        ' );
  fprintf ( 1, 'Rep #4        Rep #5\n' );
  fprintf ( 1, '\n' );
  for l_log = l_log_min : l_log_max - 1
    l = 4^( l_log );
    fprintf ( 1, '  %8d  %12f  %12f  %12f  %12f  %12f\n', ...
      l, delta(l_log,1:rep_num) );
  end

  return
end


text:
function Lambda = MakeFourierDiagonal_2D(n,deep)
% MakeFourierDiagonal_2D:
%  Usage:
%    Lambda = Inv_SeparateAngles(R,deep,MaxIts,ErrTol);
%  Inputs:
%    n       dyadic integer
%    deep    Depth of the angular splitting
%  Outputs:
%    Lambda   (2*n-1) * n/4 * 2 array of Fourier multipliers
%  Description
%    AtA has a Toeplitz structure. A Toeplitz matrix is embedded in
%    a larger circulant matrix. A circulant matrix is diagonal in a
%    Fourier basis; lambda is the vector of those diagonal
%    coefficients.
%  See Also
%    AtA_Toeplitz, MakeFourierDiagonal
%
% By Emmanuel Candes, 2003-2004

	 n2 = n/2;
         boxlen = n/2^deep;
	 boxcnt = 2^deep;

	 [ix,w] = DetailMeyerWindow([n2/4 n2/2],3);
	 lx = reverse(-ix);

	 m      =  0:(boxcnt - 1);
         ym     =  (m-boxcnt/2).*boxlen + boxlen/2;
	 slope  =  -ym./n2;

	 Lambda = zeros(2*n-1,n2/2,2);

    for r = 1:(n2/2),
      k = lx(r);
      shift  =  ym + slope.*(k+n2);
      alpha =  -k/n2;
      w = MakeSineWindow(boxlen,alpha);
      lambda = MakeFourierDiagonal(n,shift,boxlen,1,w);

      Lambda(:,r,1) = lambda;

      rsym = n2/2 + 1 - r; %shift = -shift + 1;
                           %Lambda(:,rsym,2) = MakeFourierDiagonal(n,shift,boxlen,1,w);
      Lambda(:,rsym,2) =  ifft(conj(fft(lambda)));
    end

author: thomaskuestner


text:
function updatelarvacontour(hfly,hfly_extra,pos)

if isempty(pos.xcontour),
  set(hfly,'XData',[],'YData',[]);
else
  set(hfly,'XData',[pos.xcontour(:);pos.xcontour(1)],'YData',[pos.ycontour(:);pos.ycontour(1)]);
end
if all(isfield(pos,{'xspine','yspine'})),
  set(hfly_extra(1),'XData',pos.xspine,'YData',pos.yspine);
  set(hfly_extra(2),'XData',pos.xspine(1),'YData',pos.yspine(1));
end

author: kristinbranson

repo: JAABA


text:
function  res = p2DFT(mask,imSize,ph,mode)

%res = p2DFT(mask,imSize [ ,phase,mode])
%
%
%	Implementation of partial Fourier operator.
%
%	input:
%			mask - 2D matrix with 1 in entries to compute the FT and 0 in ones tha
%				are not computed.
%			imSize - the image size (1x2)
%			phase - Phase of the image for phase correction
%			mode - 1- real, 2-cmplx
%
%	Output:
%			The operator
%
%	(c) Michael Lustig 2007

if nargin <3
	ph = 1;
end
if nargin <4
	mode = 2; % 0 - positive, 1- real, 3-cmplx
end

res.adjoint = 0;
res.mask = mask;
res.imSize = imSize;

res.dataSize = [size(mask)];
res.ph = ph;
res.mode = mode;
res = class(res,'p2DFT');


text:
function status = lclear(handles)

% tstool/lclear
% function status = lremove(handles)
% Loeschen aller Eintraege aus der Listbox und den UserData-Field der Figure, aber keine physik. Entfernung
% der Files

lhandle = handles.lboxhandle;

String = {};
set(lhandle, 'Value', 1);
set(lhandle, 'String', String);
setcurrentfile('', handles.lboxhandle, handles.currfilehandle); % jetzt ist kein File mehr der current working file
status = 0;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%io_writelcm.m
%Jamie Near, McGill University 2014.
%
% USAGE:
% RF=io_writelcm(in,outfile,te);
%
% DESCRIPTION:
% Takes MRS data in matlab structure format and writes it to a text file
% that can be read by LCModel.
%
% INPUTS:
% in         = input data in matlab structure format.
% outfile    = Desired filename of output text file.
% te         = Echo time of acquisition (in ms).
%
% OUTPUTS:
% RF         = Same as input.  Not used.  The primary output of this
%                function is a text file in LCModel raw format.

function RF=io_writelcm(in,outfile,te);
%function RF=writelcm(in,outfile,te);

if in.flags.isFourSteps
    error('ERROR:  Must first combine four subspecs using op_fourStepCombine');
end

if in.flags.subtracted
    %error('ERROR:  This operation must be done prior to combining subspecs');
end

if ~in.flags.averaged
    disp('WARNING:  Signals must be averaged first');
end

if ~in.flags.addedrcvrs
    error('ERROR:  reciever channels must be combined first');
end

% datsets=1;
% zop=0;
% t0=0;
  Bo=in.Bo;
  hzppm=42.577*Bo;
  dwelltime=in.dwelltime;
% Nuc=0;
% PatName='No Name';
% scanner='TrioTim';
% addinfo='jnear';
seq='PRESS';

RF=zeros(in.sz(in.dims.t),2);
RF(:,1)=real(in.fids(:,1));
RF(:,2)=-imag(in.fids(:,1));

%write to txt file for jmrui
fid=fopen(outfile,'w+');
fprintf(fid,' $SEQPAR');
%fprintf(fid,'\n\nFilename: %s' ,outfile);
%fprintf(fid,'\n\nPointsInDataset: %i',length(data_struct.fids));
%fprintf(fid,'\nDatasetsInFile: %i',datsets);
%fprintf(fid,'\nSamplingInterval: %1.0E',data_struct.dwelltime*1000);
%fprintf(fid,'\nZeroOrderPhase: %1.0E',zop);
%fprintf(fid,'\nBeginTime: %1.0E',t0);
%fprintf(fid,'\nTransmitterFrequency: %2.4E',data_struct.txfrq);
%fprintf(fid,'\nMagneticField: %2.1E',Bo);
%fprintf(fid,'\nTypeOfNucleus: %1.0E',Nuc);
%fprintf(fid,'\nNameOfPatient: %s',PatName);
%fprintf(fid,'\nDateOfExperiment: %i',data_struct.date);
%fprintf(fid,'\nSpectrometer: %s',scanner);
%fprintf(fid,'\nAdditionalInfo: %s\n\n\n',addinfo);
%fprintf(fid,'Signal and FFT\n');
%fprintf(fid,'sig(real)\tsig(imag)\tfft(real)\tfft(imag)\n');
%fprintf(fid,'Signal 1 out of %i in file\n',datsets);
fprintf(fid,'\n echot= %2.2f',te);
fprintf(fid,'\n seq= ''PRESS''');
fprintf(fid,'\n hzpppm= %5.6f',in.txfrq/1e6);
fprintf(fid,'\n NumberOfPoints= %i',in.sz(1));
fprintf(fid,'\n dwellTime= %5.6f' ,dwelltime);
fprintf(fid,'\n $END');
fprintf(fid,'\n $NMID');
fprintf(fid,'\n id=''ANONYMOUS '', fmtdat=''(2E15.6)''');
fprintf(fid,'\n volume=8.0');
fprintf(fid,'\n tramp=1.0');
fprintf(fid,'\n $END\n');
fprintf(fid,'  % 7.6e  % 7.6e\n',RF');
fclose(fid);

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [tfr,rtfr,hat] = tfrrgab(x,t,N,Nh,trace,K);
%TFRRGAB Reassigned Gabor spectrogram time-frequency distribution.
%	[TFR,RTFR,HAT] = TFRRGAB(X,T,N,NH,TRACE,K)
%	computes the Gabor spectrogram and its reassigned version.
%	This particular window (a Gaussian window) allows a 20 % faster
%	algorithm than the TFRRSP function.
%
%	X     : analysed signal
%	T     : the time instant(s)           (default : 1:length(X))
%	N     : number of frequency bins      (default : length(X))
%	NH    : length of the gaussian window (default : N/4))
%	TRACE : if nonzero, the progression of the algorithm is shown
%                                             (default : 0).
%	K     : value at both extremities     (default 0.001)
%	TFR,  : time-frequency representation and its reassigned
%	RTFR    version. When called without output arguments,
%	        TFRRGAB runs TFRQVIEW.
%	HAT   : Complex matrix of the reassignment vectors.
%
%	Example :
%	 sig=fmlin(128,0.1,0.4); tfrrgab(sig,1:128,128,19,1);
%
%	See also all the time-frequency representations listed in
%	 the file CONTENTS (TFR*)

%	F. Auger, May-July 1994, July 1995.
%       Copyright (c) 1996 by CNRS(France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('At least 1 parameter required');
end;
[xrow,xcol] = size(x);
if (nargin <= 2),
 N=xrow;
end;

hlength=floor(N/4);
hlength=hlength+1-rem(hlength,2);

if (nargin == 1),
 t=1:xrow;
end;

if (nargin <= 3),
 Nh=hlength; trace=0; K=0.001;
elseif (nargin == 4),
 trace = 0; K=0.001;
elseif (nargin == 5),
 K= 0.001;
end;

if (N<0),
 error('N must be greater than zero');
end;
[trow,tcol] = size(t);
if (xcol~=1),
 error('X must have only one column');
elseif (trow~=1),
 error('T must only have one row');
elseif (2^nextpow2(N)~=N & nargin==6),
 fprintf('For a faster computation, N should be a power of two\n');
end;

if (rem(Nh,2)==0),
 error('Nh must be odd');
elseif length(Nh)~=1,
 error('Nh must be a scalar');
end;

Nh2=Nh-2;
TFTBcontinue=1;
while TFTBcontinue,
 Nh2=Nh2+2;
 h=tftb_window(Nh2,'gauss',K^((Nh2-1)^2 /(Nh-1)^2));
 TFTBcontinue=(h(Nh2)*(Nh2-1)>2*K);
end;

K=K^((Nh2-1)^2 /(Nh-1)^2); Nh=Nh2; Lh=(Nh-1)/2;
h=h; Th=h.*[-Lh:Lh]';

if (tcol==1),
 Dt=1;
else
 Deltat=t(2:tcol)-t(1:tcol-1);
 Mini=min(Deltat); Maxi=max(Deltat);
 if (Mini~=Maxi),
  error('The time instants must be regularly sampled.');
 else
  Dt=Mini;
 end;
 clear Deltat Mini Maxi;
end;

tfr= zeros(N,tcol); tf2= zeros(N,tcol); tf3= zeros(N,tcol);
if trace, disp('Gabor spectrogram'); end;

for icol=1:tcol,
 if trace, disprog(icol,tcol,10); end;
 ti= t(icol);
 tau=-min([round(N/2)-1,Lh,ti-1]):min([round(N/2)-1,Lh,xrow-ti]);
 indices= rem(N+tau,N)+1;
 norm_h=norm(h(Lh+1+tau));
 tfr(indices,icol)=x(ti+tau).*conj( h(Lh+1+tau)) /norm_h;
 tf2(indices,icol)=x(ti+tau).*conj(Th(Lh+1+tau)) /norm_h;
end ;
tfr=fft(tfr); tf2=fft(tf2);
tfr=tfr(:); tf2=tf2(:);  tf3=tf3(:);
avoid_warn=find(tfr~=0.0);
tf3(avoid_warn)=round(imag(2*log(K)*N*tf2(avoid_warn)./tfr(avoid_warn)/(2.0*pi*Lh^2)));
tf2(avoid_warn)=round(real(tf2(avoid_warn)./tfr(avoid_warn)/Dt));
tfr=abs(tfr).^2;
if trace, fprintf ('\nreassignment: \n'); end;
tfr=reshape(tfr,N,tcol);
tf2=reshape(tf2,N,tcol);
tf3=reshape(tf3,N,tcol);

rtfr= zeros(N,tcol);
Ex=mean(abs(x(min(t):max(t))).^2); Threshold=1.0e-6*Ex;
for icol=1:tcol,
 if trace, disprog(icol,tcol,10); end;
 for jcol=1:N,
  if abs(tfr(jcol,icol))>Threshold,
   icolhat= icol + tf2(jcol,icol);
   icolhat=min(max(icolhat,1),tcol);
   jcolhat= jcol - tf3(jcol,icol);
   %while (jcolhat<1),jcolhat=jcolhat+N; end;
   %while (jcolhat>N),jcolhat=jcolhat-N; end;
   jcolhat=rem(rem(jcolhat-1,N)+N,N)+1;
   rtfr(jcolhat,icolhat)=rtfr(jcolhat,icolhat) + tfr(jcol,icol) ;
   tf2(jcol,icol)=jcolhat + j * icolhat;
  else
   tf2(jcol,icol)=inf*(1+j);
   rtfr(jcol,icol)=rtfr(jcol,icol) + tfr(jcol,icol) ;
  end;
 end;
end;

if trace, fprintf('\n'); end;
clear tf3;
if (nargout==0),
 TFTBcontinue=1;
 while (TFTBcontinue==1),
  choice=menu ('Choose the representation:',...
               'stop',...
               'Gabor spectrogram',...
               'reassigned Gabor spectrogram');
  if (choice==1), TFTBcontinue=0;
  elseif (choice==2),
   Q=round(tcol*N/xrow);
   tfrqview(tfr,x,t,'tfrgabor',tcol,Q,h);
  elseif (choice==3),
   tfrqview(rtfr,x,t,'tfrrgab',Nh);
  end;
 end;
elseif (nargout>2),
 hat=tf2;
end;

author: HeLiangHIT


text:
function [coordinate, rgb] = getPointsFromRGBD(RGBDframe)

X = reshape(RGBDframe(:,:,[1 2 3 5 6 4]),[],6)';
valid = X(6,:) > 0;

coordinate = X(4:6,valid);
rgb = X(1:3,valid);

author: jianxiongxiao

repo: ProfXkit


text:
function [X] = spm_en(X,p)
% Euclidean normalization
% FORMAT [X] = spm_en(X,[p]);
% X   - matrix
% p   - optional polynomial detrend [default = []]
%__________________________________________________________________________
%
% spm_en performs a Euclidean normalization setting the column-wise sum of
% squares to unity (leaving columns of zeros as zeros)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_en.m 3901 2010-05-27 16:14:36Z karl $

% detrend
%--------------------------------------------------------------------------
if nargin > 1
    X = spm_detrend(X,p);
end

% Euclidean normalization
%--------------------------------------------------------------------------
for i = 1:size(X,2)
    if any(X(:,i))
        X(:,i) = X(:,i)/sqrt(sum(X(:,i).^2));
    end
end

author: spm


text:
%% Documentation
%
%%
% Use the left sidebar to navigate through the different topics. Note that
% all documentation pages are essentially MTEX scripts. You have a copy of
% them in your MTEX download, with a file name equal to the last part of
% the URL with .html swapped for .m (so this page is |Documentation.m|).
% This means you can open them in Matlab by typing
%
%   edit fileName
%
% This allows you to execute the sections by step by step and investiage
% the effect of each single MTEX command. It also should motivate you to
% play around with different parameters, different input files or further
% analysis steps. Eventually, these documentation scripts may serve as a
% good starting point for your own analysis scripts.
%
%% Contributing to the Documentation
%
% Documenting a project like MTEX is a challenging and ongoing task.
% Therefore we are extremely happy about any spelling fixes, examples,
% theoretical explainations, special use cases, etc. As a bonus everybody
% who contributed to MTEX will automatically appear at the
% <https://github.com/mtex-toolbox/mtex/graphs/contributors contributors
% page>.
%
% The easiest way to suggest changes to the MTEX documentation is online
% via Github. This is done as follows:
%
% # sign into <https://github.com GitHub>
% # open the help page you want to change, e.g., this one
% # click
% <https://github.com/mtex-toolbox/mtex/blob/develop/doc/Documentation.m
% edit page> which leads you to the file on <https://github.com GitHub>
% # click the little pencil in the top right corner which allows you to
% edit the script online
% # enter or copy and paste your changes into the site
% # click at the bottom *Propose changes*
% # click create *pull request*
%
% Alternatively, you can send us your changes simply by <mailto:ralf.hielscher@mathematik.tu-chemnitz.de email>.
%
%% The Publish Command
%
% You may get a visual preview of your altered documentation page by using the
% matlab function
% <https://de.mathworks.com/help/matlab/matlab_prog/publishing-matlab-code.html publish>
%
%   publish filename
%
% which creates an folder, named |html|, in your current directory where you can
% find the generated html page.
%

author: mtex-toolbox

repo: mtex

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


text:
function results = vl_test_pegasos(varargin)
% VL_TEST_KDTREE
vl_test_init ;

function s = setup()
randn('state',0) ;

s.biasMultiplier = 10 ;
s.lambda = 0.01 ;

Np = 10 ;
Nn = 10 ;
Xp = diag([1 3])*randn(2, Np) ;
Xn = diag([1 3])*randn(2, Nn) ;
Xp(1,:) = Xp(1,:) + 2 + 1 ;
Xn(1,:) = Xn(1,:) - 2 + 1 ;

s.X = [Xp Xn] ;
s.y = [ones(1,Np) -ones(1,Nn)] ;
%s.w = exact_solver(s.X, s.y, s.lambda, s.biasMultiplier)
s.w = [1.181106685845652 ;
       0.098478251033487 ;
       -0.154057992404545 ] ;

function test_problem_1(s)
for conv = {@single,@double}
  vl_twister('state',0) ;
  conv = conv{1} ;
  w = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                 'NumIterations', 100000, ...
                 'BiasMultiplier', s.biasMultiplier, ...
                 'Preconditioner', conv([1 1 .1])) ;
  vl_assert_almost_equal(w, conv(s.w), 0.1) ;
end

function test_continue_training(s)
for conv = {@single,@double}
  conv = conv{1} ;

  vl_twister('state',0) ;
  w = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                 'NumIterations', 3000, ...
                 'BiasMultiplier', s.biasMultiplier) ;

  vl_twister('state',0) ;
  w1 = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                 'StartingIteration', 1, ...
                 'NumIterations', 1500, ...
                  'BiasMultiplier', s.biasMultiplier) ;
  w2 = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                  'StartingIteration', 1501, ...
                  'StartingModel', w1, ...
                  'NumIterations', 1500, ...
                  'BiasMultiplier', s.biasMultiplier) ;
  vl_assert_almost_equal(w,w2,1e-7) ;
end

function test_continue_training_with_perm(s)
perm = uint32(randperm(size(s.X,2))) ;
for conv = {@single,@double}
  conv = conv{1} ;

  vl_twister('state',0) ;
  w = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                 'NumIterations', 3000, ...
                 'BiasMultiplier', s.biasMultiplier, ...
                 'Permutation', perm) ;

  vl_twister('state',0) ;
  w1 = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                 'StartingIteration', 1, ...
                 'NumIterations', 1500, ...
                 'BiasMultiplier', s.biasMultiplier, ...
                  'Permutation', perm) ;
  w2 = vl_pegasos(conv(s.X), int8(s.y), s.lambda, ...
                  'StartingIteration', 1501, ...
                  'StartingModel', w1, ...
                  'NumIterations', 1500, ...
                  'BiasMultiplier', s.biasMultiplier, ...
                  'Permutation', perm) ;
  vl_assert_almost_equal(w,w2,1e-7) ;
end

function w = exact_solver(X, y, lambda, biasMultiplier)
N = size(X,2) ;
model = svmtrain(y', [(1:N)' X'*X], sprintf(' -c %f -t 4 ', 1/(lambda*N))) ;
w = X(:,model.SVs) * model.sv_coef ;
w(3) = - model.rho / biasMultiplier ;
format long ;
disp('model w:')
disp(w)


text:
function [id,OK] = selector(idList,NameList,Caption)
% [id,OK] = selector(idList,NameList,Caption,view)
% [id,OK] = selector(idList,NameList,Caption)
% [id,OK] = selector(idList,NameList)
% User should pick one item from a list

if(nargin==2)
  Caption = 'Select:';
end

c=0;

c=c+1;
uiStruct(c).string = Caption;
uiStruct(c).fieldName = 'main';
uiStruct(c).list = NameList;
uiStruct(c).choice = 1;
uiStruct(c).style = 'popupmenu';

height = 1;
vSkip = 0.3;
pos = [35,10,100,length(uiStruct)*(height+vSkip)+3];
x = 1;
y = length(uiStruct)*(height+vSkip)+1;
editWidth = 50;
stringWidth = pos(3)-editWidth-2;
for uiNum = 1:length(uiStruct)
  uiStruct(uiNum).stringPos = [x,y,stringWidth,height];
  uiStruct(uiNum).editPos = [x+stringWidth-1,y,editWidth,height];
  uiStruct(uiNum).value = 1;
  y = y-(height+vSkip);
end
outStruct = generaldlg(uiStruct,pos);

if ~isempty(outStruct)
  id = idList(strcmp(NameList,outStruct.main));
  OK = 1;
else
  id = 0;
  OK = 0;
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function pde = sincosNeumanndata
%% SINCOSNEUMANNDATA trigonometric  data for Poisson equation
%
%     f = 8*pi^2*sin(2*pi*x)*cos(2*pi*y);
%     u = sin(2*pi*x)*cos(2*pi*y);
%     Du = (2*pi*cos(2*pi*x)*cos(2*pi*y), -2*pi*sin(2*pi*x)*sin(2*pi*y));
%     du/dn = g_N on [0,1]^2.
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

pde = struct('f',@f,'exactu',@exactu,'Du',@Du,'g_D',@g_D,'g_N',@g_N);

% load data (right hand side function)
function rhs =  f(p)
    x = p(:,1); y = p(:,2);
    rhs =  8*pi^2*sin(2*pi*x).*cos(2*pi*y);
end
% exact solution
function u =  exactu(p)
    x = p(:,1); y = p(:,2);
    u =  sin(2*pi*x).*cos(2*pi*y);
end
% derivative of the exact solution
function uprime =  Du(p)
    x = p(:,1); y = p(:,2);
    uprime(:,1) = 2*pi*cos(2*pi*x).*cos(2*pi*y);
    uprime(:,2) = -2*pi*sin(2*pi*x).*sin(2*pi*y);
end
function u = g_D(p)
    u = exactu(p);
end
% Neumann boundary condition
function f = g_N(p,vargin)
    if nargin > 1
        f = dot(Du(p),vargin,2);
    else
        f = zeros(size(p,1),1);
        x = p(:,1); y = p(:,2);
        uprime = [2*pi*cos(2*pi*x).*cos(2*pi*y) -2*pi*sin(2*pi*x).*sin(2*pi*y)];
        leftbd = (abs(x)<eps);  % n = (-1,0);
        f(leftbd) = - uprime(leftbd,1);
        rightbd = (abs(x-1)<eps); % n = (1,0);
        f(rightbd) = uprime(rightbd,1);
        topbd = (abs(y-1)<eps);   % n = (0,1)
        f(topbd) = uprime(topbd,2);
        bottombd = (abs(y)<eps);% n = (0,-1)
        f(bottombd) = - uprime(bottombd,2);
    end
end
end

repo: ifem


text:
function y = mg(DataDuration)

% Mackey-Glass chaotic time series generator.

NumOfSamples = DataDuration;
discarded = 350;
x = zeros(NumOfSamples + discarded + 1);
y = zeros(1,NumOfSamples);    % Each column represents different time instants.
Ts = 1;
x(1) = 1.2;
tau  = 17;

for n = tau+1:(NumOfSamples + discarded)
   x(n+1) = Ts*(0.2*x(n-tau)/( 1 + (x(n-tau))^10)) + (1 - 0.1*Ts)*x(n);
end

% plot(x);
% title('Mackey-Glass Chaotic Time Series');
% xlabel('Time Sample');
% ylabel('Estimated Value');

for k=1:NumOfSamples
   y(k) = x(discarded+k);
end

author: Sable

repo: mcbench-benchmarks


text:
function [mu, sigma, c] = ...
    gmm_estimate_of_speaker(features, num_of_gaussians, num_of_iter)
    [mu, sigma, c] = GMMImpl.gmm_estimate(features, num_of_gaussians, num_of_iter);

%     obj = gmdistribution.fit(speaker_features', num_of_gaussians, ...
%         'CovType', 'diagonal', 'Regularize', 1);
%     mu = obj.mu';
%     sigma = squeeze(obj.Sigma);
%     c = obj.PComponents;
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function s = logsumexp(a, dim, method)
% Returns log(sum(exp(a),dim)) while avoiding numerical underflow.
% Default is dim = 1 (columns).
% logsumexp(a, 2) will sum across rows instead of columns.
% Unlike matlab's "sum", it will not switch the summing direction
% if you provide a row vector.
%
% Example: s = logsumexp(rand(10,3),2)

% Written by Tom Minka
% (c) Microsoft Corporation. All rights reserved.

if nargin < 2, dim = 2; end
if nargin < 3, method = 1; end

% subtract the largest in each column
y = max(a,[],dim);
if method==1
  dims = ones(1,ndims(a));
  dims(dim) = size(a,dim);
  a = a - repmat(y, dims);
else
  % Added by KPM. Just a hair faster (see timing comparison below).
  % A=rand(1000,100);N=100;
  % tic; for i=1:N, s1 = logsumexp(A,2,1); end; t1=toc;
  % tic; for i=1:N, s2=logsumexp(A,2,2); end; t2=toc;
  % assert(approxeq(s1,s2)); [t1 t2]
  a = bsxfun(@minus, a, y);
end

s = y + log(sum(exp(a),dim));
i = find(~isfinite(y));
if ~isempty(i)
  s(i) = y(i);
end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
function [w,successList] = dmc_station_meta(w,varargin)

%DMC_STATION_META adds station metadata to a waveform
%  W = DMC_STATION_META(W) looks up station metadata at the IRIS Data
%  Management Center (DMC) based on the network, station and time of each
%  waveform. These fields are added to the waveform. W may be of aritrary
%  size and dimension. If a particular network-station-time combination is
%  not found at the DMC, that particular waveform is returned as is with no
%  additional fields.
%    STATION_STARTDATE        On date of the station (Matlab date format)
%    STATION_ENDDATE          Off date of the station (Matlab date format)
%    STATION_LATITUDE         latitude in degrees (double)
%    STATION_LONGITUDE        longitude in degress (double)
%    STATION_ELEVATION        elevation above sea level in meters (double)
%    STATION_NAME             full station name (text string)
%    STATION_NUMBEROFCHANNELS no. channels from station (double)
%
%  [W,SUCCESS] = DMC_STATION_META(W) returns a matrix SUCCESS which is the
%  same size as W. SUCCESS is a mask containing ones and zeros that can be
%  used to identify which waveforms where successfully tagged with metadata
%  from the DMC. Reasons that the metadata lookup might fail include:
%       - this network-station is not on record at the DMC
%       - this network-station is not on record at the time range requested
%       - there was no internet connection at the time of connection
%       - the DMC's web services were down (rare, but does happen)
%
%  [W,SUCCESS] = DMC_STATION_META(W,'CheckEach') For speed, this function
%  only looks up a given network-station combination once and applies it
%  to all elements of W that match the network-station. This is reasonable
%  for most W that might contain many waveforms from the same station and
%  general time range (e.g. a swarm of earthquakes over a few months). In
%  cases where the user suspects station data may have changed (e.g. a
%  decade of earthquakes from an unfamiliar network), the CheckEach flag
%  rerequests metadata for each individual waveform

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

% SETTINGS
httpServer = 'http://www.iris.edu/ws/station/';

% CHECK ARGUMENTS
if numel(varargin)>1
   error('dmc_station_meta:incorrectNumberArguments','Incorrect number of arguments');
end
if ~strcmpi(class(w),'waveform')
    error('dmc_station_meta:missingWaveform','First input must be a waveform object');
end
checkEach = 0;
if numel(varargin)>=1
    if strcmpi(varargin{1},'checkEach')
        checkEach = 1;
    end
end

% GET LIST OF UNIQUE CALLS TO SERVER
indexList = {[]};
if ~checkEach
    scnl = get(w,'SCNLOBJECT');
    netSta = cellstr(strcat([char(get(scnl,'NETWORK')) char(get(scnl,'STATION')) ]));
    uniqueSta = unique(netSta);
    for n = 1:length(uniqueSta)
        indexList(n) = {find(strcmpi(netSta,uniqueSta(n)))};
    end
else
    for n =1:length(w)
        indexList(n) = {n};
    end
end

% LOOP THROUGH THE SERVER CALLS
successList = zeros(size(w));
for n = 1:numel(indexList)
    index = indexList{n};
    scnlRequest = get(w(index(1)),'SCNLOBJECT');
    station = upper(get(scnlRequest,'STATION'));
    network = upper(get(scnlRequest,'NETWORK'));
    timeStart =  datestr(floor(get(w(index(1)),'START')),'yyyy-mm-dd');
    timeEnd   =  datestr(ceil(get(w(index(1)),'END')),'yyyy-mm-dd');
    httpLink = [httpServer 'query?net=' network '&sta=' station  '&timewindow=' timeStart ',' timeEnd '&level=sta'];
    %disp(httpLink);
    [s,success] = dmc_xml2struct(httpLink);
    if success
        w(index) = addfield(w(index),'STATION_STARTDATE',datenum(s.StaMessage.Station.StationEpoch.StartDate.Text,'yyyy-mm-ddTHH:MM:SS'));
        w(index) = addfield(w(index),'STATION_ENDDATE',datenum(s.StaMessage.Station.StationEpoch.EndDate.Text,'yyyy-mm-ddTHH:MM:SS'));
        w(index) = addfield(w(index),'STATION_LATITUDE',str2double(s.StaMessage.Station.StationEpoch.Lat.Text));
        w(index) = addfield(w(index),'STATION_LONGITUDE',str2double(s.StaMessage.Station.StationEpoch.Lon.Text));
        w(index) = addfield(w(index),'STATION_ELEVATION',str2double(s.StaMessage.Station.StationEpoch.Elevation.Text));
        w(index) = addfield(w(index),'STATION_NAME',char(s.StaMessage.Station.StationEpoch.Site.Country.Text));
        w(index) = addfield(w(index),'STATION_NUMBEROFCHANNELS',str2double(s.StaMessage.Station.StationEpoch.NumberChannels.Text));
        successList(index) = 1;
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function DirT = transformRTdir(Dir,Rt)

DirT = Rt(1:3,1:3) * Dir;


text:
function mv = mv_plotVoxData(mv,dims);
%
%  mv = mv_plotVoxData(mv,dims);
%
% MultiVoxel UI:
% Plot voxel data along 2-3 dimensions,
% averaging across the others.
%
% ras, 04/05.
if ieNotDefined('mv')
    mv = get(gcf,'UserData');
end

if ieNotDefined('dims')
    dims = [2 3];
end

dimNames = {'Trial' 'Voxel' 'Condition'};
selConds = find(tc_selectedConds(mv));
condNames = mv.trials.condNames(selConds); % ignore null
nConds = length(condNames);

if length(dims) < 3
    nSubplots = 1;
else
    nSubplots = size(mv.voxAmps,dims(3));
end

if nSubplots==1
    fontsz = mv.params.fontsz;
else
    fontsz = mv.params.fontsz - 3;
end

% get data to plot; take mean across
% any non-specified dimensions
data = mv.voxData(:,:,:,selConds-1);
avgAcrossDims = setdiff(1:4,dims);
if ~isempty(avgAcrossDims)
    for i = 1:length(avgAcrossDims)
        data = nanmeanDims(data,avgAcrossDims(i));
    end
end

% permute to plotting order
data = permute(data,[dims avgAcrossDims]);

% normalize to fit in the selected color map
cmap = mv.params.cmap;
data = double(normalize(data,1,size(cmap,1)));

% delete existing axes
other = findobj('Type','axes','Parent',gcf);
delete(other);

nrows = ceil(sqrt(nSubplots));
ncols = ceil(nSubplots/nrows);

for z = 1:nSubplots
    subplot(nrows,ncols,z);
    image(data(:,:,z));
    colormap(cmap);

    % check for conditions labeling:
    % we'll treat that specially, labeling
    % each condition separately
    if dims(1)==3
        set(gca,'YTick',1:nConds,'YTickLabel',condNames);
    else
        ylabel(dimNames{dims(1)},'FontName',mv.params.font,...
                'FontSize',fontsz);
    end

    if length(dims)>=2 & dims(2)==3
        set(gca,'XTick',1:nConds,'XTickLabel',condNames);
    else
        xlabel(dimNames{dims(2)},'FontName',mv.params.font,...
                'FontSize',fontsz);
    end

    if length(dims)==3 & dims(3)==3
        title(condNames{z},'FontName',mv.params.font,...
            'FontSize',fontsz+2);
    end
end

return


text:
function f = harmean(g,m,n)
% Implemets a harmonic mean filter.
inclass = class(g);
g = im2double(g);
f = m * n ./ imfilter(1./(g + eps),ones(m,n),'replicate');
f = changeclass(inclass,f);


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   q = solve_spherical_wrist(robot, q, T, wrist)
%   Solves the inverse kinematic problem for a spherical wrist
%   robot: robot structure.
%   q: vector containing the values of the joints 1, 2 and 3.
%   T: orientation of the last reference system.
%   wrist: select -1 or 1 for two possible solutions (wrist up, wrist down)
%
%	See also DIRECTKINEMATIC.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.

function q = solve_spherical_wrist(robot, q, T, wrist, method)

switch method

     %algebraic solution
    case 'algebraic'
        T01=dh(robot, q, 1);
        T12=dh(robot, q, 2);
        T23=dh(robot, q, 3);

        Q=inv(T23)*inv(T12)*inv(T01)*T;

        %detect the degenerate case when q(5)=0, this leads to zeros
        % in Q13, Q23, Q31 and Q32 and Q33=1
        thresh=1e-12;
        %detect if q(5)==0
        % this happens when cos(q5) in the matrix Q is close to 1
        if abs(Q(3,3)-1)>thresh
            %normal solution
            if wrist==1 %wrist up
                q(4)=atan2(-Q(2,3),-Q(1,3));
                q(6)=atan2(-Q(3,2),Q(3,1));
                %q(5)=atan2(-Q(3,2)/sin(q(6)),Q(3,3));
            else %wrist down
                q(4)=atan2(-Q(2,3),-Q(1,3))-pi;
                q(6)=atan2(-Q(3,2),Q(3,1))+pi;
                %q(5)=atan2(-Q(3,2)/sin(q(6)),Q(3,3));
            end
            if abs(cos(q(6)+q(4)))>thresh
                cq5=(Q(1,1)+Q(2,2))/cos(q(4)+q(6))-1;
            end
            if abs(sin(q(6)+q(4)))>thresh
                cq5=(-Q(1,2)+Q(2,1))/sin(q(4)+q(6))-1;
            end
            if abs(sin(q(6)))>thresh
                sq5=-Q(3,2)/sin(q(6));
            end
            if abs(cos(q(6)))>thresh
                sq5=Q(3,1)/cos(q(6));
            end
            q(5)=atan2(sq5,cq5);

        else %degenerate solution, in this case, q4 cannot be determined,
             % so q(4)=0 is assigned
            if wrist==1 %wrist up
                q(4)=0;
                q(5)=0;
                q(6)=atan2(-Q(1,2)+Q(2,1),Q(1,1)+Q(2,2));
            else %wrist down
                q(4)=-pi;
                q(5)=0;
                q(6)=atan2(-Q(1,2)+Q(2,1),Q(1,1)+Q(2,2))+pi;
            end

        end

       %geometric solution
    case 'geometric'
        % T is the noa matrix defining the position/orientation of the end
        % effector's reference system
        vx6=T(1:3,1);
        vz5=T(1:3,3); % The vector a z6=T(1:3,3) is coincident with z5

        % Obtain the position and orientation of the system 3
        % using the already computed joints q1, q2 and q3
        T01=dh(robot, q, 1);
        T12=dh(robot, q, 2);
        T23=dh(robot, q, 3);
        T03=T01*T12*T23;

        vx3=T03(1:3,1);
        vy3=T03(1:3,2);
        vz3=T03(1:3,3);

        % find z4 normal to the plane formed by z3 and a
        vz4=cross(vz3, vz5);	% end effector's vector a: T(1:3,3)

        % in case of degenerate solution,
        % when vz3 and vz6 are parallel--> then z4=0 0 0, choose q(4)=0 as solution
        if norm(vz4) <= 0.000001
            if wrist == 1 %wrist up
                q(4)=0;
            else
                q(4)=-pi; %wrist down
            end
        else
            %this is the normal and most frequent solution
            cosq4=wrist*dot(-vy3,vz4);
            sinq4=wrist*dot(vx3,vz4);
            q(4)=atan2(sinq4, cosq4);
        end
        %propagate the value of q(4) to compute the system 4
        T34=dh(robot, q, 4);
        T04=T03*T34;
        vx4=T04(1:3,1);
        vy4=T04(1:3,2);

        % solve for q5
        cosq5=dot(vy4,vz5);
        sinq5=dot(-vx4,vz5);
        q(5)=atan2(sinq5, cosq5);

        %propagate now q(5) to compute T05
        T45=dh(robot, q, 5);
        T05=T04*T45;
        vx5=T05(1:3,1);
        vy5=T05(1:3,2);

        % solve for q6
        cosq6=dot(vx6,vx5);
        sinq6=dot(vx6,vy5);
        q(6)=atan2(sinq6, cosq6);

    otherwise
        disp('no method specified in solve_spherical_wrist');
end

repo: ARTE


text:
function cs_star = surfaceConcentration(cs_barrato,jflux,Q,T,param)
% surfaceConcentration evaluates the concentration of Li-ions at the electrode surfaces.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

% Diffusion coefficients for the solid phase
[Dps_eff, Dns_eff] = param.SolidDiffusionCoefficientsFunction(T,param);

% Check what kind of solid diffusion model has been chosen.
if(param.SolidPhaseDiffusion==1) % Two parameters model
    % Evaluates the average surface concentration in both the electrodes.
    % Cathode
    cs_star_p = cs_barrato(1:param.Np)-(param.Rp_p./(Dps_eff.*5)).*jflux(1:param.Np);
    % Anode
    cs_star_n = cs_barrato(param.Np+1:end)-(param.Rp_n./(Dns_eff.*5)).*jflux(param.Np+1:end);
elseif(param.SolidPhaseDiffusion==2) % Three parameters model
    % Cathode
    cs_star_p = cs_barrato(1:param.Np)+(param.Rp_p./(Dps_eff.*35)).*(-jflux(1:param.Np)+8*Dps_eff.*Q(1:param.Np));
    % Anode
    cs_star_n = cs_barrato(param.Np+1:end)+(param.Rp_n./(Dns_eff.*35)).*(-jflux(param.Np+1:end)+8*Dns_eff.*Q(param.Np+1:end));
elseif(param.SolidPhaseDiffusion==3) % Full model
    p_indices = param.Nr_p:param.Nr_p:param.Nr_p*param.Np;
    n_indices = param.Nr_n:param.Nr_n:param.Nr_n*param.Nn;
    % If the full model has been used, just take the concentration data at r=Rp
    cs_star_p = cs_barrato(p_indices);
    cs_star_n = cs_barrato(n_indices+p_indices(end));
end
% Return the residuals
cs_star = [cs_star_p;cs_star_n];

end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
function camd_make
%CAMD_MAKE to compile camd for use in MATLAB
%
% Example:
%   camd_make
%
% See also camd.

% Copyright 1994-2007, Tim Davis, University of Florida,
% Patrick R. Amestoy, Iain S. Duff, and Yanqing Chen.

details = 0 ;	    % 1 if details of each command are to be printed

d = '' ;
if (~isempty (strfind (computer, '64')))
    d = '-largeArrayDims' ;
end

i = sprintf ('-I..%sInclude -I..%s..%sUFconfig', filesep, filesep, filesep) ;
cmd = sprintf ('mex -O %s -DDLONG -output camd %s camd_mex.c', d, i) ;
files = {'camd_order', 'camd_dump', 'camd_postorder', ...
    'camd_aat', 'camd_2', 'camd_1', 'camd_defaults', 'camd_control', ...
    'camd_info', 'camd_valid', 'camd_global', 'camd_preprocess' } ;
for i = 1 : length (files)
    cmd = sprintf ('%s ..%sSource%s%s.c', cmd, filesep, filesep, files {i}) ;
end
if (details)
    fprintf ('%s\n', cmd) ;
end
eval (cmd) ;

fprintf ('CAMD successfully compiled.\n') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [success] = saveNeurodataWithoutBorders(image_masks,roi_response_data,algorithm,outputFilePath,varargin)
	% Takes cell extraction outputs and saves then in NWB format per format in https://github.com/schnitzer-lab/nwb_schnitzer_lab.
	% Biafra Ahanonu
	% started: 2020.04.03 [16:01:45]
	% Based on mat2nwb in https://github.com/schnitzer-lab/nwb_schnitzer_lab.
	% inputs
		% image_masks - [x y z] matrix
		% roi_response_data - {1 N} cell with N = number of different signal traces for that algorithm. Make sure each signal trace matrix is in form of [nSignals nFrames].
		% algorithm - Name of the algorithm.
		% outputFilePath - file path to save NWB file to.
	% outputs
		%

	% changelog
		% 2020.07.01 [09:40:20] - Convert roi_response_data to cell if user inputs only a matrix.
		% 2020.09.15 [20:30:32] - Automatically creates directory where file is to be stored if it is not present.
		% 2021.02.01 [?15:14:40] - Function checks that yaml, matnwb, and nwb_schnitzer_lab loaded, else tries to load to make sure all dependencies are present and active.
		% 2021.02.01 [15:19:40] - Update `_external_programs` to call ciapkg.getDirExternalPrograms() to standardize call across all functions.
		% 2021.02.03 [12:34:06] - Added a check for inputs with a single signal and function returns as it is not supported.
		% 2021.03.20 [19:35:28] - Update to checking if only a single signal input.
	% TODO
		%

	[success] = ciapkg.nwb.saveNeurodataWithoutBorders(image_masks,roi_response_data,algorithm,outputFilePath,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function []=panel2disp(n,N,m,p,k,T,Ymat,Xmat,Units,endo,exo,const,beta_gibbs,B_median,beta_median,beta_std,beta_lbound,beta_ubound,sigma_gibbs,sigma_median,D_estimates,gamma_estimates,ar,lambda1,lambda3,lambda4,startdate,enddate,forecast_record,forecast_estimates,Fcperiods,stringdates3,Fstartdate,Fcenddate,Feval,Fcomp,data_endo_c,data_endo_c_lags,data_exo_c,It,Bu,IRF,IRFt,pref,names,PriorExcel)

% obtain first a point estimate betatilde of the VAR coefficients
% this is simply the median
betatilde=beta_median;
Btilde=reshape(betatilde,k,n);

% check whether the model is stationary
[stationary eigmodulus]=bear.checkstable(betatilde,n,p,k);

% the other measures are unit-specific, hence, loop over units
for ii=1:N

% obtain predicted values
Yp(:,:,ii)=Xmat(:,:,ii)*B_median;
% then produce the corresponding residuals, using (1.9.4)
EPS(:,:,ii)=Ymat(:,:,ii)-Yp(:,:,ii);

% Compute then the sum of squared residuals
% compute first the RSS matrix, defined in (1.9.5)
RSS(:,:,ii)=EPS(:,:,ii)'*EPS(:,:,ii);
% retain only the diagonal elements to get the vector of RSSi values
rss(:,:,ii)=diag(RSS(:,:,ii));

% Go on calculating R2
% generate Mbar
Mbar=eye(T)-ones(T,T)/T;
% then compute the TSS matrix, defined in (1.9.8)
TSS(:,:,ii)=Ymat(:,:,ii)'*Mbar*Ymat(:,:,ii);
% generate the R2 matrix in (1.9.9)
R2(:,:,ii)=eye(n)-RSS(:,:,ii)./TSS(:,:,ii);
% retain only the diagonal elements to get the vector of R2 values
r2(:,:,ii)=diag(R2(:,:,ii));

% then calculate the adjusted R2, using (1.9.11)
R2bar(:,:,ii)=eye(n)-((T-1)/(T-k))*(eye(n)-R2(:,:,ii));
% retain only the diagonal elements to get the vector of R2bar values
r2bar(:,:,ii)=diag(R2bar(:,:,ii));

end

% now start displaying and saving the results

% preliminary task: create and open the txt file used to save the results

filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'wt');

% print toolbox header

fprintf('%s\n','');
fprintf(fid,'%s\n','');

% print the list of contributors
bear.printcontributors(fid);

% print then estimation results

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

toolboxinfo='BEAR toolbox estimates';
fprintf('%s\n',toolboxinfo);
fprintf(fid,'%s\n',toolboxinfo);

time=clock;
datestring=datestr(time);
dateinfo=['Date: ' datestring(1,1:11) '   Time: ' datestring(1,13:17)];
fprintf('%s\n',dateinfo);
fprintf(fid,'%s\n',dateinfo);

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

VARtypeinfo='Panel VAR: pooled estimator';
fprintf('%s\n',VARtypeinfo);
fprintf(fid,'%s\n',VARtypeinfo);

if IRFt==1
SVARinfo='structural decomposition: none';
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);
elseif IRFt==2
SVARinfo='structural decomposition: choleski factorisation';
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);
elseif IRFt==3
SVARinfo='structural decomposition: triangular factorisation';
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);
end

temp='units: ';
for ii=1:N
temp=[temp ' ' Units{ii,1} ' '];
end
unitinfo=temp;
fprintf('%s\n',unitinfo);
fprintf(fid,'%s\n',unitinfo);

temp='endogenous variables: ';
for ii=1:n
temp=[temp ' ' endo{ii,1} ' '];
end
endoinfo=temp;
fprintf('%s\n',endoinfo);
fprintf(fid,'%s\n',endoinfo);

temp='exogenous variables: ';
if const==0 && m==0
temp=[temp ' none'];
elseif const==1 && m==1
temp=[temp ' constant '];
elseif const==0 && m>0
   for ii=1:m-1
   temp=[temp ' ' exo{ii,1} ' '];
   end
elseif const==1 && m>1
temp=[temp ' constant '];
   for ii=1:m-1
   temp=[temp ' ' exo{ii,1} ' '];
   end
end
exoinfo=temp;
fprintf('%s\n',exoinfo);
fprintf(fid,'%s\n',exoinfo);

sampledateinfo=['estimation sample: ' startdate '-' enddate];
fprintf('%s\n',sampledateinfo);
fprintf(fid,'%s\n',sampledateinfo);

samplelengthinfo=['sample size (omitting initial conditions): ' num2str(T)];
fprintf('%s\n',samplelengthinfo);
fprintf(fid,'%s\n',samplelengthinfo);

laginfo=['number of lags included in regression: ' num2str(p)];
fprintf('%s\n',laginfo);
fprintf(fid,'%s\n',laginfo);

hyperparam1='hyperparameters:';
fprintf('%s\n',hyperparam1);
fprintf(fid,'%s\n',hyperparam1);

if PriorExcel==1
    arprint=[];
    for ii=1:n
        arprint=[arprint num2str(ar(ii,1)) '  '];
    end
 hyperparam2=['autoregressive coefficients (ar):                ' arprint];
else
 hyperparam2=['autoregressive coefficients (ar):                ' num2str(ar(1,1))];
end
fprintf('%s\n',hyperparam2);
fprintf(fid,'%s\n',hyperparam2);

hyperparam3=['overall tightness (lambda1):              ' num2str(lambda1)];
fprintf('%s\n',hyperparam3);
fprintf(fid,'%s\n',hyperparam3);

hyperparam4=['lag decay (lambda3):                      ' num2str(lambda3)];
fprintf('%s\n',hyperparam4);
fprintf(fid,'%s\n',hyperparam4);

hyperparam5=['exogenous variable tightness (lambda4):   ' num2str(lambda4(1,1))];
fprintf('%s\n',hyperparam5);
fprintf(fid,'%s\n',hyperparam5);

% display coefficient estimates

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

coeffinfo=['VAR coefficients (Common to all units):'];
fprintf('%s\n',coeffinfo);
fprintf(fid,'%s\n',coeffinfo);

for ii=1:n
fprintf('%s\n','');
fprintf(fid,'%s\n','');
if ii~=1
fprintf('%s\n','');
fprintf(fid,'%s\n','');
end

endoinfo=['Endogenous: ' endo{ii,1}];
fprintf('%s\n',endoinfo);
fprintf(fid,'%s\n',endoinfo);

coeffheader=fprintf('%25s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');
coeffheader=fprintf(fid,'%25s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');

% handle the endogenous
   for jj=1:n
      for kk=1:p
      values=[beta_median((ii-1)*k+n*(kk-1)+jj,1) beta_std((ii-1)*k+n*(kk-1)+jj,1) beta_lbound((ii-1)*k+n*(kk-1)+jj,1) beta_ubound((ii-1)*k+n*(kk-1)+jj,1)];
      fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1},'(-',int2str(kk),')'),values);
      fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1},'(-',int2str(kk),')'),values);
      end
   end

% handle the exogenous
   % if there is no constant:
   if const==0
      % if there is no exogenous at all, obvioulsy, don't display anything
      if m==0
      % if there is no constant but some other exogenous, display them
      else
         for jj=1:m
         values=[beta_median(ii*k-m+jj,1) beta_std(ii*k-m+jj,1) beta_lbound(ii*k-m+jj,1) beta_ubound(ii*k-m+jj,1)];
         fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
         fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
         end
      end
   % if there is a constant
   else
   % display the results related to the constant
         values=[beta_median(ii*k-m+1,1) beta_std(ii*k-m+1,1) beta_lbound(ii*k-m+1,1) beta_ubound(ii*k-m+1,1)];
         fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n','Constant',values);
         fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n','Constant',values);
      % if there is no other exogenous, stop here
      if m==1
      % if there are other exogenous, display their results
      else
         for jj=1:m-1
         values=[beta_median(ii*k-m+jj+1,1) beta_std(ii*k-m+jj+1,1) beta_lbound(ii*k-m+jj+1,1) beta_ubound(ii*k-m+jj+1,1)];
         fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
         fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
         end
      end
   end

fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display evaluation measures

   % loop over units
   for jj=1:N
   unitinfo=['unit: ' Units{jj,1}];
   fprintf('%s\n',unitinfo);
   fprintf(fid,'%s\n',unitinfo);

   rssinfo=['Sum of squared residuals: ' num2str(rss(ii,1,jj),'%.2f')];
   fprintf('%s\n',rssinfo);
   fprintf(fid,'%s\n',rssinfo);

   r2info=['R-squared: ' num2str(r2(ii,1,jj),'%.3f')];
   fprintf('%s\n',r2info);
   fprintf(fid,'%s\n',r2info);

   adjr2info=['adj. R-squared: ' num2str(r2bar(ii,1,jj),'%.3f')];
   fprintf('%s\n',adjr2info);
   fprintf(fid,'%s\n',adjr2info);

   fprintf('%s\n','');
   fprintf(fid,'%s\n','');

   end

end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display VAR stability results
eigmodulus=reshape(eigmodulus,p,n);
stabilityinfo1=['Roots of the characteristic polynomial (modulus):'];
fprintf('%s\n',stabilityinfo1);
fprintf(fid,'%s\n',stabilityinfo1);
for jj=1:p
temp=num2str(eigmodulus(jj,1),'%.3f');
   for kk=2:n
   temp=[temp,'  ',num2str(eigmodulus(jj,kk),'%.3f')];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end
if stationary==1;
stabilityinfo2=['No root lies outside the unit circle.'];
stabilityinfo3=['The estimated VAR model satisfies the stability condition'];
fprintf('%s\n',stabilityinfo2);
fprintf(fid,'%s\n',stabilityinfo2);
fprintf('%s\n',stabilityinfo3);
fprintf(fid,'%s\n',stabilityinfo3);
else
stabilityinfo2=['Warning: at leat one root lies on or outside the unit circle.'];
stabilityinfo3=['The estimated VAR model will not be stable'];
end
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display posterior for sigma
sigmainfo=['sigma (residual covariance matrix): posterior estimates'];
fprintf('%s\n',sigmainfo);
fprintf(fid,'%s\n',sigmainfo);
% calculate the (integer) length of the largest number in sigma, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(sigma_median))))));
% add a separator, a potential minus sign, and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:n
temp=[];
   for jj=1:n
   % convert matrix entry into string
   number=num2str(sigma_median(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

% then display the results for D and gamma, if a structural decomposition was selected

if IRF==1 && IRFt~=1
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

svarinfo1=['D (structural decomposition matrix): posterior estimates'];
fprintf('%s\n',svarinfo1);
fprintf(fid,'%s\n',svarinfo1);

% recover D
D=reshape(D_estimates,n,n);
% calculate the (integer) length of the largest number in D, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(D))))));
% add a separator, a potential minus sign and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:n
temp=[];
   for jj=1:n
   % convert matrix entry into string
   number=num2str(D(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

svarinfo2=['gamma (structural disturbances covariance matrix): posterior estimates'];
fprintf('%s\n',svarinfo2);
fprintf(fid,'%s\n',svarinfo2);

% recover gamma
gamma=reshape(gamma_estimates,n,n);
% calculate the (integer) length of the largest number in D, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(gamma))))));
% add a separator, a potential minus sign and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:n
temp=[];
   for jj=1:n
   % convert matrix entry into string
   number=num2str(gamma(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% finally, estimate and print the unit-specific elements: for this model, only forecast evaluation remains to be printed
% first note that forecast evaluation can only be conducted if it is activated, and if there is some observable data after the beginning of the forecast
if Feval==1 && Fcomp==1

   % loop over units
   for ii=1:N
   unitinfo=['unit-specific components: ' Units{ii,1}];
   fprintf('%s\n',unitinfo);
   fprintf(fid,'%s\n',unitinfo);
   fprintf('%s\n','');
   fprintf(fid,'%s\n','');
   % compute forecast evaluation
   [RMSE MAE MAPE Ustat CRPS_estimates S1_estimates S2_estimates]=bear.panelfeval(n,p,k,beta_gibbs,sigma_gibbs,forecast_record(:,:,ii),forecast_estimates(:,:,ii),Fcperiods,data_endo_c(:,:,ii),data_endo_c_lags(:,:,ii),data_exo_c,const,It,Bu);
   % then display and save the results
   bear.panelfprint(n,endo,RMSE,MAE,MAPE,Ustat,CRPS_estimates,S1_estimates,S2_estimates,stringdates3,Fstartdate,Fcenddate,Fcperiods,fid);
   fprintf('%s\n','');
   fprintf(fid,'%s\n','');
   fprintf('%s\n','');
   fprintf(fid,'%s\n','');
   fprintf('%s\n','');
   fprintf(fid,'%s\n','');
   fprintf('%s\n','');
   fprintf(fid,'%s\n','');
   end

% if forecast evaluation is activated but not possible, return a message to signal it
elseif Feval==1 && Fcomp==0

finfo1=['Forecast evaluation cannot be conducted.'];
fprintf('%s\n',finfo1);
fprintf(fid,'%s\n',finfo1);
finfo2=['Forecasts start in ' Fstartdate ', while observable data is available only until ' names{end,1} '.'];
fprintf('%s\n',finfo2);
fprintf(fid,'%s\n',finfo2);
finfo3=['To obtain forecast evaluation, the forecast start date must be anterior to the end of the data set.'];
fprintf('%s\n',finfo3);
fprintf(fid,'%s\n',finfo3);

% if forecast evaluation is not activated altogether, do not do anything
end

fclose(fid);

author: european-central-bank

repo: BEAR-toolbox


text:
% function H = histo2D(D,[Xlo Xhi],Xn,[Ylo Yhi],Yn,Xlab,Ylab,Title)
%
% 2 Dimensional Histogram (size(H) == [Yn Xn])
% Counts number of points in the bins defined by
% X = linspace(Xlo,Xhi,Xn) and
% Y = linspace(Ylo,Yhi,Yn)

function H = histo2D(D,Xrange,Xn,Yrange,Yn,Xlab,Ylab,Title)

Xlo = Xrange(1) ; Xhi = Xrange(2) ;
Ylo = Yrange(1) ; Yhi = Yrange(2) ;
X = linspace(Xlo,Xhi,Xn)' ;
Y = linspace(Ylo,Yhi,Yn)' ;

Dx = D(:,1) ; Dy = D(:,2) ;
n = length(D) ;

H = zeros(Yn,Xn) ;

for i = 1:n
    x = dsearchn(X,Dx(i)) ;
    y = dsearchn(Y,Dy(i)) ;
    H(y,x) = H(y,x) + 1 ;
end ;

figure , pcolor(X,Y,H) ;
% Xmid = 0.5*(X(1:end-1)+X(2:end)) ;
% Ymid = 0.5*(Y(1:end-1)+Y(2:end)) ;
% figure , pcolor(Xmid,Ymid,H) ;
colorbar ; shading flat ; axis square tight ; grid on ;
xlabel(Xlab) ; xlabel(Ylab) ; title(Title) ;


text:
% The COBRAToolbox: testTransformModel2KEGG.m
%
% Purpose:
%     - test the transformModel2KEGG function
%
% Authors:
%     - Jacek Wachowiak
global CBTDIR
% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testTransformModel2KEGG'));
cd(fileDir);

% test variables
model = getDistributedModel('ecoli_core_model.mat');
Dictionary = {'0', '1'; '2', '3'};
ref_KEGGID = cell(72, 1);

% function outputs
modelKEGG = transformModel2KEGG(model, Dictionary);
try
    modelKEGG_2 = transformModel2KEGG(model);
catch ME
    assert(length(ME.message) > 0);
end

% test
assert(isequal(modelKEGG.metsAbr, model.mets));
assert(isequal(modelKEGG.metKEGGID, ref_KEGGID));

% change to old directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
function [s] = getsubfield(s, f);

% GETSUBFIELD returns a field from a structure just like the standard
% Matlab GETFIELD function, except that you can also specify nested fields
% using a '.' in the fieldname. The nesting can be arbitrary deep.
%
% Use as
%   f = getsubfield(s, 'fieldname')
% or as
%   f = getsubfield(s, 'fieldname.subfieldname')
%
% See also GETFIELD, ISSUBFIELD, SETSUBFIELD

% Copyright (C) 2005, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: getsubfield.m 2885 2011-02-16 09:41:58Z roboos $

if ~ischar(f)
  error('incorrect input argument for fieldname');
end

t = {};
while (1)
  [t{end+1}, f] = strtok(f, '.');
  if isempty(f)
    break
  end
end
s = getfield(s, t{:});

author: PatternRecognition

repo: OpenBMI


text:
%% OPT_BATCH_ITHACA uses the BATCH command to run the OPT code on Ithaca.
%
%  Discussion:
%
%    The OPT code is a function, so first we must write a script
%    called OPT_SCRIPT that runs the function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2010
%
%  Author:
%
%    John Burkardt
%
  clear

  fprintf ( 1, '\n' );
  fprintf ( 1, 'OPT_BATCH_ITHACA\n' );
  fprintf ( 1, '  Run OPT_SCRIPT locally.\n' );
%
%  BATCH sends the script for execution.
%
  my_job = batch ( 'opt_script', ...
    'CaptureDiary', true, ...
    'FileDependencies', { 'quad_fun' }, ...
    'Configuration', 'ithaca_2009b', ...
    'CurrentDirectory', '/home/burkardt/matlab', ...
    'matlabpool', 4 );
%
%  WAIT pauses the MATLAB session til the job completes.
%
  wait ( my_job );
%
%  DIARY displays any messages printed during execution.
%
  diary ( my_job );
%
%  LOAD makes the script's workspace available:
%
%  z_star = optimal control.
%  PAR = problem parameters.
%  exit = positive if calculation was successful.
%
  load ( my_job );
%
%  On successful conclusion, display the plot.
%
  if ( 0 <= exit && isfield ( PAR, 'plot' ) )
    feval ( PAR.plot, z_star, PAR )
  end

%
%  DESTROY cleans up data about the job we no longer need.
%
  destroy ( my_job );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'OPT_BATCH_ITHACA\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

author: johannesgerer

repo: jburkardt-m


text:
%% load data

load('FisherIris.mat')
rng('default');

%% Run umap

umap_coords = run_umap(meas);

% see also 'graphic' cluster output

%%
f1 = create_figure('umap on iris dataset', 1, 3);

plot(umap_coords(:, 1), umap_coords(:, 2), 'ko');

xlabel('umap(1)'); ylabel('umap(2)');
title('UMAP with colors indicating true classes');

[indic, names, condf] = string2indicator(species);
colors = scn_standard_colors(length(names));

for i = 1:length(names)

    wh = condf == i;

    plot(umap_coords(wh, 1), umap_coords(wh, 2), 'ko', 'MarkerFaceColor', colors{i});

end

drawnow

%% Cluster data in derived umap space and re-plot

% note: clusterIdentifiers from run_umap didn't produce great results for me (tor)

% note: you really want enough (> 2) umap dimensions for this to be meaningful.
% Run again with 3 umap dims
% Will run with same number of dims as original dimensions, but not sure
% yet if this is a good idea.

[umap_coords3d, umap, umap_clusterIdentifiers, extras] = run_umap(meas, 'cluster_output', 'numeric', 'n_components', 3);

clusterIdentifiers = clusterdata(umap_coords3d, 'linkage','ward','savememory','on','maxclust', 4);

%% re-plot

figure(f1)
subplot(1, 3, 2)

plot(umap_coords(:, 1), umap_coords(:, 2), 'ko');

xlabel('umap(1)'); ylabel('umap(2)');
title('UMAP with colors indicating estimated clusters');

n = length(unique(clusterIdentifiers));

colors = scn_standard_colors(n);

for i = 1:n

    wh = clusterIdentifiers == i;

    plot(umap_coords(wh, 1), umap_coords(wh, 2), 'ko', 'MarkerFaceColor', colors{i});

end

drawnow

%%
figure(f1)
subplot(1, 3, 3)

plot(umap_coords3d(:, 1), umap_coords3d(:, 2), 'ko');

xlabel('umap(1)of 3'); ylabel('umap(2) of 3');
title('UMAP with colors indicating estimated clusters');

n = length(unique(clusterIdentifiers));

colors = scn_standard_colors(n);

for i = 1:n

    wh = clusterIdentifiers == i;

    plot(umap_coords3d(wh, 1), umap_coords3d(wh, 2), 'ko', 'MarkerFaceColor', colors{i});

end

drawnow

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function plot_miss_rate_min(plot_obj,plot_args,legend_string)
% Plots the contribution of the misses to the minimum dcf curve.
% Inputs:
%   plot_args: A cell array of arguments to be passed to 'plot' that control
%     the appearance of the curve. See Matlab's help on 'plot' for information.
%   legend_string: Optional. A string to describe this curve in the legend.
%     This string will be prepended with the system name.

if ~exist('legend_string','var')
    legend_string = '';
end

plot_obj.plot_miss_rate(plot_obj.minPmiss,plot_args,legend_string);
end

author: nesl


text:
classdef MimToolCallback < handle
    % MimToolCallback.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        Axes
        Axes3d
        ImageDisplayParameters
        Reporting
        ViewerPanel
    end

    methods

        function obj = MimToolCallback(viewing_panel, image_display_parameters, axes, axes_3d, reporting)
            obj.ImageDisplayParameters = image_display_parameters;
            obj.ViewerPanel = viewing_panel;
            obj.Axes = axes;
            obj.Axes3d = axes_3d;
            obj.Reporting = reporting;
        end

        function EnablePan(obj, enabled)
            obj.GetAxes.EnablePan(enabled);
            obj.GetAxes3d.EnablePan(enabled);
        end

        function EnableZoom(obj, enabled)
            obj.GetAxes.EnableZoom(enabled);
            obj.GetAxes3d.EnableZoom(enabled);
        end

        function EnableRotate3d(obj, enabled)
            obj.GetAxes3d.EnableRotate3d(enabled);
        end

        function [min_coords, max_coords] = GetImageLimits(obj)
            % Gets the current limits of the visible image axes

            [min_coords, max_coords] = obj.GetAxes.GetImageLimits;
        end

        function SetImageLimits(obj, min_coords, max_coords)
            % Adjusts the image axes to make the image visible between the specified
            % coordinates

            x_lim = [min_coords(1), max_coords(1)];
            y_lim = [min_coords(2), max_coords(2)];
            obj.GetAxes.SetLimits(x_lim, y_lim);
        end

        function SetWindowWithinLimits(obj, window)
            % Sets the window subject to the current constraints

            window_limits = obj.ViewerPanel.WindowLimits;
            if ~isempty(window_limits)
                window = max(window, window_limits(1));
                window = min(window, window_limits(2));
                obj.ImageDisplayParameters.Window = window;
            end
        end

        function SetLevelWithinLimits(obj, level)
            % Sets the level subject to the current constraints

            level_limits = obj.ViewerPanel.LevelLimits;
            if ~isempty(level_limits)
                level = max(level, level_limits(1));
                level = min(level, level_limits(2));
                obj.ImageDisplayParameters.Level = level;
            end
        end

        function axes_handle = GetAxes(obj)
            axes_handle = obj.Axes;
            if isempty(axes_handle)
                obj.Reporting.Error('MimToolCallback:AxesDoNotExist', 'Axes have not been created');
            end
        end

        function axes_handle = GetAxes3d(obj)
            axes_handle = obj.Axes3d.GetRenderAxes;
            if isempty(axes_handle)
                obj.Reporting.Error('MimToolCallback:AxesDoNotExist', '3D Axes have not been created');
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [y_in_x,x_in_y]=cosmo_overlap(xs,ys)
% compute overlap between vectors or cellstrings in two cells
%
% [y_in_x,x_in_y]=cosmo_overlap(xs,ys)
%
% Inputs:
%    xs         Nx1 cell with all elements either numeric arrays or
%               cells with strings
%    ys         Mx1 cell with all elements either numeric arrays or
%               cells with strings
%
% Output:
%    y_in_x     MxN cell with the (i,j)-th element indicating the ratio of
%               elements in ys{j} that are present in xs{i}
%    x_in_i     MxN cell with the (i,j)-th element indicating the ratio of
%               elements in ys{i} that are present in ys{i}
%
% Examples:
%     % Compute overlap between two cells with cellstrings
%     xs={{'a'},{'a','b'},{'a','b','c'},{}};
%     ys={{'b'},{'c','b','a'}};
%     [x_in_y,y_in_x]=cosmo_overlap(xs,ys);
%     cosmo_disp(x_in_y);
%     %|| [ 0     0.333
%     %||   1     0.667
%     %||   1         1
%     %||   0         0 ]
%     %||
%     cosmo_disp(y_in_x);
%     %|| [     0         1
%     %||     0.5         1
%     %||   0.333         1
%     %||     NaN       NaN ]
%
%     % Compute overlap between two cells with numeric arrays
%     xs={1,[1 2],1:3,[]};
%     ys={2,[3,2,1]};
%     [x_in_y,y_in_x]=cosmo_overlap(xs,ys);
%     cosmo_disp(x_in_y);
%     %|| [ 0     0.333
%     %||   1     0.667
%     %||   1         1
%     %||   0         0 ]
%     %||
%     cosmo_disp(y_in_x);
%     %|| [     0         1
%     %||     0.5         1
%     %||   0.333         1
%     %||     NaN       NaN ]
%
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    [xs_vec,nx,nxs,xi]=get_counts(xs);
    [ys_vec,ny,nys,yi]=get_counts(ys);

    xy_cell=[xs_vec(:);ys_vec(:)];
    xyc=cat(1,xy_cell{:});

    % index xs from 1 to nx, and ys from (nx+1) to (nx+ny)
    xyi=[xi;(nx+yi)];

    % space for histogram
    h=zeros(nx,ny);

    % position of index of last value in xs
    xs_pos_last=sum(nxs);

    % because string comparisons are slow, use their indices instead
    [unq,unused,idxs]=unique(xyc);
    nunq=numel(unq);
    [idxs_sorted,i_sorted]=sort(idxs);

    msk=[true; any(diff(idxs_sorted,1),2)];
    unq_start_end_pos=[find(msk); 1+numel(msk)];

    for j=1:nunq
        % more readible, but slower:
        %    start_pos=unq_start_end_pos(j);
        %    end_pos=unq_start_end_pos(j+1)-1;
        %    if i_sorted(start_pos)>xs_pos_last [...]
        if i_sorted(unq_start_end_pos(j))<=xs_pos_last && ...
                i_sorted(unq_start_end_pos(j+1)-1)>xs_pos_last

            start_pos=unq_start_end_pos(j);
            end_pos=unq_start_end_pos(j+1)-1;

            i=i_sorted(start_pos:(end_pos));

            first_y=find_first_greater_than(i,xs_pos_last);
            px=xyi(i(1:(first_y-1)));
            py=xyi(i(first_y:end))-nx;

            h(px,py)=h(px,py)+1;
        end
    end

    x_in_y=bsxfun(@rdivide,h,nxs);
    y_in_x=bsxfun(@rdivide,h,nys');

function i=find_first_greater_than(sorted_vs,thr)
    % using binary search, find the first position i in sorted_vs
    % so that sorted(vs)>thr
    % it is assumed that sorted_vs is sorted
    if sorted_vs(end)<=thr
        i=numel(sorted_vs)+1;
        return
    end

    first=1;
    last=numel(sorted_vs);

    while first<last
        mid=floor((first+last)/2);
        if sorted_vs(mid)<=thr
            first=mid+1;
        else
            last=mid;
        end
    end
    i=first;
    %assert(sorted_vs(i)>thr);
    %assert(i==1 || sorted_vs(i-1)<=thr);

function [xs_vec,n,c,i]=get_counts(xs)
    n=numel(xs);
    c=zeros(n,1);
    xs_vec=cell(n,1);
    for k=1:n
        xsk=xs{k};
        if ~isnumeric(xsk) && ~iscellstr(xsk)
            error('only cells with numeric or cellstr input is supported');
        end
        c(k)=numel(xsk);
        xs_vec{k}=xsk(:);
    end

    nc=sum(c);
    i=zeros(nc,1);
    pos=0;
    for k=1:n
        ck=c(k);
        idxs=pos+(1:ck);
        i(idxs)=k;
        pos=pos+ck;
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
classdef prtOptionsClusterPlot
    % xxx Need Help xxx

    properties

        nSamplesPerDim = [500 100 20]; % Number of samples to use for plotting
        colorsFunction = @prtPlotUtilClassColors; % Colors function handle

        twoClassColorMapFunction = @prtPlotUtilTwoClassColorMap; % Two class colormap function handle
    end

    methods
        function obj = prtOptionsClusterPlot(varargin)
            obj = prtUtilAssignStringValuePairs(obj,varargin{:});
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function prob_test094 ( )

%*****************************************************************************80
%
%% TEST094 tests INVERSE_GAUSSIAN_MEAN, INVERSE_GAUSSIAN_SAMPLE, INVERSE_GAUSSIAN_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST094\n' );
  fprintf ( 1, '  For the Inverse Gaussian PDF:\n' );
  fprintf ( 1, '  INVERSE_GAUSSIAN_MEAN computes the mean;\n' );
  fprintf ( 1, '  INVERSE_GAUSSIAN_SAMPLE samples;\n' );
  fprintf ( 1, '  INVERSE_GAUSSIAN_VARIANCE computes the variance.\n' );

  a = 2.0;
  b = 3.0;

  check = inverse_gaussian_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST094 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = inverse_gaussian_mean ( a, b );
  variance = inverse_gaussian_variance ( a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF parameter B =             %14f\n', b );
  fprintf ( 1, '  PDF mean =                    %14f\n', mean );
  fprintf ( 1, '  PDF variance =                %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = inverse_gaussian_sample ( a, b, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end


text:
function color = getColorFromColorScale(value, colorScale)
%getColorFromColorScale Obtains color from color scale
%
% color = getColorFromColorScale(value, colorScale)
%
%INPUT
% value         Vector containing values
%
%OPTIONAL INPUT
% colorScale    Color scale used to dermine color for values
%               (Default = cool(100))
%
%OUTPUT
% color         n x 3 matrix containing RGB color values for each input
%               value
%
%Richard Que 12/2009
%
if ~exist('colorScale','var')
    colorScale = (cool(100));
    colorScale = round(colorScale*255);
end
color = zeros(size(value,1),3);
%normalize value
if max(value)~=0 & max(value)>1
    value = abs(value/max(value));
end
value(isinf(value))=0;
index = round(value*(size(colorScale,1)-1));
for i=1:size(value,1)
    color(i,1:3) = colorScale(index(i)+1,1:3);
end


text:
function varargout = subsref(this, S)
% Allows subscript referencing with dot notation of dimLabels
%
%   Y = MrDimInfo()
%   Y.subsref(S)
%
% This is a method of class MrDimInfo.
% It checks whether the provided S(1).subs is
%   a) a valid dimLabel and returns the reduced MrDimInfo for this dimension (via get_dims(dimLabel)).
%       e.g. dimInfo.z.samplingPoints
%   b) a valid property of MrDimInfo and S(2).subs a valid dimLabel, which
%       is transformed into a dim-index to allow referencing
%       e.g. dimInfo.nSamples('z')
% For anything but these first-level dot notations, it uses the builtin subsref.
%
% IN
%   S   structure with two fields:
%           type is a char vector containing '()', '{}', or '.', indicating the type of indexing used.
%           NOTE: we only overload the dot notation here to index dimLabels
%           subs is a cell array or character array containing the actual subscripts.% OUT
%
% EXAMPLE
%   This enables the usage of dimLabels in the dot notation, e.g.
%
%   dimInfo = MrDimInfo('dimLabels', {'x', 'y', 'z'},
%                       'samplingPoints', {1:10, -10:-1, 5:14});
%   dimInfo.z.samplingPoints
%       => ans =
%      5     6     7     8     9    10    11    12    13    14
%
%   dimInfo.nSamples('z');
%
%   See also MrDimInfo builtin.subsref

% Author:   Lars Kasper & Saskia Bollmann
% Created:  2017-06-28
% Copyright (C) 2017 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (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/>.

switch S(1).type
    case '.'
        % do custom dot-referencing  for valid dimLabels (e.g. dimInfo.x etc)
        if ismember(S(1).subs, this.dimLabels)

            if numel(S) == 1 % return reduced dimInfo object
                varargout = {this.get_dims(S(1).subs)};
            else
                % retrieve reduced dimInfo and continue with classical subsref
                % from there
                varargout = {builtin('subsref',this.get_dims(S(1).subs),S(2:end))};
            end
        elseif ismember(S(1).subs, properties(this)) && numel(S) > 1 && numel(S(2).subs) == 1
                % do custom dot-referencing allowing for property(dimLabel), e.g. resolutions('x')
                % by converting char/cell indices to numerical ones and run normal
                % subsref
                % note: allows only 1D indexing, i.e. not for ranges(:, 1:2)
                S(2).subs = {this.get_dim_index(S(2).subs{:})};
                varargout = {builtin('subsref',this,S)};
        elseif ismember(S(1).subs, methods(this))
            % a method call in dot notation, e.g. dimInfo.select() shall be
            % treated correctly, even for variable number of output
            % arguments
            if nargout
                [varargout{1:nargout}] = builtin('subsref',this,S);
            else
                varargout = {builtin('subsref',this,S)};
            end
        else % hope for the best
                varargout = {builtin('subsref',this,S)};
        end
    case '()' % allow to retrieve dimInfo('x'), dimInfo({'x','y'}), dimInfo(1:3)
        % for valid dimension labels/ indices
        sub = S(1).subs{:};
        isValidDimLabel = (ischar(sub) || iscellstr(sub)) && all(ismember(sub, this.dimLabels));
        % dimIndices have to be integer and within range of validdimensions
        isValidDimIndex = isnumeric(sub) && isequal(fix(sub),sub) ...
            && all(sub >=1) && all(sub <= this.nDims);
        if isValidDimLabel || isValidDimIndex
            varargout = {this.get_dims(S(1).subs{:})};
        else
            % use builtin indexing and hope for the best
            varargout = {builtin('subsref',this,S)};
        end
    otherwise
        % use builtin indexing for Y{i,j}
        varargout = {builtin('subsref',this,S)};
end

author: translationalneuromodeling

repo: tapas


text:
function [o,count,SSQ,S4M] = sumskipnan(i,DIM)
% SUMSKIPNAN adds all non-NaN values.
%
% All NaN's are skipped; NaN's are considered as missing values.
% SUMSKIPNAN of NaN's only  gives O; and the number of valid elements is return.
% SUMSKIPNAN is also the elementary function for calculating
% various statistics (e.g. MEAN, STD, VAR, RMS, MEANSQ, SKEWNESS,
% KURTOSIS, MOMENT, STATISTIC etc.) from data with missing values.
% SUMSKIPNAN implements the DIMENSION-argument for data with missing values.
% Also the second output argument return the number of valid elements (not NaNs)
%
% Y = sumskipnan(x [,DIM])
% [Y,N,SSQ] = sumskipnan(x [,DIM])
%
% DIM	dimension
%	1 sum of columns
%	2 sum of rows
%	default or []: first DIMENSION with more than 1 element
%
% Y	resulting sum
% N	number of valid (not missing) elements
% SSQ	sum of squares
%
% The mean & standard error of the mean and
%	Y./N & sqrt((SSQ-Y.*Y./N)./(N.*max(N-1,0)));
% the mean square & the standard error of the mean square and
% 	SSQ./N & sqrt((S4M-SSQ.^2./N)./(N.*max(N-1,0)))
%
% features:
% - can deal with NaN's (missing values)
% - implements dimension argument.
% - compatible with Matlab and Octave
%
% see also: SUM, NANSUM, MEAN, STD, VAR, RMS, MEANSQ,
%      SSQ, MOMENT, SKEWNESS, KURTOSIS, SEM

%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

%	$Revision: 1.23 $
%	$Id: sumskipnan.m,v 1.23 2003/10/31 18:15:38 schloegl Exp $
%    Copyright (C) 2000-2003 by Alois Schloegl <a.schloegl@ieee.org>

if nargin<2,
        DIM = [];
end;

% an efficient implementation in C of the following lines
% could significantly increase performance
% only one loop and only one check for isnan is needed.
% A MEX-Implementation is available in sumskipnan.cpp.
%
% Outline of the algorithm:
% for { k=1,o=0,count=0; k++; k<N}
% 	if ~isnan(i(k))
% 	{ 	o     += i(k);
%               count += 1;
%		tmp    = i(k)*i(k)
%		o2    += tmp;
%		o3    += tmp.*tmp;
%       };

if isempty(DIM),
        DIM=min(find(size(i)>1));
        if isempty(DIM), DIM = 1; end;
end;
if nargout>1,
        count = sum(~isnan(i),DIM);
end;

%if flag_implicit_skip_nan, %%% skip always NaN's
i(isnan(i)) = 0;
%end;
o = sum(i,DIM);
if nargout>2,
        i = real(i).^2 + imag(i).^2;
        SSQ = sum(i,DIM);
        if nargout>3,
                S4M = sum(i.^2,DIM);
        end;
end;