text:
classdef GemCinePanel < GemVirtualPanel
    % GemCinePanel  Image axes with scrollbar for displaing a slice of a 3D volume
    %
    %     GemCinePanel contains a 2D image and a slider, which allows the user to
    %     cine through slices of a 3D images
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2014.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties (Access = protected)
        Axes
        ImageSource
        ImageParameters
        Slider
    end

    events
        MousePositionChanged
    end

    methods
        function obj = GemCinePanel(parent, image_source, image_parameters, image_overlay_axes)
            obj = obj@GemVirtualPanel(parent);

            obj.ImageSource = image_source;
            obj.ImageParameters = image_parameters;

            obj.Slider = GemSlider(obj);
            obj.AddChild(obj.Slider);

            obj.Axes = image_overlay_axes;
            obj.AddChild(obj.Axes);

            obj.AddEventListener(obj.Slider, 'SliderValueChanged', @obj.SliderValueChanged);
        end

        function Resize(obj, position)
            Resize@GemUserInterfaceObject(obj, position);

            % Resize axes
            axis_width = max(1, position(3) - obj.Slider.SliderWidth);
            axes_position = [position(1), position(2), axis_width, position(4)];
            obj.Axes.Resize(axes_position);

            % Resize slider
            slider_x = position(1) + axis_width;
            slider_position = [slider_x, position(2) - 1, obj.Slider.SliderWidth, position(4)];
            obj.Slider.Resize(slider_position);
        end

        function SetSliceNumber(obj, slice_number)
            current_slice_value = obj.Slider.SliderValue;
            if (current_slice_value ~= slice_number)
                obj.Slider.SetSliderValue(slice_number);
            end
        end

        function SetSliderLimits(obj, min, max)
            obj.Slider.SetSliderLimits(min, max);
        end

        function SetSliderSteps(obj, steps)
            obj.Slider.SetSliderSteps(steps);
        end

        function EnableSlider(obj, enabled)
            if enabled
                obj.Slider.Enable;
            else
                obj.Slider.Disable;
            end
        end

        function frame = Capture(obj, image_size, orientation)
            drawnow;

            origin = [0.5, 0.5, 0.5];
            image_limit = origin + image_size;
            switch orientation
                case GemImageOrientation.XZ
                    xlim_image = [origin(2), image_limit(2)];
                    ylim_image = [origin(3), image_limit(3)];
                case GemImageOrientation.YZ
                    xlim_image = [origin(1), image_limit(1)];
                    ylim_image = [origin(3), image_limit(3)];
                case GemImageOrientation.XY
                    xlim_image = [origin(2), image_limit(2)];
                    ylim_image = [origin(1), image_limit(1)];
            end

            frame = obj.Axes.Capture(xlim_image, ylim_image);
        end

        function ZoomTo(obj, i_limits, j_limits, k_limits)
            % Changes the current axis limits to the specified global coordinates
            % i_limits = [minimum_i, maximum_i] and the same for j, k.

            obj.Axes.ZoomTo(obj, i_limits, j_limits, k_limits);

        end

        function ClearAxesCache(obj)
            obj.Axes.ClearAxesCache;
        end

        function UpdateAxes(obj)
            obj.Axes.UpdateAxes;
        end

        function screen_coords = GetScreenCoordinates(obj)
            coords = obj.Axes.GetCurrentPoint;
            if (~isempty(coords))
                i_screen = coords(2,1);
                j_screen = coords(2,2);
                screen_coords = [i_screen, j_screen];
            else
                screen_coords = [0, 0];
            end
        end

        function current_point = GetCurrentPoint(obj)
            current_point = obj.Axes.GetCurrentPoint;
        end

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

        function global_coords = GetImageCoordinates(obj)
            coords = round(obj.GetCurrentPoint);
            if (~isempty(coords))
                orientation = obj.ImageParameters.Orientation;
                i_screen = coords(2,1);
                j_screen = coords(2,2);
                k_screen = obj.ImageParameters.SliceNumber(orientation);

                switch orientation
                    case GemImageOrientation.XZ
                        i = k_screen;
                        j = i_screen;
                        k = j_screen;
                    case GemImageOrientation.YZ
                        i = i_screen;
                        j = k_screen;
                        k = j_screen;
                    case GemImageOrientation.XY
                        i = j_screen;
                        j = i_screen;
                        k = k_screen;
                end
            else
                i = 1;
                j = 1;
                k = 1;
            end
            global_coords = obj.ImageSource.Image.LocalToGlobalCoordinates([i, j, k]);
        end
    end

    methods (Access = protected)
        function input_has_been_processed = MouseDown(obj, click_point, selection_type, src, eventdata)
            % This method is called when the mouse is clicked inside the control

            obj.UpdateStatus(true);
            input_has_been_processed = true;
        end

        function input_has_been_processed = MouseUp(obj, click_point, selection_type, src, eventdata)
            % This method is called when the mouse is released inside the control

            obj.UpdateStatus(true);
            input_has_been_processed = true;
        end

        function input_has_been_processed = MouseHasMoved(obj, click_point, selection_type, src, eventdata)
            % Mouse has moved over the figure

            obj.UpdateStatus(true);
            input_has_been_processed = true;
        end

        function input_has_been_processed = MouseDragged(obj, click_point, selection_type, src, eventdata)
            % Mouse dragged over the figure

            obj.UpdateStatus(true);
            input_has_been_processed = true;
        end

        function input_has_been_processed = MouseExit(obj, click_point, selection_type, src, eventdata)
            % This method is called when the mouse exits a control which previously
            % processed a MouseHasMoved event

            obj.UpdateStatus(false);
            input_has_been_processed = true;
        end

        function UpdateStatus(obj, in_image)
            image_coordinates = obj.GetImageCoordinates;
            notify(obj, 'MousePositionChanged', CoreEventData(GemCoordsInImage(image_coordinates, in_image)));
        end

        function SliderValueChanged(obj, ~, ~)
            obj.ImageParameters.SliceNumber(obj.ImageParameters.Orientation) = round(obj.Slider.SliderValue);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function test08 ( dim_num, n, z, ns, sample_routine, seed_init )

%*****************************************************************************80
%
%% TEST08 tests NU_MEASURE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST08\n' );
  fprintf ( 1, '  NU_MEASURE computes the NU measure of quality.\n' );
  fprintf ( 1, '  The cell volume deviation       Nu = %14f\n', ...
    nu_measure ( dim_num, n, z, ns, sample_routine, seed_init ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%kkgentbindex 'generate index for a toolbox manual'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kgentbindex.pane file
%
% Parameters:
% Toggle: alltb 'Generate Index for All toolboxes', default: 0: 'Generates Index for all toolbox manuals'
% String: tb 'Generate Single Toolbox Index', default: ' ': 'Generate a manual index for a single toolbox'
%
% Example: kkgentbindex( {'alltb',0;'tb',' '})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kgentbindex - generate index for a toolbox manual
%
%  DESCRIPTION
% .I kgentbindex
% is used to generate the index chapter of a manual, by combining the
% index information from each chapter within the manual.  Run the
% command from the directory <toolbox_dir>/manual/<Toolbox_dir>.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkgentbindex(varargin)
Inputs={};
if nargin ==0
  arglist={'',''};
elseif nargin ==1
  arglist=varargin{1};
else error('Usage: [out1,..] = kkgentbindex(arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'alltb', 0;'tb', ' '};
maxval={0,0};
minval={0,0};
istoggle=[1,1];
was_set=istoggle * 0;
paramtype={'Toggle','String'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'kgentbindex"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [mResult, fMls, fMc, fMu, fSigma, mDatPredBest, vPredBest] = plot_McEMR(mCatalog, fBinning)
% Same as calc_McCdfnormal with plotting the fitting steps and the final result
%
% [mResult, fMls, fMc, fMu, fSigma, mDatPredBest, vPredBest] = plot_McEMR(mCatalog, fBinning);
% -----------------------------------------------------------------------------------------------------
%
% Incoming variables:
% mCatalog   : EQ catalog
% fBinning   : Binning interval, usually 0.1
%
% Outgoing variables:
% mResult     : Solution matrix including
%               vProbability: maximum likelihood score
%               vMc         : Mc values
%               vX_res      : mu (of normal CDF), sigma (of normal CDF), residuum, exitflag
%               vNmaxBest   : Number of events in lowest magnitude bin considered complete
%               vABValue    : a and b-value
% fMls       : likelihood score --> best Mc
% fMc        : Best estimated magnitude of completeness
% mDatPredBest   : Matrix of non-cumulative FMD [Prediction, magnitudes, original distribution]
% vPredBest      : Matrix of non-cumulative FMD below Mc [magnitude, prediction, uncertainty of prediction]
%
% J. Woessner: jochen.woessner@sed.ethz.ch
% updated: 29.09.04

% Initialize
vProbability = [];
vMc = [];
vABValue =[];
mFitRes = [];
vX_res = [];
vNCumTmp = [];
mDataPred = [];
vPredBest = [];
vDeltaBest = [];
vX_res = [];
vNmaxBest = [];
mResult=[];
mDatPredBest = [];

% Determine exact time period
fPeriod1 = max(mCatalog.Date) - min(mCatalog.Date);

% Determine max. and min. magnitude
fMaxMag = ceil(10 * max(mCatalog.Magnitude)) / 10;

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

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

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

%% Shift to positive values
% if fMinMag ~= 0
%     fMcBound = fMcTry-fMinMag;
% end
% Loop over Mc-values
for fMc = fMcBound-0.4:0.1:fMcBound+0.8
    fMc = round(fMc, -1);
    vFMD = vFMDorg;
    vNonCFMD = vNonCFMDorg;
    vNonCFMD = fliplr(vNonCFMD);
    % Select magnitudes to calculate b- anda-value
    vSel = mCatalog.Magnitude > fMc-fBinning/2;
    if sum(vSel) >= 20
        [ fBValue, fStdDev, fAValue] =  calc_bmemag(mCatalog.Magnitude(vSel), fBinning);
        % Normalize to time period
        vFMD(2,:) = vFMD(2,:)./fPeriod1; % ceil taken out
        vNonCFMD(2,:) = vNonCFMD(2,:)./fPeriod1; % ceil removed
        % Compute quantity of earthquakes by power law
        fMaxMagFMD = max(vNonCFMD(1,:));
        fMinMagFMD = min(vNonCFMD(1,:));
        vMstep = [fMinMagFMD:0.1:fMaxMagFMD];
        vNCum = 10.^(fAValue-fBValue.*vMstep); % Cumulative number

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

        % Normalize vN
        vN = vN./fPeriod1;
        % Data selection
        % mData = Non-cumulative FMD values from GR-law and original data
        mData = [vN' vNonCFMD'];
        vSel = (mData(:,2) >= fMc);
        mDataTest = mData(~vSel,:);
        mDataTmp = mData.subset(vSel);
%         % Check for zeros in observed data
        vSelCheck = (mDataTest(:,3) == 0);
        mDataTest = mDataTest(~vSelCheck,:);
        % Choices of normalization
        fNmax = mDataTmp(1,3); % Frequency of events in Mc bin
        %fNmax = max(mDataTest(:,3));  % Use maximum frequency of events in bins below Mc
        %fNmax = mDataTest(length(mDataTest(:,1)),3); % Use frequency of events at bin Mc-0.1 -> best fit
        if (~isempty(isempty(fNmax)) &&  ~isnan(fNmax) & fNmax ~= 0 & length(mDataTest(:,1)) > 4)
            mDataTest(:,3) = mDataTest(:,3)/fNmax; % Normalize datavalues for fitting with CDF
            % Move to M=0 to fit with lsq-algorithm
            fMinMagTmp = min(mDataTest(:,2));
            mDataTest(:,2) = mDataTest(:,2)-fMinMagTmp;
            % Curve fitting: Non cumulative part below Mc
            options = optimset;
            %options = optimset('Display','off','Tolfun',1e-7,'TolX',0.0001,'MaxFunEvals', 100000,'MaxIter',10000);
            options = optimset('Display','off','Tolfun',1e-5,'TolX',0.001,'MaxFunEvals', 1000,'MaxIter',1000);
            [vX, resnorm, resid, exitflag, output, lambda, jacobian]=lsqcurvefit(@calc_normalCDF,[fSmu  fSSigma], mDataTest(:,2), mDataTest(:,3),[],[],options);
            mDataTest(:,1) = normcdf(mDataTest(:,2), vX(1), vX(2))*fNmax;
            if (length(mDataTest(:,2)) > length(vX(1,:)))
                %% Confidence interval determination
                % vPred : Predicted values of lognormal function
                % vPred+-delta : 95% confidence level of true values
                [vPred,delta] = nlpredci(@calc_normalCDF,mDataTest(:,2),vX, resid, jacobian);
            else
                vPred = NaN;
                delta = NaN;
            end % END: This section is due for errors produced with datasets less long than amount of parameters in vX
            % Results of fitting procedure
            mFitRes = [mFitRes; vX resnorm exitflag];
            % Move back to original magnitudes
            mDataTest(:,2) = mDataTest(:,2)+fMinMagTmp;
            % Set data together
            mDataTest(:,3) = mDataTest(:,3)*fNmax;
            mDataPred = [mDataTest; mDataTmp];
            % Denormalize to calculate probabilities
            mDataPred(:,1) = round(mDataPred(:,1).*fPeriod1);
            mDataPred(:,3) = mDataPred(:,3).*fPeriod1;
            vProb_ = calc_log10poisspdf2(mDataPred(:,3), mDataPred(:,1)); % Non-cumulative

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

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

            % Keep best fitting model
            if (fProbability == min(vProbability))
                vDeltaBest = delta;
                vPredBest = [mDataTest(:,2) vPred*fNmax*fPeriod1 delta*fNmax*fPeriod1]; % Gives back uncertainty
                %fMc+fMinMag : Test procedure
                mDatPredBest = [mDataPred];
           end
        else
            %disp('Not enough data');
            % Setting values
            fProbability = NaN;
            fMc = NaN;
            vX(1) = NaN;
            vX(2) = NaN;
            resnorm = NaN;
            exitflag = NaN;
            delta = NaN;
            vPred = [NaN NaN NaN];
            fNmax = NaN;
            fAValue = NaN;
            fBValue = NaN;
            vProbability = [vProbability; fProbability];
            vMc = [vMc; fMc];
            vX_res = [vX_res; vX resnorm exitflag];
%             vDeltaBest = [vDeltaBest; NaN];
%             vPredBest = [vPredBest; NaN NaN NaN];
            vNmaxBest = [vNmaxBest; fNmax];
            vABValue = [vABValue; fAValue fBValue];
        end
    end

    % Clear variables
    vNCumTmp = [];
    mModelDat = [];
    vNCum = [];
    vSel = [];
    mDataTest = [];
    mDataPred = [];
end % END of FOR fMc
% Result matrix
mResult = [mResult; vProbability vMc vX_res vNmaxBest vABValue];

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

try % Try catch block used as mDatPreBest not always calculated
    % Reconstruct vector of magnitudes from model for Period 1
    vMag = [];
    mModelFMD = [round(mDatPredBest(:,1)) mDatPredBest(:,2)];
    vSel = (mModelFMD(:,1) ~= 0); % Remove bins with zero frequency of zero events
    mData = mModelFMD(vSel,:);
    for nCnt=1:length(mData(:,1))
        fM = repmat(mData(nCnt,2),mData(nCnt,1),1);
        vMag = [vMag; fM];
    end
    % Calculate KS-Test
    [bH,fPval,fKsstat] = kstest2(roundn(mCatalog.Magnitude,-1),roundn(vMag,-1),0.05,0)
catch
    bH = NaN;
    fPval = NaN;
    fKsstat = NaN;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Plot result

% Plot best fitting CDF vs. magnitude
figure_w_normalized_uicontrolunits('tag','best_cdfit','Name','CDF fit','Units','normalized','Nextplot','add',...
                    'Numbertitle','off','visible','on');
p1=plot(mDatPredBest(:,2), mDatPredBest(:,3),'Marker','d','MarkerEdgeColor',[0 0 0],'MarkerFaceColor',[0.4 0.4 0.4],'Markersize',8,'Linewidth',2,'Linestyle','none');
set(gca,'NextPlot','add');
p2=plot(mDatPredBest(:,2),mDatPredBest(:,1),'Marker','o','MarkerEdgeColor',[0.25 0.25 0.25],'MarkerFaceColor',[0.75 0.75 0.75],'Markersize',8,'Linewidth',2,'Linestyle','none');
%sTitlestr = ['mu = ' num2str(fMu) ', sigma = ' num2str(fSigma)];
%title(sTitlestr)
p3=plot(vPredBest(:,1),vPredBest(:,2),'-','Color',[0.2 0.2 0.2],'Linewidth',2);
p4=plot(vPredBest(:,1),vPredBest(:,2)+vPredBest(:,3),'--','Color',[0.2 0.2 0.2],'Linewidth',2);
plot(vPredBest(:,1),vPredBest(:,2)-vPredBest(:,3),'--','Color',[0.2 0.2 0.2],'Linewidth',2);
xlabel('Magnitude','FontSize',14,'FontWeight','bold');
ylabel('Number of events','FontSize',14,'FontWeight','bold');
xlim = ([min(mDatPredBest(:,2)) max(mDatPredBest(:,2))]);
ylim = ([0 max(mDatPredBest(:))]);
hl1=legend([p1,p2,p3,p4],'Noncum. FMD','Noncum. model FMD ','Model','Model uncertainty');
set(hl1,'FontWeight','bold','FontSize',12')
set(gca,'visible','on','FontSize',12,'FontWeight','bold','LineWidth',2,'Box','on');
set(gca,'NextPlot','replace');

% Plot Non-/Cumulative distribution, original and predicted
figure_w_normalized_uicontrolunits('tag','ncumdist','Name','Best model','Units','normalized','Nextplot','add',...
                    'Numbertitle','off','visible','on');
p1 = semilogy(vNonCFMDorg(1,:)', vNonCFMDorg(2,:)', 'Marker','d','MarkerEdgeColor',[0 0 0],'MarkerFaceColor',[0.4 0.4 0.4],'Markersize',8,'Linewidth',2,'Linestyle','none');
set(gca,'NextPlot','add');
p2 = semilogy(vFMDorg(1,:)', vFMDorg(2,:)', 'Marker','^','MarkerEdgeColor',[0 0 0],'MarkerFaceColor',[0.75 0.75 0.75],'Markersize',8,'Linewidth',2,'Linestyle','none');
%p2=semilogy(vCFMD(1,:)', vNBest.*fPeriod1, '*',
p3=semilogy(mDatPredBest(:,2),mDatPredBest(:,1),'Marker','o','MarkerEdgeColor',[0.25 0.25 0.25],'MarkerFaceColor',[0.75 0.75 0.75],'Markersize',8,'Linewidth',2,'Linestyle','none');
%fBvalue = mResult(vSel,9)
%fAValue = mResult(vSel,8)
vPoly = [-1*fBvalue fAvalue]
vMagnitudes = [0:0.01:max(mDatPredBest(:,2))];
fBFunc = 10.^(polyval(vPoly, vMagnitudes));
p4 = semilogy(vMagnitudes, fBFunc,'Linewidth',3,'Linestyle','--','Color',[0.4 0.4 0.4])
xlim = ([min(mDatPredBest(:,2)) max(mDatPredBest(:,2))]);
ylim = ([0 100000]);
%sTitlestr = ['Mc = ' num2str(fMc) ' using Normal CDF fitting'];
%title(sTitlestr)
sText1 = ['KS-Test: H = ' num2str(bH)];
sText2 = ['Mc(EMR) = ' num2str(fMc)];
hT1 = text(max(vFMDorg(1,:))-1, max(vFMDorg(2,:))*0.6,sText1);
hT2 = text(max(vFMDorg(1,:))-1, max(vFMDorg(2,:))*0.4,sText2);
set(hT1,'FontSize',12,'FontWeight','bold');
set(hT2,'FontSize',12,'FontWeight','bold');
xlabel('Magnitude','FontSize',14,'FontWeight','bold');
ylabel('Number of events','FontSize',14,'FontWeight','bold');
hl1=legend([p1 p2 p3 p4],'Noncum. FMD','Cum. FMD','Noncum. model FMD','Cum. model FMD');
set(hl1,'FontWeight','bold','FontSize',12')
set(gca,'visible','on','FontSize',12,'FontWeight','bold','LineWidth',2,'Box','on');

%%%% Plotting the likelihood scores
figure
%vSel2 = (vMc >= 0.6 & vMc<= 2);
%vMc = vMc(vSel2);
%vProbability = vProbability(vSel2);
plot(vMc, vProbability,'Marker','s','MarkerFaceColor',[0.5 0.5 0.5],'MarkerEdgeColor',[0 0 0],'Markersize',10','Linewidth',2,'Color',[0 0 0],'visible','on');%sTitlestr = ['mu = ' num2str(fMu) ', sigma = ' num2str(fSigma)];
hl2=legend('MLE');
set(hl2,'FontWeight','bold','FontSize',12');
xlim = ([min(mDatPredBest(:,2)) max(mDatPredBest(:,2))]);
set(gca,'visible','on','FontSize',12,'FontWeight','bold','LineWidth',2,'Box','on');
xlabel('Magnitude','FontSize',14,'FontWeight','bold');
ylabel('Likelihood','FontSize',14,'FontWeight','bold');

author: CelsoReyes


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ============================================
% patch_light.m
%
% Goal: surface rendering with lighting
%       need to open figure or subplot before use
%
% Li Shen
% 02/21/2006 - create

function h = patch_lightmesh(vertices,faces)

% different color options
fcolor = [0.55 0.5 0.55];
fcolor = [0.3 0.2 0.2]*2;
fcolor = [0.6 0.4 0.4];
fcolor = [0.5 0.2 0.2];
fcolor = [0.5 0.5 0.5];
fcolor = [0.6 0.3 0.3];
h = patch('faces', faces, 'vertices', vertices,'FaceColor', fcolor, 'EdgeAlpha', 1, 'EdgeColor','k', 'FaceLighting', 'phong', ...
      'AmbientStrength', 0.9, 'SpecularColorReflectance', 0.5, 'DiffuseStrength',0.5,'SpecularExponent',5);

c = [1 1 1]*0.5;
light('Position',[ 1  0 0],'Style','infinite', 'Color', c);
light('Position',[-1  0 0],'Style','infinite', 'Color', c);
light('Position',[ 0  1 0],'Style','infinite', 'Color', c);
light('Position',[ 0 -1 0],'Style','infinite', 'Color', c);
lighting gouraud;

axis image; box on; view(3);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function value = year_length_bahai ( y )

%*****************************************************************************80
%
%% YEAR_LENGTH_BAHAI returns the number of days in a Bahai year.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year to be checked.
%
%    Output, integer VALUE, the number of
%    days in the year.
%
  if ( year_is_leap_bahai ( y ) )
    value = 366;
  else
    value = 365;
  end

  return
end


text:
function test_example_phase_reset

% MEM 4gb
% WALLTIME 00:10:00

%
%% Simulate an oscillatory signal with phase resetting
%
%% # Narrow-band oscillation
%
% The following code simulates a clean oscillating signal with a phase reset at time zero.
%
clear all
close all

f = 5;
t = (-1000:1:1000)/1000;

figure
for i=1:10 % repeat 10 times
s = nan(size(t));
s(t<0 ) = cos(2*pi*f*t(t<0 ) + 2*pi*rand(1)); % random phase
s(t>=0) = cos(2*pi*f*t(t>=0)               ); % constant phase
s = 0.2*s+i;
d(i,:) = s; % remember the signal on each repetition
plot(t, s);
hold on;
end
axis auto

figure;
plot(t, mean(d,1))

%
%% # Broad-band oscillation
%
% The following code also includes a small "random walk" in the phase, i.e. the signal is a little bit broad-band and over time there is some phase dispersion. At t=0 there is still a phase reset. The phase dispersion causes the average ERF again to disappear some time following the phase reset.
%
clear all
close all
f = 5;
t = (-2000:1:5000)/1000;

figure
for i=1:10 % repeat 10 times

p = 2*pi*f*t; % linear phase increase
p(t<0) = p(t<0) + 2*pi*rand(1); % random phase prior to TMS pulse
p = p + cumsum(2*pi*(rand(size(p))-0.5)*0.01); % random walk
p = p - p(t==0);

s = cos(p);
s = 0.2*s+i;
d(i,:) = s; % remember the signal on each repetition
plot(t, s);
hold on;
end
axis auto

figure;
plot(t, mean(d,1))

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function b = isParallel(v1, v2, varargin)
%ISPARALLEL Check parallelism of two vectors.
%
%   B = isParallel(V1, V2)
%   where V1 and V2 are two row vectors of length ND, ND being the
%   dimension, returns 1 if the vectors are parallel, and 0 otherwise.
%
%   Also works when V1 and V2 are two N-by-ND arrays with same number of
%   rows. In this case, return a N-by-1 array containing 1 at the positions
%   of parallel vectors.
%
%   Also works when one of V1 or V2 is N-by-1 and the other one is N-by-ND
%   array, in this case return N-by-1 results.
%
%   B = isParallel(V1, V2, ACCURACY)
%   specifies the accuracy for numerical computation. Default value is
%   1e-14.
%
%
%   Example
%   isParallel([1 2], [2 4])
%   ans =
%       1
%   isParallel([1 2], [1 3])
%   ans =
%       0
%
%   See also
%   vectors2d, isPerpendicular, lines2d
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2006-04-25
% Copyright 2006-2022 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)

% default accuracy
acc = 1e-14;
if ~isempty(varargin)
    acc = abs(varargin{1});
end

% normalize vectors
v1 = normalizeVector(v1);
v2 = normalizeVector(v2);

% adapt size of inputs if needed
n1 = size(v1, 1);
n2 = size(v2, 1);
if n1 ~= n2
    if n1 == 1
        v1 = v1(ones(n2,1), :);
    elseif n2 == 1
        v2 = v2(ones(n1,1), :);
    end
end

% performs computation
if size(v1, 2) == 2
    % computation for plane vectors
    b = abs(v1(:, 1) .* v2(:, 2) - v1(:, 2) .* v2(:, 1)) < acc;
else
    % computation in greater dimensions
    b = vectorNorm(cross(v1, v2, 2)) < acc;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function h=editpm1(fig,x,y,dx,dy,str,clb)

clr=[0.8313725490196079   0.8156862745098039   0.7843137254901961]; % background color

h=uicontrol('parent',fig,...
    'Style','edit','units','normalized','fontunits','normalized',...
    'BackgroundColor',[1 1 1],...
    'HorizontalAlignment','center',...
    'fontsize',0.8,...
    'position',[x y-dy/2 dx dy],'string',str,...
    'callback',clb);

st=0.1;

p=uicontrol('parent',fig,...
    'Style','pushbutton','units','normalized','fontunits','normalized',...
    'BackgroundColor',clr,...
    'HorizontalAlignment','center',...
    'fontsize',1,...
    'position',[x+dx y-dy/2+dy/2 dy/2 dy/2],'string','+',...
    'callback',['stp1(' num2str(h,'%20.20e') ',' num2str(st,'%20.20e') ',''' clb ''',' num2str(fig,'%20.20e') ')' ]);

m=uicontrol('parent',fig,...
    'Style','pushbutton','units','normalized','fontunits','normalized',...
    'BackgroundColor',clr,...
    'HorizontalAlignment','center',...
    'fontsize',1,...
    'position',[x+dx y-dy/2 dy/2 dy/2],'string','-',...
    'callback',['stp1(' num2str(h,'%20.20e') ',' num2str(-st,'%20.20e') ',''' clb ''',' num2str(fig,'%20.20e') ')' ]);

author: Sable

repo: mcbench-benchmarks


text:
function myribbon = riverplot_draw_ribbon(topleft, topright, bottomleft, bottomright, varargin)

%
% This is a subfunction of the riverplot toolbox.
%  - Draws a ribbon given four [x y] coordinate pairs specifying the
%  corners of the ribbon in clockwise order.
%  - Used in riverplot_ribbon.m
%
% :Usage:
% ::
%
%     myribbon = riverplot_draw_ribbon((topleft, topright, bottomleft, bottomright, [optional inputs])
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2016  Tor Wager
%
%     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/>.
% ..
%
% :Inputs:
%
%   **rect1:**
%        Structure with handles and decription of imrect class object,
%        created with riverplot_rect.m.  This is always the LEFTMOST
%        rect of the pair.
%
%   **rect2:**
%        Same as rect1, but this is the RIGHTMOST rect of the pair
%
% :Optional Inputs:
%   **coverage:**
%        'coverage' followed by 2-element vector for [rect1 rect2], where
%        each value is between 0 and 1 and specifies the percentage of the
%        plot the ribbon will cover.
%
%   **color:**
%        followed by [r g b] color spec for ribbon
%
%   **ax:**
%        followed by axis handle to plot
%
%   **from_bottom:**
%        align ribbon with bottom edges of rects
%
%   **steepness**
%        followed by value for sigmoid steepness, on [0 1] inteval
%
% :Outputs:
%
%   **myribbon:**
%        structure with handles and info for ribbon
%
%
% :Examples:
% ::
%
% rect1 = riverplot_rect('color', [.7 .3 .3], 'position', [1 2 1 3]);
% rect2 = riverplot_rect('color', [.3 .3 .6], 'position', [5 4 1 3]);
% newbr = rect2.topleft;
% newbr(2) = newbr(2) - 1.5;
% myribbon = riverplot_draw_ribbon(rect1.topright, rect2.topleft, rect1.bottomright, newbr, 'color', [.7 .3 .3]);
%
% :References:
%   CITATION(s) HERE
%
% :See also:
%   - list other functions related to this one, and alternatives*
%

% ..
%    Programmers' notes:
%    Created by Tor Wager, July 2, 2016
% ..

% ..
%    DEFAULTS AND INPUTS
% ..

ax = gca;
%myposition = [1 3 2 5];
mycolor = [.3 .3 .8];
steepness = .05;

% optional inputs with default values
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case 'ax', ax = varargin{i+1}; varargin{i+1} = [];
            %case 'position', myposition = varargin{i+1}; varargin{i+1} = [];
            case 'color', mycolor = varargin{i+1}; varargin{i+1} = [];

            case 'steepness', steepness = varargin{i+1}; varargin{i+1} = [];

            case {'coverage', 'from_bottom', 'y_offset'} % pass it on, irrelevant but passed from higher level calling fcn

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

axes(ax)

% each line coordinate is an [x,y] pair, e.g., topleft = [x,y] for top left
% point in absolute plot coordinates

line1 = riverplot_line(topleft, topright, 'w', 1, steepness);

line2 = riverplot_line(bottomleft, bottomright, 'w', 1, steepness);

ph = patch([line1.xcoords; line2.xcoords(end:-1:1)], [line1.ycoords; line2.ycoords(end:-1:1)], mycolor);

myribbon = struct('line1', line1, 'line2', line2, 'patchh', ph);

set(line1.h, 'Color', mycolor);
set(line2.h, 'Color', mycolor);
set(ph, 'FaceAlpha', .4);

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef ZXH_CF10 < PROBLEM
% <multi/many> <real> <large/none> <constrained>
% Constrained benchmark MOP proposed by Zhou, Xiang, and He

%------------------------------- Reference --------------------------------
% Y. Zhou, Y. Xiang, and X. He, Constrained multiobjective optimization:
% Test problem construction and performance evaluations, IEEE Transactions
% on Evolutionary Computation, 2021, 25(1): 172-186.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        k;  % Number of constrained variables
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = obj.M+10;  end
            obj.lower    = zeros(1,obj.D) + 1e-10;
            obj.upper    = ones(1,obj.D)  - 1e-10;
            obj.encoding = ones(1,obj.D);
            if obj.M <= 3
                obj.k = obj.M - 1;
            elseif obj.M > 3 && obj.M <= 8
                obj.k = floor(obj.M/2);
            else
                obj.k = 3;
            end
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            PopDec = varargin{1};
            OptX   = 0.2;
            [N,D]  = size(PopDec);
            M      = obj.M;
            % Step 1: Compute cumsum
            Sx = cumsum(PopDec(:,1:M).^2,2,'reverse');
            % Step 2: Compute theta
            THETA = 2/pi*atan(sqrt(Sx(:,2:end))./PopDec(:,1:M-1));
            % Step 3: Calculate Griewank function
            h = 5*(sum((PopDec(:,M+1:end)-OptX).^2,2)-prod(cos(10*pi*(PopDec(:,M+1:end)-OptX)./repmat(sqrt(1:(D-M)),N,1)),2)+1);
            % Step 4: Compute T_
            T = (1 - Sx(:,1)).^2 + h;
            % Step 5: Objectives (convex)
            G      = 1-[ones(N,1) cumprod(sin(pi/2*THETA),2)] .* [cos(pi/2*THETA) ones(N,1)];
            PopObj = G .* repmat((1+T),1,M);
            % Step 6: Constraints
            PopCon(:,1) = Sx(:,1) + h - 1;
            PopCon(:,2) = -(Sx(:,1) + h - 1/4);
            for i = 1 : obj.k
                PopCon(:,i+2) = min(THETA(:,i)-1/4,3/4-THETA(:,i));
            end
            Population = SOLUTION(varargin{1},PopObj,PopCon,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,obj.M);
            R = 1 - R./repmat(sqrt(sum(R.^2,2)),1,obj.M);
            T = zeros(size(R));
            for i = obj.M-1 : -1 : 1
                T(:,i) = atan((1-R(:,i+1))./(1-R(:,i))./cos(T(:,i+1)));
            end
            THETA = T(:,1:obj.k)*2/pi;
            Valid = all(THETA<=1/4|THETA>=3/4,2);
            R     = R(Valid,:);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a  = linspace(0,pi/2,30)';
                R  = {1-sin(a)*cos(a'),1-sin(a)*sin(a'),1-cos(a)*ones(size(a'))};
                T2 = atan((1-R{3})./(1-R{2}));
                T1 = asin((1-R{2})./cos(T2));
                THETA = cat(3,T1,T2)*2/pi;
                Valid = all(THETA<=1/4|THETA>=3/4,3);
                R{1}(~Valid) = nan;
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% Undistort points using camera parameters
%
% This function is much faster than Matlab's undistortPoints and works on large
% point sets. The code is based on cvUndistortPoints function from OpenCV library.
% Also this function works with NaN values.
%
%  USAGE
%   undistortedPoints = undistortPoints(points, cameraParams)
%   points              Nx2 ([x y]) matrix of points.
%   cameraParams        Camera calibration information.
%   undistortedPoints   Nx2 ([x y]) matrix with undistorted points.
%
function undistortedPoints = undistortPoints(points, cameraParams)
    origin = [cameraParams.IntrinsicMatrix(3, 1:2)];
    cx = origin(1);
    cy = origin(2);
    fx = cameraParams.IntrinsicMatrix(1, 1);
    fy = cameraParams.IntrinsicMatrix(2, 2);
    ifx = 1 / fx;
    ify = 1 / fy;
    iters = 50;
    k = cameraParams.RadialDistortion;
    p_coef = cameraParams.TangentialDistortion;

    if length(k) == 2
        k(3) = 0;
    end

    x = points(:, 1);
    y = points(:, 2);
    x = (x - cx) .* ifx;
    y = (y - cy) .* ify;
    x0 = x;
    y0 = y;

    for j = 1:iters
        r2 = x.*x + y.*y;
        icdist = 1 ./ (1 + ((k(3).*r2 + k(2)).*r2 + k(1)) .* r2);
        deltaX = 2.*p_coef(1).*x.*y + p_coef(2).*(r2 + 2.*x.*x);
        deltaY = p_coef(1) .* (r2 + 2.*y.*y) + 2.*p_coef(2).*x.*y;
        x = (x0 - deltaX) .* icdist;
        y = (y0 - deltaY) .* icdist;
    end

    undistortedPoints(:, 1) = round(x .* fx + cx);
    undistortedPoints(:, 2) = round(y .* fy + cy);
end

author: kavli-ntnu


text:
function lcvt_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests CVT, R8MAT_LATINIZE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2006
%
%  Author:
%
%    John Burkardt
%
  m = 2;
  n = 25;
  latin_steps = 3;
  sample_function_cvt = 0;
  sample_function_init = 3;
  sample_num_cvt = 100000;
  sample_num_steps = 50;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  CVT computes a Centroidal Voronoi Tessellation.\n' );
  fprintf ( 1, '  R8MAT_LATINIZE makes it a Latin Hypersquare.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this test, we initialize the generators to\n' );
  fprintf ( 1, '  grid points; this is an unstable CVT solution.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension M =        %d\n', m );
  fprintf ( 1, '  Number of generators N =     %d\n', n );
  fprintf ( 1, '  Initial random number seed = %d\n', seed );
  fprintf ( 1, '\n' );

  if ( sample_function_init == -1 )
    fprintf ( 1, '  Initialize using RAND (MATLAB intrinsic).\n' );
  elseif ( sample_function_init == 0 )
    fprintf ( 1, '  Initialize using UNIFORM.\n' );
  elseif ( sample_function_init == 1 )
    fprintf ( 1, '  Initialize using HALTON.\n' );
  elseif ( sample_function_init == 2 )
    fprintf ( 1, '  Initialize using GRID.\n' );
  elseif ( sample_function_init == 3 )
    fprintf ( 1, '  USER will initialize data.\n' );
  end

  if ( sample_function_cvt == -1 )
    fprintf ( 1, '  Sample using RAND (MATLAB intrinsic).\n' );
  elseif ( sample_function_cvt == 0 )
    fprintf ( 1, '  Sample using UNIFORM.\n' );
  elseif ( sample_function_cvt == 1 )
    fprintf ( 1, '  Sample using HALTON.\n' );
  elseif ( sample_function_cvt == 2 )
    fprintf ( 1, '  Sample using GRID.\n' );
  end

  fprintf ( 1, '  Number of sample points = %d\n', sample_num_cvt );
  fprintf ( 1, '  Number of sample steps =  %d\n', sample_num_steps );

  ngrid = 5;

  for rank = 0 : n-1
    tuple = tuple_next_fast ( ngrid, m, rank );
    generator(1:m,rank+1) = ( 2 * tuple(1:m)' - 1 ) / ( 2 * ngrid );
  end

  r8mat_transpose_print ( m, n, generator, '  Initial generators (rows):' );

  for i = 1 : latin_steps

    [ generator, seed ] = cvt ( m, n, sample_function_init, ...
      sample_function_cvt, sample_num_cvt, sample_num_steps, seed, generator );

    r8mat_transpose_print ( m, n, generator, '  After CVT steps:' );

    generator = r8mat_latinize ( m, n, generator );

    r8mat_transpose_print ( m, n, generator, '  After Latin step:' );

  end

  return
end


text:
classdef Event < Data
    properties (Access = public)
        sample;
        label;
    end

    methods (Access = public)
        function obj = Event(sample,label)
            if nargin == 2
                obj.sample = sample;
                obj.label = label;
            end
            obj.type = DataType.kEvent;
        end
    end

    methods (Access = public, Static)

        %creates an array of events with the provided locations and an
        %empty label
        function events = EventsArrayWithEventLocations(eventLocations)
            nEvents = length(eventLocations);
            events = repmat(Event,1,nEvents);

            for i = 1 : nEvents
                eventLocation = eventLocations(i);
                events(i) = Event(eventLocation,[]);
            end
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function [mXSVolcanoes] = volcano_Transform2CrossSection(mVolcanoes, fWidth, fLength, fLat1, fLon1, fLat2, fLon2)

vSel = (mVolcanoes(:,3) == 1);
mType1 = mVolcanoes(vSel,:);

vSel = (mVolcanoes(:,3) == 2);
mType2 = mVolcanoes(vSel,:);

[vType1X, vType1Y] = lc_xsec2(mType1(:,2)', mType1(:,1)', mType1(:,1)*0, fWidth, fLength, fLat1, fLon1, fLat2, fLon2);

mType1 = [vType1X' vType1Y];
mType1(:,3) = 1;

[vType2X, vType2Y] = lc_xsec2(mType2(:,2)', mType2(:,1)', mType2(:,1)*0, fWidth, fLength, fLat1, fLon1, fLat2, fLon2);

mType2 = [vType2X' vType2Y];
mType2(:,3) = 2;

mXSVolcanoes = [mType1; mType2];


text:
function imgFrame(fig,file,ax)

% IMGFRAME  Exports figure graphics to image file.
%   IMGFRAME(FIG,FILE) creates from the graphics contents of FIG
%   a image file FILE.
%
%   IMGFRAME(FIG,FILE,AX) uses the contents of axes AX in
%   figure FIG instead.
%
%   Examples:
%       imgFrame(1,'image.png',gca) exports the current axes in figure 1.
%       imgFrame(gcf,'image.png') exports the current figure.
%       imgFrame(1,sprintf('image-%03d.png',t)) for loops governed by index
%           t, exports the figure 1 at time t. Use it to make movies.
%
%   See also GETFRAME, FRAME2IM, IMWRITE.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

% figure(fig)
if nargin < 3
    frmh = fig;
else
    frmh = ax;
end

% drawnow
frm = getframe(frmh);
im  = frame2im(frm);
imwrite(im,file);

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
clear;
		% Setarea formatului de scriere a randurilor
format compact
        % Definirea functiei de integrat
F=inline('1./x');
		% Incarcarea vectorului x cu elementele cuprinse
		% intre limitele de integrare
x=1:1e-2:2;
		% Evaluarea functiei de integrat in aceste puncte
y=feval(F,x);
		% Incarcarea solutiei exacte
sol_exact=log(2);
        % Calcularea aproximativa a integralei cu
		% ajutorul metodei lui Simpson
I_S=quad(F,1,2);
		% Calcularea erorii relative la aceasta metoda
eroare_S=abs(I_S-sol_exact)/sol_exact*100;
		% Afisarea rezultatelor
disp(['Solutia exacta a integralei= ',num2str(sol_exact)]);
disp(' ');
disp('Solutia aproximativa a integralei');
disp(['prin metoda lui Simpson=', num2str(I_S)]);
disp(['Eroarea relativa= ',num2str(eroare_S),' %']);
		% Calcularea aproximativa a integralei cu
		% ajutorul metodei lui Lobatto
I_L=quadl(F,1,2);
		% Calcularea erorii relative la aceasta metoda
eroare_L=abs(I_L-sol_exact)/sol_exact*100;
		% Afisarea rezultatelor
disp(' ');
disp('Solutia aproximativa a integralei ');
disp(['prin metoda lui Lobatto=', num2str(I_L)]);
disp(['Eroarea relativa= ',num2str(eroare_L),' %']);
		% Calcularea aproximativa a integralei cu
		% ajutorul metodei trapezelor
I_t=trapz(x,y);
		% Calcularea erorii relative la aceasta metoda
eroare_t=abs(I_t-sol_exact)/sol_exact*100;
		% Afisarea rezultatelor
disp(' ');
disp('Solutia aproximativa a integralei ');
disp(['prin metoda trapezelor=', num2str(I_t)]);
disp(['Eroarea relativa= ',num2str(eroare_t),' %']);
		% Generarea unui vector cu erorile obtinute
		% la diferite metode
eror=[eroare_S, eroare_L,eroare_t];
		% Determinarea erorii minime
		%  - val_eroare contine valoarea minima
		%    din vectorul eror
		%  - nr_eroare contine numarul elementului
		%    minim din vectorul eror
[val_eroare, nr_eroare]=min(eror);
		% Afisarea erorii minime si a metodei cu ajutorul
		% careia s-a obtinut
disp(' ');
disp(['Eroarea cea mai mica de ',num2str(val_eroare)]);
switch nr_eroare
    case 1
	    disp('s-a obtinut cu ajutorul metodei Simpson');
    case 2
	    disp('s-a obtinut cu ajutorul metodei Lobatto');
    case 3
	    disp('s-a obtinut cu ajutorul metodei trapezelor');
end;


text:
function a = helmert_inverse ( n )

%*****************************************************************************80
%
%% HELMERT_INVERSE returns the inverse of the HELMERT matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the inverse matrix.
%
  a = ( helmert ( n ) )';

  return
end


text:
function [G_D2,MSE]=GenerateGD2Table(numiter,N,alphaNT,G_D2init,Range_prior,Range_post)
%function [G_D2,MSE]=GenerateGD2Table(numiter,N,alphaNT,Range_prior,Range_post);
% Generates gain table for MMSE noise power estimation. The recursive
% nature of the estimation problem is taken into account by an iterative
% scheme. For more details see:
%
% J.S. Erkelens and R. Heusdens, "Tracking of nonstationary noise based on
% data-driven recursive noise power estimation", IEEE Trans. Audio,
% Speech & Lang. Proc., vol. 16, no. 6, pp. 1112-1123, 2008.
%
% J.S. Erkelens and R. Heusdens, "Fast noise tracking based on recursive
% smoothing of MMSE noise power estimates", ICASSP 2008, pp. 4873-4876.
%
% INPUT VARIABLES:
% numiter : number of optimization iterations.
% N : frame length.
% alphaNT : smoothing parameter in prior SNR parameter for noise tracking.
% Range_prior : range of the prior SNR parameter, e.g, [-19 40] (as used
% in the papers) means that the prior SNR is quantized between -19 dB and
% 40 dB in steps of 1 dB.
% Range_post: range of the posterior SNR. In the papers [-30 40] is used,
% meaning that the gain table G_D2 has dimensions 60-by-71.
% G_D2init : Initial gain table for MMSE estimation of noise power. May be
% found with the basic data-driven method. Simply initializing with
% ones(60,71) should also work (for Range_prior=[-19 40] and
% Range_post=[-30 40]).
%
% OUTPUT VARIABLES:
% G_D2 : gain table for MMSE estimation of noise power.
% MSE : normalized mean-square error per iteration (NOTE: since we are
% initializing with the noisy powers, this MSE may INCREASE with each
% iteration.
% If we apply the gain functions to the training data RECURSIVELY, however,
% the mse decreases with each iteration step. Both MSEs (on the traindata
% and when applied recursively) do converge to the same end value. See the
% 2nd paper mentioned above for a more thorough explanation and
% illustration.
%
% Copyright 2008: Delft University of Technology, Information and
% Communication Theory Group. The software is free for non-commercial use.
% This program comes WITHOUT ANY WARRANTY.
%
% Last modified: 17-3-2008.

global NoiseSpectrum
load NoiseSpectrum

Gd2new=G_D2init;
% Intialize MSE
MSE=zeros(1,numiter);
% SNRs (dB scale)
SNRs=[-12.5:5:27.5];
% Directory with training data
DIR='C:/data/Timit_TBW/Matrices/';
% List of training data files
fid=fopen('C:/data/Timit_TBW/trainlist.txt','r');
% read the first filename
s=fgetl(fid);
for k1=1:numiter
    % Initialize statistics needed to compute the optimal gain values
    sumD2R2=zeros(Range_prior(2)-Range_prior(1)+1,Range_post(2)-Range_post(1)+1);
    sumR4=sumD2R2;sumD4=sumD2R2;
    while s~=-1 % while not EOF
        % Use 2 Timit directories for training
        if strcmpi(s(16:18),'dr4') | strcmpi(s(16:18),'dr5')
            % Strip some path info and convert to lowercase
            s=lower(s(16:end-3));
            % Some difference between Windows and Linux
            I=find(s=='\');s(I)='/';
            % Collect statistics for all SNRs
            for snrnr=1:length(SNRs)
                % strip path info, retain file name
                filestr=s(11:end-1);
                % directory name for arrays of frame numbers with
                % sufficient speech energy.
                dirstr=[DIR s(1:10)];
                % load indices of frames with sufficient speech energy.
                % Only those are used for optimization of the gain.
                eval(['load ' dirstr filestr '_IsufE'])
                % directory name for the data at current SNR
                dirstr=[DIR s(1:10) 'Snr' num2str(snrnr) '/'];
                % load noise powers (needed to compute the statistics
                % for gain optimization).
                iterstr='Init';matstrD='Dmatrix';
                filestrD=[dirstr filestr '_' matstrD 'Snr' num2str(snrnr) iterstr];
                eval(['load ' filestrD ' ' matstrD])
                % "Dclean" means the "true" noise powers. The final
                % estimates will of course not be "clean" but contaminated
                % by some speech power.
                Dclean=Dmatrix;
                % Load noisy amplitudes (Rmatrix)
                filestrR=[dirstr filestr '_Rmatrix' 'Snr' num2str(snrnr) 'Init'];
                eval(['load ' filestrR ' Rmatrix'])
                % Load estimates from the previous iteration (Dmatrix) ...
                iterstr='Iter';
                filestrD=[dirstr filestr '_' matstrD 'Snr' num2str(snrnr) iterstr];
                eval(['load ' filestrD ' ' matstrD])
                % ... but initialize with 'Rmatrix', i.e., noisy amplitudes
                % in the first iteration.
                if k1==1,Dmatrix=Rmatrix;end
                % Collect statistics and calculate new data
                [sumd4,sumd2r2,sumr4,Dmatrix]=collectdata(N,alphaNT,IsufE,Dclean,Dmatrix,Rmatrix,Gd2new,Range_prior,Range_post,SNRs(snrnr));
                sumD2R2=sumD2R2+sumd2r2;sumR4=sumR4+sumr4;sumD4=sumD4+sumd4;
                % Save updated Dmatrix
                iterstr='Iter';
                filestrD=[dirstr filestr '_' matstrD 'Snr' num2str(snrnr) iterstr];
                eval(['save ' filestrD ' ' matstrD])
            end % for snrnr=1:length(SNRs)
        end % strcmpi(s(16:18),'dr4') | strcmpi(s(16:18),'dr5')
        % read the next file name
        s=fgetl(fid);
    end % while s~=-1
    % Update gain function and calculate MSE
    sumR41=sumR4;
    sumR41(sumR41==0)=1;% avoid dividing by 0
    Gd2new=sumD2R2./sumR41;
    % Calculate MSE
    MSEd=sumD4-2*Gd2new.*sumD2R2+(Gd2new.^2).*sumR4;
    % Normalize
    MSEd=sum(sum(MSEd))/sum(sum(sumD4));
    MSE(k1)=MSEd;
    % Save some stuff
    alphap=0.1;alphad=0.85;T=4;% These values are used below in the function collectdata
    eval(['save ' DIR2 'GainTable_iteration' num2str(k1) ' Gd2new MSE alphaNT alphad alphap N T']);
    % display number of iterations to go and MSE
    disp(numiter-k1),disp(MSEd)
    % Go through all the data again in the next iteration
    frewind(fid);s=fgetl(fid);
end % for k1=1:numiter
fclose(fid);
G_D2=Gd2new;

function [sumd4,sumd2r2,sumr4,Dnew]=collectdata(N,alphaNT,IsufE,Dclean,Dmatrix,Rmatrix,Gd2new,Range_prior,Range_post,dBsnr)
global NoiseSpectrum
% Initialization of statistics for computing the gain table
sumd2r2=zeros(Range_prior(2)-Range_prior(1)+1,Range_post(2)-Range_post(1)+1);
sumr4=sumd2r2;
sumd4=sumd2r2;
alpha2=0.85*ones(1,N/2+1)';
alphap=0.1;alphad=alpha2(1);
T=4;
% Frequency bins used to collect the statistics from
freqbins=10:110;
% speech presence probability
Pspi=zeros(1,N/2+1)';Pspii=Pspi;
% speech presence index
Isp=zeros(1,N/2+1)';
% minimum value prior SNR
ksi_min=10^(-1.9);% -19 dB
% Initialization of noise variance for iteration i and i+1 (denoted by ii)
labda_di=NoiseSpectrum/(10^(dBsnr/10));%True noise variance.
% Appropriate when the speech is normalized as in GenerateTraindata.m
% Note that statistics are collected for bins in the passband only.
labda_dii=labda_di;
S=size(Dmatrix);
Dnew=zeros(S);Dnew(:,1)=Dmatrix(:,1);
% frame length
N=2*S(1)-2;
S=S(2);
% length of a column of the gain matrix
Gsize=size(Gd2new);Gsize=Gsize(1);
for m=2:S
    R=Rmatrix(:,m);D=Dclean(:,m);
    % Update prior SNR (priorNT) and posterior SNR (postSNR) using data
    % from the last iteration
    postSNR=R.^2./labda_di;
    B0=Rmatrix(:,m-1).^2;
    priorNT=max(alphaNT*B0./labda_di+(1-alphaNT)*postSNR,ksi_min);
    dBprior=10*log10(priorNT);
    dBpost=10*log10(postSNR+eps);
    % Indices in gain table
    Iprior=min(max(round(dBprior),Range_prior(1)),Range_prior(2))-Range_prior(1)+1;
    Ipost=min(max(round(dBpost),Range_post(1)),Range_post(2))-Range_post(1)+1;
    % Update labda_di
    D20=Dmatrix(:,m).^2;
    Mpost=freqsmooth(postSNR,1);Isp=(Mpost>T);
    Pspi=alphap.*Pspi+(1-alphap).*Isp;
    alpha2i=alphad+(1-alphad).*Pspi;
    labda_di=alpha2i.*labda_di+(1-alpha2i).*D20;
    % Calculate new estimated noise powers
    D2=GainR2(Gd2new,Iprior,Ipost,R,N);
    Dnew(:,m)=sqrt(D2);
    % Calculate new priorNT en postSNR
    postSNR=R.^2./labda_dii;
    priorNT=max(alphaNT*B0./labda_dii+(1-alphaNT)*postSNR,ksi_min);
    dBprior=10*log10(priorNT);
    dBpost=10*log10(postSNR+eps);
    Iprior=min(max(round(dBprior),Range_prior(1)),Range_prior(2))-Range_prior(1)+1;
    Ipost=min(max(round(dBpost),Range_post(1)),Range_post(2))-Range_post(1)+1;
    % Collect data from bins 10 t/m 110 (because telephone bandwidth speech
    % is limited to about 300-3400 Hz) for frames with sufficient speech
    % energy.
    if ~isempty(intersect(m,IsufE))
        %CHECK FOR MULTIPLE ELEMENTS
        index2=Iprior(freqbins)+Gsize*(Ipost(freqbins)-1);
        if length(unique(index2))==length(index2)
            sumd4(index2)=sumd4(index2)+D(freqbins).^4;
            sumr4(index2)=sumr4(index2)+R(freqbins).^4;
            sumd2r2(index2)=sumd2r2(index2)+(R(freqbins).*D(freqbins)).^2;
        else
            for k1=1:length(freqbins)
                index2=Iprior(freqbins(k1))+Gsize*(Ipost(freqbins(k1))-1);
                sumd4(index2)=sumd4(index2)+D(freqbins(k1)).^4;
                sumr4(index2)=sumr4(index2)+R(freqbins(k1)).^4;
                sumd2r2(index2)=sumd2r2(index2)+(R(freqbins(k1)).*D(freqbins(k1))).^2;
            end
        end
    end
    % Update labda_d with new estimate of D^2
    D20=Dnew(:,m).^2;
    Mpost=freqsmooth(postSNR,1);Isp=(Mpost>T);
    Pspii=alphap.*Pspii+(1-alphap).*Isp;
    alpha2ii=alphad+(1-alphad).*Pspii;
    labda_dii=alpha2ii.*labda_dii+(1-alpha2ii).*D20;
end

function Amplitude2=GainR2(Gmatrix,Iprior,Ipost,R,N)
L=size(Gmatrix);L=L(1);
index=1:N/2+1;index2=Iprior(index)+L*(Ipost(index)-1);
Amplitude2=Gmatrix(index2).*(R.^2);

function mx=freqsmooth(x,w)
L=length(x);mx=zeros(size(x));
for k=1:L
    i1=max(k-w,1);i2=min(k+w,L);
    mx(k)=sum(x(i1:i2))/(i2-i1+1);
end


text:
function homeDirName=GetUserHomeDirName()

if isunix() || ismac()
  homeDirName=getenv('HOME');
else
  % Windows
  homeDriveName=getenv('%HOMEDRIVE%');
  homeDirName=getenv('%HOMEDIR%');
  homeDirName=fullfile(homeDriveName,homeDirName);
end

end

author: kristinbranson

repo: JAABA


text:
function result = polygon_1 ( n, v )

%*****************************************************************************80
%
%% POLYGON_INTEGRAL_1 integrates the function 1 over a polygon.
%
%  Discussion
%
%    The polygon is bounded by the points (X(1:N), Y(1:N)).
%
%    INTEGRAL = 0.5 * SUM ( 1 <= I <= N )
%      ( X(I) + X(I-1) ) * ( Y(I) - Y(I-1) )
%
%    where X(0) and Y(0) should be replaced by X(N) and Y(N).
%
%    Note that the integral of 1 over a polygon is the area of the polygon.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    S F Bockman,
%    Generalizing the Formula for Areas of Polygons to Moments,
%    American Mathematical Society Monthly,
%    1989, pages 131-132.
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%    N should be at least 3 for a nonzero result.
%
%    Input, real V(2,N), the coordinates of the vertices
%    of the polygon.  These vertices should be given in counter-clockwise order.
%
%    Output, real RESULT, the value of the integral.
%
  result = 0.0;

  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'POLYGON_INTEGRAL_1 - Warning!\n' );
    fprintf ( 1, '  The number of vertices must be at least 3.\n' );
    fprintf ( 1, '  The input value of N = %d\n', n );
    error ( 'POLYGON_INTEGRAL_1 - Fatal error!' );
  end

  for i = 1 : n

    if ( i == 1 )
      im1 = n;
    else
      im1 = i - 1;
    end

    result = result + 0.5 * ( v(1,i) + v(1,im1) ) * ( v(2,i) - v(2,im1) );

  end

  return
end


text:
function [v s] = ricefit(mudat, vr)
%RICEFIT Estimate parameters for Rice/Rician distribution from data.
%   Uses simple moment-matching approach, with numerical optimisation.
%   [v s] = ricefit(dat) estimates v and s from samples in dat
%   [v s] = ricefit(mu, vr) estimates v and s from given mean and variance
%
%   R ~ Rice(v, s) if R = sqrt(X^2 + Y^2), where X ~ N(v*cos(a), s^2) and
%   Y ~ N(v*sin(a), s^2) are independent normal distributions (any real a).
%   Note that v and s are *not* the mean and standard deviation of R!
%
%   Reference: http://en.wikipedia.org/wiki/Rice_distribution (!)
%
%   Example:
%     % Sample data, fit model, compare expected, fitted & observed moments
%     V = 5; S = 4; N = 1000;
%     [MN VR] = ricestat(V, S)      % expected mean and variance
%     r = ricernd(V*ones(1, N), S);
%     [v s] = ricefit(dat)          % fitted Rician distribution parameters
%     [mn vr] = ricestat(v, s)      % fitted mean and variance
%     mean(dat), var(dat)           % observed mean and variance
%
%   See also RICEPDF, RICERND, RICESTAT

%   Missing (?) 'See also's RICECDF, RICEINV, RICELIKE

%   Inspired by normfit from the MATLAB statistics toolbox
%   Copyright 2008 Ged Ridgway (Ged at cantab dot net)

if nargin == 1
    dat = mudat;
    mu = mean(dat(:));
    vr = var(dat(:));
else
    mu = mudat;
end

% Optimise cost based on difference of mu and vr from ricestat(v, s) as a
% function of v and s, using mu and sqrt(vr) as initial values for v and s
cost = @(x) moment_cost(x(1), x(2), mu, vr);
x = fminsearch(cost, [mu sqrt(vr)]);
if nargout == 1
    v = x;
else
    v = x(1);
    s = x(2);
end

function cost = moment_cost(v, s, MU, VR)
% Very naive approach...
% Not clear how to wait relative errors in mean and variance...
[mu vr] = ricestat(v, s);
cost = norm([mu - MU; vr - VR;]);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% function [res loc] = maxkmex(list, k)
%
% Matlab C-Mex
% Purpose: Same as MAXK, i.e.,
% Return in RES the K largest elements of LIST
% LOC is Location of the largest values: RES=LIST(LOC)
% This MEX works on double only, and output RES is unsorted
% Algorithm according to http://en.wikipedia.org/wiki/Selection_algorithm
% Compilation: mex -O -v maxkmex.c
% Author Bruno Luong <brunoluong@yahoo.com>
% Last update: 07/April/2009
%

error('Mex file not yet compiled. Action: mex -O -v maxkmex.c');

author: andrewssobral

repo: lrslibrary


text:
function [PDAGs, nodes] = mk_nbrs_of_pdag_add(cpdag,engine)
% MK_NBRS_OF_PDAG_ADD Make the superior inclusion boundary of CPDAG.
% [PDAGs, nodes] = mk_nbrs_of_pdag_add(CPDAG)
%
% PDAGs{i} is the i'th neighbor of CPDAG0 generated by INSERT(X,Y,T)  with
% nodes{i,1:2}=[X Y]
% nodes{i,3}=T
%
% See D.M. Chickering 2002 : "Optimal Structure Identification with Greedy Search".
%
% philippe.leray@univ-nantes.fr, francois.olivier.c.h@gmail.com
% 24 july 2003

compteur=0 ;
N=length(cpdag);
G=pdag_to_dag(cpdag);
if nargin==1,
 bnet_tmp=mk_bnet(pdag_to_dag(cpdag),2*ones(N,1));
 bnet_tmp=mk_bnet(G,2*ones(N,1));
 engine_tmp=struct(jtree_inf_engine(bnet_tmp));
 clear bnet_tmp
else
 engine_tmp=struct(engine);
end
cliques=engine_tmp.cliques;
nbcliques=length(cliques);
clear engine_tmp;

% find in the PDAG all the X Y not connected
[LX LY]=find((cpdag|cpdag')+eye(N)==0);
nlinks=length(LX); % here is a bug when nlinks is zeros, for i=1:nlinks fail 31-37 added by hanbin
if nlinks==0
   [i,j]=find(cpdag);k=unidrnd(size(i,1),1);
   cpdag(i(k),j(k))=0;
   [LX LY]=find((cpdag|cpdag')+eye(N)==0);
    nlinks=length(LX);
end

for i=1:nlinks
    X=LX(i);
    Y=LY(i);
    % Neighbors of Y
    NY = myintersect(find(cpdag(:,Y)), find(cpdag(Y,:)));
    % Adjacents of X
    AX =  myunion(find(cpdag(:,X)), find(cpdag(X,:)));
    % Neighbors of Y adjacent to X
    NAYX = myintersect(NY,AX);
    % Neighbors of Y NOT adjacent to X
    NNAYX = mysetdiff(NY, NAYX);

    % this function recursively "walks" (dfs) in the graph representation of NNA powerset
    liste=NNAYX;
    if ~isempty(liste)
        premier=liste(1);
        dernier=liste(end);
    end;

    current_set=[];
    fini=0;
    evite2 = 0 ;

    while ~fini
        isclique=0;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Test 1
        % is [NAYX current_set] a clique ?
        if isempty(current_set)
            NAYXT=NAYX;
        elseif isempty(NAYX)
            NAYXT=current_set;
        else
            NAYXT = union(NAYX,current_set);
        end
        % if isempty(NAYXT),isclique=1;end

        isclique = isempty(NAYXT) | ismemberclique(NAYXT,cliques) ;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End Test 1

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Test 2
        % is there exist a partially directed path Y...X in PDAG \ NAYXT ?
        if isclique
            if evite2
                test2=1 ;
            else
                %%%% calcul test 2
                L2 = setdiff(1:N,NAYXT);
                test2=~partialconnected(cpdag(L2,L2),find(L2==Y),find(L2==X));

            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End Test 2

            if test2
                evitet2=1;

                % avoid testing INSERT(X,Y,[]) and INSERT(Y,X,[]) if pa(Y)=pa(X)
                test0=1 ;
                if (X>Y)
                    if length(current_set)==0
                        PaX=setdiff(find(cpdag(X,:)),find(cpdag(:,X)));
                        PaY=setdiff(find(cpdag(Y,:)),find(cpdag(:,Y)));
                        test0=(length(setdiff(PaX,PaY))~=0);
                    end
                end
                if test0
                    %fprintf(' INSERT(%d,%d,',X,Y); fprintf('%d',current_set);
                    %fprintf(')\n');
                    compteur=compteur+1;
                    nodes{compteur,1}=X;
                    nodes{compteur,2}=Y;
                    nodes{compteur,3}=current_set;
                    ptmp=cpdag;
                    ptmp(X,Y)=1;
                    ptmp(current_set,Y)=1;
                    ptmp(Y,current_set)=0;
                    PDAGs{compteur}=ptmp;
                end
            end
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Next Set to test ...
        % what is the net set in the powerset ?
        if length(liste)==0
            fini=1 ; next_set=[];
        elseif length(current_set)==0
            next_set=[premier];     % first node after the root []
        else
            actuel=current_set(end);
            if actuel==dernier
                if length(current_set)==1
                    fini=1;            % no more node ...
                else
                    ancien=current_set(end-1);  % new "branch"
                    next_set=[current_set(1:end-2) liste(find(liste==ancien)+1)];
                end
            else                            % new node in the "branch"
                if ~isclique
                    if length(current_set)==1
                        fini=1;            % no more node ...
                    else
                        ancien=current_set(end-1);  % new "branch"
                        next_set=[current_set(1:end-2) liste(find(liste==ancien)+1)];
                        evite2=0;
                    end
                else
                    next_set=[current_set liste(find(liste==actuel)+1)];
                end
            end
        end
        current_set=next_set;
    end
end

%%%%%%%%%
function resu = ismemberclique(v,cliques)
finiclique = 0 ; resu=0 ;
cl=1;  ncl=length(cliques) ;
while (~finiclique) & (cl<=ncl);
    if ismember(v,cliques{cl})
        resu=1;
        finiclique=1 ;
    end
    cl=cl+1;
end

%%%%%%%%%
function resu = partialconnected(G,Y,X)
tmp=expm(G);
resu = (tmp(Y,X)~=0);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function imchi=sskkimbook2(omega,rechi,omega1,imchi1)
%The program inputs are the vector of the frequency
%(or energy) components, the vector of the real
% part of the susceptibility
%under examination, anchor point, the value of the
%imaginary part at the anchor point.
%The two vectors must have the same length
%and the frequency vector omega must be equispaced.
%If not, apply MATLAB functions such as interp.
%Note that the anchor point must be situated in one of the values
%of the vector omega.
%This function uses the function kkimbook2
%The output is the estimate of the
%real part as obtained by using SSKK relations.
%This software is distributed under the GNU licence agreement
%by Valerio Lucarini
%email: lucarini@alum.mit.edu
%University of Camerino
%Department of Mathematics and Computer Science
%Camerino, Italy

if size(omega,1)>size(omega,2);
omega=omega';
end; if size(rechi,1)>size(rechi,2);
rechi=rechi';
end;
%Here the program rearranges the two vectors so that,
%whichever their initial shape, they become row vectors.
g=size(omega,2);
%Size of the vectors.%
k=0;
for j=1:g;
if omega(j)==omega1;
k=j;
end;
end;
%Determination of the anchor point.
imchi=kkimbook2(omega,rechi,0);
%Application of K-K relations
imchi=imchi+omega1^(-1)*omega.^(1)*(imchi1-imchi(k));
%The subtracted relation upgrades the estimate obtained
%with K-K relations.


text:
function [varargout]=voxelDistMarch(varargin)

%% Parse input

%Default option structure
defaultOptionStruct.toleranceLevel=0;
defaultOptionStruct.waitBarOn=false(1,1);
defaultOptionStruct.unitEdgeOn=false(1,1);
defaultOptionStruct.W=[];
defaultOptionStruct.Ws=[];

switch nargin
    case 1
        M=varargin{1};
        indStart=1;
        voxelSize=[1 1 1];
        optionStruct=defaultOptionStruct;
    case 2
        M=varargin{1};
        indStart=varargin{2};
        voxelSize=[1 1 1];
        optionStruct=defaultOptionStruct;
    case 3
        M=varargin{1};
        indStart=varargin{2};
        voxelSize=varargin{3};
        optionStruct=defaultOptionStruct;
    case 4
        M=varargin{1};
        indStart=varargin{2};
        voxelSize=varargin{3};
        optionStruct=varargin{4};
end

if numel(voxelSize)==1
    voxelSize=voxelSize*ones(1,3);
end
defaultOptionStruct.numSeeds=numel(indStart);

[optionStruct]=structComplete(optionStruct,defaultOptionStruct,1);

if nargout==2
    computeSeedIndex=true(1,1);
else
    computeSeedIndex=false(1,1);
end

% Get optionional inputs
toleranceLevel=optionStruct.toleranceLevel;
numSeeds=optionStruct.numSeeds;
waitBarOn=optionStruct.waitBarOn;
unitEdgeOn=optionStruct.unitEdgeOn;

numVoxels=numel(M); %Number of vertices

W=optionStruct.W;
if isempty(W)
    W=ones(numVoxels,1);
end

Ws=optionStruct.Ws;
if isempty(Ws)
    Ws=ones(numel(indStart),1);
end

numStart=numel(indStart);
if numSeeds<numStart
    warning('Number of seeds is smaller than number of start points, assuming these should be additional points')
    numSeeds=numSeeds+numStart;
    optionStruct.numSeeds=numSeeds; %override input in case function is called recursively
end
numSteps=(numSeeds-numStart)+1;

%% Variables and computations outside of iterative loop

if computeSeedIndex
    indAll=(1:1:numVoxels)'; %Row indices for all points
end

% k=sphereIm(1);
% ind=find(k);
% [ii,jj,kk]=ind2sub(size(k),ind);
% ii=ii([1 2 3 5 6 7]);
% jj=jj([1 2 3 5 6 7]);
% kk=kk([1 2 3 5 6 7]);
% ii=ii-2;
% jj=jj-2;
% kk=kk-2;

k=ones(3,3,3);
ind=find(k);
[ii,jj,kk]=ind2sub(size(k),ind);
ii=ii([1:13 15:end]);
jj=jj([1:13 15:end]);
kk=kk([1:13 15:end]);
ii=ii-2;
jj=jj-2;
kk=kk-2;

dd=sqrt((ii.*voxelSize(1)).^2+(jj.*voxelSize(2)).^2+(kk.*voxelSize(3)).^2);

C=maskfind(M,1:numel(M),ii,jj,kk);
siz=size(C); %Size of C
L=C>0; %Logic for valid indices in C

%Calculate edge lengths
if unitEdgeOn
    DE=nan(siz);
    DE(L)=1; %Unit edge lenghts
else
    DE=repmat(dd,[1 numel(M)])';
    DE(~L)=0;
end

%Initialize distance vector
d=nan(numVoxels,1);
d(indStart)=0;
DD=nan(siz);
IND_d=C(L);

%Initialize seed index vector
if computeSeedIndex
    i=nan(numVoxels,1);
    i(indStart)=indStart;
    II=nan(siz);
end

d_previous=d;
i_previous=i;

%% Propagate and iteratively update distances
indSeed=nan(1,numSeeds);
indSeed(1:numStart)=indStart;

if numSteps>1
    numLoopSteps=numSteps+1;
else
    numLoopSteps=1;
end

if waitBarOn
    hw=waitbar(1/numLoopSteps,['meshDistMarch...',num2str(round(100.*1/numLoopSteps)),'%']);
end

for q=1:1:numLoopSteps
    indStart=indSeed(~isnan(indSeed));
    while 1

        %Update distance data
        DD(L)=d(IND_d); %The distance data currently at neighbours
        if isempty(Ws)
            D_check=DD+DE; %Distance data plus edge lenghts
        else
            WS=ones(numVoxels,1);
            WS(indStart)=Ws;
            WS(~isnan(i))=WS(i(~isnan(i)));
            D_check=DD+DE./WS(:,ones(size(DE,2),1)); %Distance data plus edge lenghts
        end

        if computeSeedIndex
            [d,indMin]=min(D_check,[],2,'omitnan'); %Assign minimum distance

            %Update index data
            II(L)=i(IND_d); %The seed indices currently at neighbours
            ind=sub2ind(siz,indAll,indMin); %Linear indices for minimum used
            i=II(ind); %Get seed index at point chosen
            i(indStart)=indStart; %Override start seed indices
        else
            d=min(D_check,[],2,'omitnan'); %Assign minimum distance
        end
        d(indStart)=0; %Override starts to be zero

        logicFix=d_previous<d;
        d(logicFix)=d_previous(logicFix);
        if computeSeedIndex
            i(logicFix)=i_previous(logicFix);
        end

%         sum((d_previous-d).^2,'omitnan')
%         [~,indMax]=nanmax((d_previous-d).^2)

        %Check convergence
        if nnz(isnan(d))==nnz(isnan(d_previous)) %Check if number of NaN's changed
            if sum((d_previous-d).^2,'omitnan')<=toleranceLevel %Check if converged to within tolerance
                break %Stop while loop
            end
        end

        d_previous=d; %Store previous
        if computeSeedIndex
            i_previous=i;
        end

    end
    if q>1 && q<=numSteps
        [~,indSeed(numStart+(q-1))]=max(d.*W,[],'omitnan');
    end

    if waitBarOn
        waitbar(q/numLoopSteps,hw,['meshDistMarch...',num2str(round(100.*q/numLoopSteps)),'%']);
    end
end

if waitBarOn
    close(hw);
end

%% Store output
varargout{1}=reshape(d,size(M));
if nargout==2
    varargout{2}=reshape(i,size(M));
end

end
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% 预测世界石油产量
%% 历史数据
clear;clc;
x = (1994 : 2003)';    % 年份
year = 1990;
x = x - year;
y = [67.052, 68.008, 69.803, 72.024, 73.400, 72.063, 74.669, 74.487, 74.065, 76.777]';   % 百万桶每日
plot(x,y,'rp');
%% 直线拟合
C = [x.^0, x.^1];
A = C' * C;  % 法方程系数矩阵
b = C' * y;  % 法方程向量
p_1 = A \ b   % 一次多项式系数，升幂排列
RMSE_1 = sqrt(sum((y - C * p_1).^2)) / sqrt(length(C))
%% 抛物线拟合
C = [x.^0, x.^1, x.^2];
A = C' * C;  % 法方程系数矩阵
b = C' * y;  % 法方程向量
p_2 = A \ b   % 二次多项式系数，升幂排列
RMSE_2 = sqrt(sum((y - C * p_2).^2)) / sqrt(length(C))
%% 三次曲线拟合
C = [x.^0, x.^1, x.^2, x.^3];
A = C' * C;  % 法方程系数矩阵
b = C' * y;  % 法方程向量
p_3 = A \ b   % 三次多项式系数，升幂排列
RMSE_3 = sqrt(sum((y - C * p_3).^2)) / sqrt(length(C))
%% 开始预测
x1=x(end) + (1:5)';   % 待预测时间, 10年
xi=[x;x1];
p = p_1;   % 指定模型
plot(x + year,y,'*',xi + year,polyval(flip(p),xi),x1 + year,polyval(flip(p),x1),'rp');
% 结论：预测结果仅供参考

repo: NumericalAnalysis


text:
function jrt=jacprojRT(j, i, rt, xyz, r0, a)
% symbolic projection function Jacobian
% code automatically generated with maple

  qr0=r0(j*4+1:(j+1)*4);

  t1 = (rt(1) ^ 2);
  t2 = (rt(2) ^ 2);
  t3 = (rt(3) ^ 2);
  t5 = sqrt((1 - t1 - t2 - t3));
  t6 = 0.1e1 / t5;
  t7 = t6 * qr0(2);
  t9 = -t7 * rt(1) + qr0(1);
  t11 = t6 * qr0(3);
  t13 = -t11 * rt(1) - qr0(4);
  t15 = t6 * qr0(4);
  t17 = -t15 * rt(1) + qr0(3);
  t19 = -t9 * xyz(1) - t13 * xyz(2) - t17 * xyz(3);
  t24 = -t5 * qr0(2) - qr0(1) * rt(1) - rt(2) * qr0(4) + rt(3) * qr0(3);
  t31 = t5 * qr0(3) + qr0(1) * rt(2) + rt(3) * qr0(2) - rt(1) * qr0(4);
  t37 = t5 * qr0(4) + qr0(1) * rt(3) + rt(1) * qr0(3) - rt(2) * qr0(2);
  t39 = t24 * xyz(1) - t31 * xyz(2) - t37 * xyz(3);
  t41 = t6 * qr0(1);
  t43 = -t41 * rt(1) - qr0(2);
  t48 = t5 * qr0(1) - rt(1) * qr0(2) - rt(2) * qr0(3) - rt(3) * qr0(4);
  t52 = t48 * xyz(1) + t31 * xyz(3) - t37 * xyz(2);
  t57 = t43 * xyz(1) + t13 * xyz(3) - t17 * xyz(2);
  t62 = t43 * xyz(2) + t17 * xyz(1) - t9 * xyz(3);
  t67 = t48 * xyz(2) + t37 * xyz(1) + t24 * xyz(3);
  t72 = t43 * xyz(3) + t9 * xyz(2) - t13 * xyz(1);
  t77 = t48 * xyz(3) - t24 * xyz(2) - t31 * xyz(1);
  t89 = -t19 * t31 - t39 * t13 + t43 * t67 + t48 * t62 + t72 * t24 - t77 * t9 + t57 * t37 + t52 * t17;
  t99 = -t19 * t37 - t39 * t17 + t43 * t77 + t48 * t72 - t57 * t31 - t52 * t13 - t62 * t24 + t67 * t9;
  t106 = -t39 * t37 + t48 * t77 - t52 * t31 - t67 * t24 + rt(6);
  t107 = 0.1e1 / t106;
  t119 = -t39 * t31 + t48 * t67 + t77 * t24 + t52 * t37 + rt(5);
  t123 = t106 ^ 2;
  t124 = 0.1e1 / t123;
  t125 = (a(1) * (t39 * t24 + t48 * t52 - t67 * t37 + t77 * t31 + rt(4)) + a(2) * t119 + a(3) * t106) * t124;
  t129 = -t7 * rt(2) + qr0(4);
  t132 = -t11 * rt(2) + qr0(1);
  t135 = -t15 * rt(2) - qr0(2);
  t137 = -t129 * xyz(1) - t132 * xyz(2) - t135 * xyz(3);
  t141 = -t41 * rt(2) - qr0(3);
  t146 = t141 * xyz(1) + t132 * xyz(3) - t135 * xyz(2);
  t151 = t141 * xyz(2) + t135 * xyz(1) - t129 * xyz(3);
  t157 = t141 * xyz(3) + t129 * xyz(2) - t132 * xyz(1);
  t170 = -t137 * t31 - t39 * t132 + t141 * t67 + t48 * t151 + t157 * t24 - t77 * t129 + t146 * t37 + t52 * t135;
  t180 = -t137 * t37 - t39 * t135 + t141 * t77 + t48 * t157 - t146 * t31 - t52 * t132 - t151 * t24 + t67 * t129;
  t187 = -t7 * rt(3) - qr0(3);
  t190 = -t11 * rt(3) + qr0(2);
  t193 = -t15 * rt(3) + qr0(1);
  t195 = -t187 * xyz(1) - t190 * xyz(2) - t193 * xyz(3);
  t199 = -t41 * rt(3) - qr0(4);
  t204 = t199 * xyz(1) + t190 * xyz(3) - t193 * xyz(2);
  t209 = t199 * xyz(2) + t193 * xyz(1) - t187 * xyz(3);
  t215 = t199 * xyz(3) + t187 * xyz(2) - t190 * xyz(1);
  t228 = -t195 * t31 - t39 * t190 + t199 * t67 + t48 * t209 + t215 * t24 - t77 * t187 + t204 * t37 + t52 * t193;
  t238 = -t195 * t37 - t39 * t193 + t199 * t77 + t48 * t215 - t204 * t31 - t52 * t190 - t209 * t24 + t67 * t187;
  t255 = (a(4) * t119 + a(5) * t106) * t124;
  jrt(1) = (a(1) * (t19 * t24 - t39 * t9 + t43 * t52 + t48 * t57 - t62 * t37 - t67 * t17 + t72 * t31 + t77 * t13) + a(2) * t89 + a(3) * t99) * t107 - t125 * t99;
  jrt(2) = (a(1) * (t137 * t24 - t39 * t129 + t141 * t52 + t48 * t146 - t151 * t37 - t67 * t135 + t157 * t31 + t77 * t132) + a(2) * t170 + a(3) * t180) * t107 - t125 * t180;
  jrt(3) = (a(1) * (t195 * t24 - t39 * t187 + t199 * t52 + t48 * t204 - t209 * t37 - t67 * t193 + t215 * t31 + t77 * t190) + a(2) * t228 + a(3) * t238) * t107 - t125 * t238;
  jrt(4) = a(1) * t107;
  jrt(5) = a(2) * t107;
  jrt(6) = a(3) * t107 - t125;
  jrt(7) = (a(4) * t89 + a(5) * t99) * t107 - t255 * t99;
  jrt(8) = (a(4) * t170 + a(5) * t180) * t107 - t255 * t180;
  jrt(9) = (a(4) * t228 + a(5) * t238) * t107 - t255 * t238;
  jrt(10) = 0.0e0;
  jrt(11) = a(4) * t107;
  jrt(12) = a(5) * t107 - t255;

author: vrabaud


text:
function [Q, R, E] = qr(f, varargin)
%QR   QR factorisation of an array-valued BNDFUN.
%   [Q, R] = QR(F) returns a QR factorisation of F such that F = Q*R, where the
%   BNDFUN Q is orthogonal (with respect to the continuous L^2 norm on the
%   domain of F) and of the same size as F and R is an m x m upper-triangular
%   matrix when F has m columns.
%
%   [Q, R, E] = QR(F) produces unitary Q, upper-triangular R, and a permutation
%   matrix E so that F*E = Q*R. The column permutation E is chosen to reduce
%   fill-in in R.
%
%   [Q, R, E] = QR(F, 'vector') returns the permutation information as a vector
%   instead of a matrix.  That is, E is a row vector such that F(:,E) = Q*R.
%   Similarly, [Q, R, E] = QR(F, 'matrix') returns a permutation matrix E. This
%   is the default behavior.

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

% Deal with empty case:
if ( isempty(f) )
    Q = [];
    R = [];
    E = [];
    return
end

% Initialise Q to be a BNDFUN:
Q = f;

% Rescaling factor, (b - a)/2.
rescaleFactor = .5*diff(f.domain);

% Call QR on the ONEFUN of f:
if ( nargout == 3 )
    [Q.onefun, R, E] = qr(f.onefun, varargin{:});
else
    [Q.onefun, R] = qr(f.onefun, varargin{:});
end

% Rescale so that columns of Q will be orthonormal (rather than orthogonal):
Q = Q/sqrt(rescaleFactor);

% Rescale R so that f = QR:
R = R*sqrt(rescaleFactor);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
cur_dir = pwd;
cd(fileparts(mfilename('fullpath')));

try
    fprintf('Downloading model_ResNet-101L...\n');
    urlwrite('https://onedrive.live.com/download?resid=F371D9563727B96F!91963&authkey=!AM-EuzuUJelv9Po', ...
        'models_ResNet-101L.zip');

    fprintf('Unzipping...\n');
    unzip('models_ResNet-101L.zip', '..');

    fprintf('Done.\n');
    delete('models_ResNet-101L.zip');
catch
    fprintf('Error in downloading, please try links in README.md https://github.com/daijifeng001/R-FCN');
end

cd(cur_dir);

repo: R-FCN


text:
% foo.m
%
% This m-file shows an example of how to use ApEn function
%
% More specifically, it generates three simulated data with different
% complexity (sine, chirp, and Gaussian noise), and plots the ApEn values
% with varying r

%-------------------------------------------------------------------------
% coded by Kijoon Lee, kjlee@ntu.edu.sg
% Mar 21st, 2012
%-------------------------------------------------------------------------
clear all

m = 2;      % embedded dimension
tau = 1;    % time delay for downsampling
N = 1000;
t = 0.001*(1:N);

% generate simulated data
sint = sin(2*pi*10*t);      % sine curve
chirpt = chirp(t,0,1,150);  % chirp signal
whitet = randn(1,N);        % white noise

% calculate standard deviations
sd1 = std(sint);
sd2 = std(chirpt);
sd3 = std(whitet);

% specify the range of r
rnum = 30;
result = zeros(3,rnum);

% main calculation and display
figure
for i = 1:rnum
    r = i*0.02;
    result(1,i) = ApEn(m, r*sd1, sint, tau);
    result(2,i) = ApEn(m,r*sd2,chirpt, tau);
    result(3,i) = ApEn(m,r*sd3,whitet, tau);
end

r = 0.02*(1:rnum);
plot(r,result(1,:),'o-',r,result(2,:),'o-',r,result(3,:),'o-')
axis([0 rnum*0.02 0 1.05*max(result(:))])
legend('sin','chirp','white noise')
title(['ApEn, m=' num2str(m) ', \tau=' num2str(tau)],'fontsize',14)
xlabel r


text:
function [h,p,ci,stats, p_square, t_square, q_square, OUT] = ttest(bs, wh_group, varargin)
%TTEST Conducts a one- or two-sampled ttest between subjects at every node
%
% FDR correction: positive false discovery rate (pFDR) from
%   the p-values P of multiple-hypothesis testing using the procedure
%   described by Storey (2002).
%
%   wh_group: a logical array, n_subjects x 1, that defines the two groups.
%             Results show wh_group == T > wh_group == F
%             If empty, performs one-sample t-test
%
%   n_subjects must be equal to size(bs.connectivity.regions.r, 3). For
%   now, this function defaults to testing the matrices in
%   bs.connectivity.regions.r, but could be expanded later.
%
%   Input:
%       'doplot': show a plot
%       'labels': followed by a 2-element cell array, containing labels for
%                 the top two plots
%
%   Output:
%       [h,p,ci,stats] from matlab's ttest2
%       p_square, t_square: p and t values, in the square format
%       q_square: fdr corrected p-values, from mafdr()
%       OUT: structure of summary stats compatible with plot_correlation_matrix.m
%
% examples:
%
%  [h,p,ci,stats, p_square, t_square, q_square, OUT] = ttest(obj, []);
%  plot_correlation_matrix(OUT, 'nofigure', varargin{:});
%  title('Mean connectivity, q < 0.05 pFDR thresholded');
%
%   Yoni Ashar, March 2020

% Programmers' notes:
% Tor Wager updated Dec 2020 to add one-sample t-test

% vectorize. mat is then subjects x nodes in lower triangle

if isempty(which('mafdr')), error('Sorry, you need the Matlab bioinformatics toolbox on your Matlab path to use the Storey 2002 mafdr function called here.'); end

if nargin < 2
    twosample = false;
else
    twosample = ~isempty(wh_group);
end

if twosample && ~islogical(wh_group), error('wh_group must be a logical array'), end

doplot = 0;
labels = {'Mean r, Group == 1', 'Mean r, Group==0'};

for i=1:length(varargin)

    if ischar(varargin{i})
        switch varargin{i}
            case 'doplot'
                doplot = 1;
            case 'labels'
                labels = varargin{i+1};

        end
    end
end

mat = bs.flatten_conn_matrices();

% test at each node

if twosample
    % -- ttest2 tests columnwise
    [h,p,ci,stats] = ttest2(mat(wh_group, :), mat(~wh_group, :));
else
    [h,p,ci,stats] = ttest(mat);
end

% revert back to square form
t_square = squareform(stats.tstat, 'tomatrix');

p_square = squareform (p, 'tomatrix');
p_square(logical(eye(size(p_square)))) = 1; % set the diagonal p values to be 1

[fdr] = mafdr(p');%, 'BHFDR', true);
q_square = squareform (fdr, 'tomatrix');
q_square(logical(eye(size(q_square)))) = 1; % set the diagonal p values to be 1

if twosample

    OUT = struct('descrip', 'Two-sample t-test; r is connectivity diff group 1 - 2; p field is pFDR q-values', ...
        'r', mean(bs.connectivity.regions.r(:,:,wh_group), 3) - mean(bs.connectivity.regions.r(:,:, ~wh_group), 3), ...
        'p', q_square, 'sig', q_square < 0.05, 'wh_group', wh_group);

else

    OUT = struct('descrip', 'One-sample t-test; r is mean correlation, p field is pFDR q-values', 'r', mean(bs.connectivity.regions.r, 3), 'p', q_square, 'sig', q_square < 0.05);

end

if doplot

    if twosample

        create_figure('brainpathways multi two-sample t-test', 2, 2)
        mat = bs.connectivity.regions.r;
        imagesc(mean(mat(:,:,wh_group), 3)); colorbar, title(labels{1})
        subplot(2,2,2)
        imagesc(mean(mat(:,:,~wh_group), 3)); colorbar, title(labels{2})
        subplot(2,2,3)
        imagesc(t_square); colorbar, title('T-statistic: group difference')
        subplot(2,2,4)
        imagesc(q_square); colorbar, title('FDR corrected p values: group difference')

    else

        plot_correlation_matrix(OUT, 'nofigure', varargin{:});

        title('Mean connectivity, q < 0.05 pFDR thresholded');

%         create_figure('brainpathways multi one-sample t-test', 1, 3)
%         mat = bs.connectivity.regions.r;
%         m = mean(mat, 3);
%
%         imagesc(m);
%         colorbar, title('Mean correlation')
%         subplot(1, 3, 2)
%
%         imagesc(t_square); colorbar, title('T-statistic: group')
%
%         subplot(1, 3, 3)
%         mt = m;
%         mt(q_square >= 0.05) = 0;
%         imagesc(mt); colorbar, title('FDR corrected p values: group difference')

        drawnow
    end
end

end


text:
function pass = test_innerProduct(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end

singPref = pref;
singPref.blowup = true;

% Seed for random number:
seedRNG(6178);

%% Functions on [-inf inf]:

% Set the domain:
dom = [-Inf Inf];

opf = @(x) 2-exp(-x.^2);
opg = @(x) exp(-x.^2);
f = unbndfun(opf, struct('domain', dom));
g = unbndfun(opg, struct('domain', dom));

I = innerProduct(f, g);
IExact = (sqrt(pi)*(4 - sqrt(2)))/2;
err = abs(I - IExact);
pass(1) = err < 2e7*max(eps*get(f,'vscale'), ...
    eps*get(g,'vscale'));

%% Functions on [a inf]:

% Set the domain:
dom = [1 Inf];

opf = @(x) x;
opg = @(x) exp(-x);
pref = chebfunpref();
f = unbndfun(opf, struct('domain', dom, 'exponents', [0 1]), pref);
g = unbndfun(opg, struct('domain', dom));
warning('off', 'CHEBFUN:UNBNDFUN:sum:slowDecay');
I = innerProduct(f, g);
warning('off', 'CHEBFUN:UNBNDFUN:sum:slowDecay');
IExact = 2*exp(-1);
err = abs(I - IExact);
pass(2) = err < 2e8*max(eps*get(f,'vscale'), ...
    eps*get(g,'vscale'));

%% Functions on [-inf b]:

% Set the domain:
dom = [-Inf -3*pi];

opf = @(x) 1./x;
opg = @(x) 2./x;
f = unbndfun(opf, struct('domain', dom));
g = unbndfun(opg, struct('domain', dom));
I = innerProduct(f, g);
IExact = 2/(3*pi);
err = abs(I - IExact);
pass(3) = err < 1e5*eps*get(f,'vscale');

end


text:
%TESTFEATDOM Test feature domains
%
%	   [N,I,J] = TESTFEATDOM(A,K,M)
%
% The feature domains of the dataset A are tested. When given,
% the vector K should contain the indices of the features to
% be tested. Optionally, a vector M can be given which indicates which
% of the objects should be taken for testing.
%
% N = 0 if the dataset values are within the domains to be tested.
% N = 1 if somewhere a value is outside a domain. I and J return
%       the indices to the erroneous objects and features.
%
%	TESTFEATDOM(A,K,M)
%
% Performs the test and prints an error message if a feature value
% outside a domain is encountered.

author: marianux

repo: ecg-kit

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


text:
function chanMaps = loadChanMaps()

ksroot = fileparts(fileparts(mfilename('fullpath')));
chanMapFiles = dir(fullfile(ksroot, 'configFiles', '*.mat'));

idx = 1;
chanMaps = [];
orig_state = warning;
warning('off','all')

for c = 1:numel(chanMapFiles)

    q = load(fullfile(ksroot, 'configFiles', chanMapFiles(c).name));

    cm = createValidChanMap(q, chanMapFiles(c).name);
    if ~isempty(cm)
        if idx==1; chanMaps = cm; else; chanMaps(idx) = cm; end
        idx = idx+1;
    end

end

warning(orig_state);

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function coeffs = alias(coeffs, m)
%ALIAS   Alias Chebyshev coefficients on the 1st kind Chebyshev grid.
%   ALIAS(C, M) aliases the Chebyshev coefficients stored in the column vector C
%   to have length M. If M > LENGTH(C), the coefficients are padded with zeros.
%   If C is a matrix of coefficients, each of the columns is aliased to length
%   M.
%
% See also PROLONG.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note that the formula for aliasing on the 1st-kind Chebyshev grid is
% different from that for the 2nd-kind grid, even though the coefficients
% being aliased are for 1st-kind Chebyshev polynomials in both cases.
%
% Useful References:
%   Fox, L. and Parker, I. B., Chebyshev polynomials in Numerical Analysis,
%   Oxford University Press, 1972.  (pp. 67)
%
%   Mason, J. C. and Handscomb, D. C., Chebyshev polynomials, Chapman &
%   Hall/CRC, Boca Raton, FL, 2003.  (pp. 153)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = size(coeffs, 1);

% Pad with zeros:
if ( m > n )
    coeffs = [ coeffs ; zeros(m-n, size(coeffs, 2)) ];
    return
end

% Alias coefficients:
if ( m == 1 )
    % Reduce to a single point:
    e = ones(1, ceil(n/2));
    e(2:2:end) = -1;
    coeffs = e*coeffs(1:2:end,:);
elseif ( m > n/2 )
    % If m > n/2, only single coefficients are aliased, and we can vectorise.
    j = ((m + 1):n).';
    k = abs(mod(j + m - 2, 2*m) - m + 1) + 1;
    p = floor((j-1+m)/(2*m));
    t = (-1).^p;
    coeffs(k,:) = coeffs(k,:) + bsxfun(@times, t, coeffs(j,:));
else
    % Otherwise we must do everything in a tight loop. (Which is slower!)
    for j = (m + 1):n
        k = abs(mod(j + m - 2, 2*m) - m + 1) + 1;
        sgn = 1 - 2*mod(floor((j - 1 + m)/(2*m)), 2);
        coeffs(k,:) = coeffs(k,:) + sgn*coeffs(j,:);
    end
end

% Truncate:
coeffs = coeffs(1:m,:);

end


text:
function z = minus( x, y, cheat )

%   Disciplined convex programming information for MINUS:
%      Terms in a difference must have opposite curvature. Real affine
%      expressions are both convex and concave, so they can be involved
%      in a difference with any nonlinear expression. Complex affine (or
%      constant) expressions, however, are neither, so they can only be
%      involved in differences with other affine expressions. So, for
%      example, the following differences are valid:
%         {convex}-{concave}   {concave}-{convex}   {affine}-{affine}
%      The following are not:
%         {convex}-{concave}  {convex}-{complex constant}
%      For vectors, matrices, and arrays, these rules are verified
%      independently for each element.
%
%   Disciplined geometric programming information for MINUS:
%      Non-constant expressions (log-convex or log-concave) may not be
%      involved in a subtraction in disciplined geometric programs.

if nargin < 3, cheat = false; end
z = plus( x, y, true, cheat );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [retval] = getbatch

% GETBATCH returns a number that can be used to distinguish the subsequent
% batches. The number automatically increases over subsequent calls.

% Copyright (C) 2011-2012, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% 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$

persistent batch
if isempty(batch)
  batch = 0;
end
batch = batch + 1;
retval = batch;

author: fieldtrip

repo: fieldtrip