text:
function [proj, angles, blk, projinfo] = BatchReadXim(datafolder,varargin)
% [proj, angle, blk, projinfo] = BatchReadXim(foldername, varargin)
% Varian .xim files exported from TrueBeam on-board imager/EPID, Halcyon imager
% .xim files are exported under service mode, as a single .zip file
% Once unzipped, all relevant files are organized under the same root folder
% Input:
%       datafolder    : root folder of the upzipped files,
%                       or the subfolder Acquisitions where all xim files are stored
%       thd           : 0(default), no motion lag correcion
%                       thd,  motion correcion threadhold
% Output:
%       proj          : .xim pixel_images frame by frame
%       angles        : KVSourceRtn at each frame
%       projinfo      : proj_*.xim info cell arrays
% Date: 2020-04-11
% Author: Yi Du, yi.du@hotmail.com

thd = 0;
if(nargin == 2)
    thd = varargin{1};
end

read_xim_info = 0;
if(nargout==4)
    read_xim_info = 1;
end

fprintf("Dataset in: %s \n", datafolder);

ximfilelist = dir([datafolder filesep '**' filesep 'Proj_*.xim']);

proj = [];
angles = [];

%
proj_no = length(ximfilelist);
%
count = 1;
for ii = 1:proj_no
	ximfilename = fullfile(ximfilelist(ii).folder, ximfilelist(ii).name);
    if(~mod(ii,50))
        fprintf("Loading: %d / %d \n", ii, proj_no);
    end
	[page, rtn] = mexReadXim(ximfilename);
	if(~isempty(page))
        % load first page
        if(count==1)
            projinfo{count} = ReadXim(ximfilename, 0);
    		angles(count) = rtn;
        	proj(:,:,count) = double(page');
            count = count + 1;
        else
            % proj info
            if(read_xim_info)
                projinfo{count} = ReadXim(ximfilename, 0);
            end

            % whether threshold is applied
            if(thd)
                if(abs(rtn-angles(end))>thd)
                    angles(count) = rtn;
                    proj(:,:,count) = double(page');
                    count = count + 1;
                else
                    continue;
                end
            else
                angles(count) = rtn;
                proj(:,:,count) = double(page');
                count = count + 1;
            end
        end
	end
end
fprintf("Loading complete \n");

%% Load Blk
tag_bowtie = strfind(projinfo{1}.properties.KVCollimatorShape, 'None');
% wi Bowtie1
if(isempty(tag_bowtie))
    blkfilestr = dir([datafolder filesep 'Calibrations' filesep 'AIR-*' filesep '**' filesep 'FilterBowtie.xim']);
else
    blkfilestr = dir([datafolder filesep 'Calibrations' filesep 'AIR-*' filesep '**' filesep 'Filter.xim']);
end

blkfilename = fullfile(blkfilestr.folder, blkfilestr.name);
blk = mexReadXim(blkfilename);
blk = double(blk');

%% Don't know why mAs correction doesn't work
%{
% mAs for projection sampling
mAsproj = projinfo{1}.properties.KVMilliAmperes * projinfo{1}.properties.KVMilliSeconds;
% mAs for blk field
blkinfo = ReadXim(blkfilename, 0);
mAsblk = blkinfo.properties.KVMilliAmperes * blkinfo.properties.KVMilliSeconds;
% mAs output correction
gain = mAsproj/mAsblk;
%}
% Emperical value
gain = 2.5;
blk = blk * gain;
% denoising blk
blk = medfilt2(blk, [5 5]);

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function Plots1D(Params, Model, YTM, ShortRates)
% =========================================================================
% Plots yield curve for estimated model. All rates are annualy compounded.
%
% INPUT:
%
%         Params - Estimated parameters
%          Model - 'NS' for Nelson-Siegel model
%                  'Svensson' for Svensson model
%            YTM -
%         Settle -
% USES:
%        NScurve
%
% Kamil Kladivko
% email: kladivko@gmail.com
% December 2010
% Cite as:
% Kladivko Kamil (2010). The Czech Treasury Yield Curve from 1999 to the Present,
% Czech Journal of Economics and Finance, 60(4): 307-335
% =========================================================================

% Figure 1 (Spot, FWD and Par rates)
FontSize = 14;
Mat = 0:0.1:ceil(max(YTM.TimeFraction));
%Mat = 0:0.1:40; % 27.1.2010
[ZeroRates InstFwdRates] = NScurve(Params, Mat, Model);
ZeroRates = exp(ZeroRates) - 1; % Annual compounding
InstFwdRates = exp(InstFwdRates) - 1; % Annual compounding
% Par Rates
MatForPar = 1:1:ceil(max(YTM.TimeFraction));
[ZeroForParRates junk] = NScurve(Params, MatForPar, Model);
ZeroForParRates = exp(ZeroForParRates) - 1;  % Annual compounding
ParRates = zero2par(ZeroForParRates);
% Plot Zero, Observed YTM, Fit YTM, Par, Fwd
figure
plot(Mat, ZeroRates.*100, '-b', 'LineWidth', 2) % Zero
hold on
set(gca, 'FontSize', FontSize, 'FontName', 'Arial');
plot(MatForPar, ParRates.*100, '--b', 'MarkerSize', 8, 'LineWidth', 2) % Par rates
plot(Mat, InstFwdRates.*100, ':b', 'MarkerSize', 8, 'LineWidth', 2) % FWD rates
if ~isempty(ShortRates)
    SRcount = length(ShortRates.IR);
else
    SRcount = 0;
end
plot(YTM.TimeFraction(SRcount+1:end), YTM.Obs(SRcount+1:end).*100, 'ko', 'MarkerSize', 5, 'MarkerFaceColor','k') % YTM of observed bonds
plot(YTM.TimeFraction(SRcount+1:end), YTM.Fit(SRcount+1:end).*100, 'ks', 'MarkerSize', 7) % YTM of fitted bonds
plot(YTM.TimeFraction(1:SRcount), YTM.Obs(1:SRcount).*100, 'k^', 'MarkerSize', 6, 'MarkerFaceColor','k')
% Bonds ID (issue) number
BondsCount = length(YTM.TimeFraction);
minYTM = min(YTM.Obs, YTM.Fit);
j = 1;
for i=SRcount+1:BondsCount
    h(j) = text(YTM.TimeFraction(i), minYTM(i).*100 - 0.15, YTM.Emission{i});
    j = j + 1;
end
set(h, 'FontSize', FontSize, 'FontName', 'Arial');
xlim([0 max(Mat)]);
ymin = min([ZeroRates, ParRates, InstFwdRates]);
ymax = max([ZeroRates, ParRates, InstFwdRates]);
ylim([ymin*100-0.5 ymax*100+0.5]);
% ======================================================
xlabel('Maturity in years')
ylabel('Percents')
grid off
hold off
if SRcount ~= 0
    legend('Spot Rates', 'Par Rates', 'Instantaneous Forward Rates',...
           'Observed Yields to Maturity', 'Fitted Yields to Maturity', 'Short-Rates')
else
    legend('Spot Rates', 'Par Rates', 'Instantaneous Forward Rates',...
           'Observed Yields to Maturity', 'Fitted Yields to Maturity')

end
legend('Location', 'NorthWest')
legend('boxoff')

end

author: Sable

repo: mcbench-benchmarks


text:
function h = fun2hist(fun,S)
%FUN2HIST Generates a histogram from a given digital function.
%   H = FUN2HIST(FUN,S) generates histogram H from a 1-D input
%   function FUN (a vector). The number of bins in H is equal to the
%   number of elements of FUN. If only FUN is provided in the input,
%   then H is normalized so that the sum of its components equals 1.
%   If S (a scalar) is provided, H is unnormalized, in the sense
%   that the sum of its components is equal to S. This is useful,
%   for example, if it is required that the sum of the components of
%   the histogram equal the number of pixels in an image. In this
%   case, the value of S would be S=M*N, where M and N are the
%   number of rows and columns in the image, respectively.
%
%   When S is specified, the elements of H are converted to
%   integers, so there is likely to be small roundoff errors between
%   the shapes of FUN and H. The reason for converting to integers
%   is that the sum of the components of an unnormalized histogram
%   has to equal the number of pixels (an integer) in the
%   corresponding image, as noted above.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Preliminaries
% Set defaults.
mode = 'u';
if nargin == 1
   mode = 'n';
end
% Number of bins (intensity levels in the histogram).
L = numel(fun);
% Initialize histogram.
h(1:L) = 0;
% Get index of the elements of fun that are not 0. These will be the
% populated bins in H.
idx = find(fun ~= 0);
IDXL = numel(idx);

% Generate the histogram from input function fun.
switch mode
   case 'n'
      h = fun/sum(fun);
   case 'u'
      % Convert to integers.
      h = round(S*(fun/sum(fun)));
      % Most likely, sum(H) will not equal S because of rounding.
      % If sum(H) is less than S, distribute the difference
      % between S and sum(H) equally (i.e., by increments of 1)
      % among elements of the populated histogram bins, starting
      % from the left and proceesing to the right. If sum(H) > S,
      % then take away 1 pixel from the populated bins.
      D = sum(h) - S;
      if D < 0
         count = abs(D);
         % Loop through the histogram, adding elements to the
         % populated bins until all D elements have been added.
         while count
            K = IDXL;
            if count < IDXL
               K = count;
            end
            for I = 1:K
               % Add counts to the populated bins only.
               h(idx(I)) = h(idx(I)) + 1;
               count = count - 1;
            end
         end
      elseif D > 0
         count = D;
         % Loop through the histogram, subtracting elements until
         % all D elements have been subtracted.
         while count
            K = IDXL;
            if count < IDXL
               K = count;
            end
            for I = 1:K
               % Subtract counts from the populated bins only.
               % But make sure they don't go negative.
               h(idx(I)) = h(idx(I)) - 1;
               if h(idx(I)) < 0
                  % Restore the amount subtracted.
                  h(idx(I)) = h(idx(I)) - 1;
                  % And reduce count so that the count will be
                  % subtracted elsewhere.
                  count = count + 1; % As if nothing had happened.
               end
            count = count - 1;
            end
         end
      end
end

author: dipum

repo: dipum-toolbox

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


text:
function sol = Inverse_Kinematics(q,pd,l)
%
px = pd(1);
py = pd(2);
l1 = l(1);
l2 = l(2);
distance = px^2 + py^2;
if distance > (l1+l2)^2
    k = pd/norm(pd);
    p = (l1+l2)*k;
    px = p(1);
    py = p(2);
elseif distance < (l1-l2)^2
    k = pd/norm(pd);
    p = (l1-l2)*k;
    px = p(1);
    py = p(2);
end

cosq2 = (distance - l1^2 - l2^2)/2/l1/l2;
sinq2 = sqrt(1-cosq2^2);
sincosq2 = zeros(2,2);
sincosq2(:,1) = [sinq2;cosq2];
sincosq2(:,2) = [-sinq2;cosq2];
solutions = zeros(2,4);
for i=1:2
    sinq2 = sincosq2(1,i);
    cosq2 = sincosq2(2,i);
    A = l1 + l2*cosq2;
    B = l2*sinq2;
    sinq1 = (px*A - py*B)/(A^2+B^2);
    cosq1 = (px*B + py*A)/(A^2+B^2);
    q1 = atan2(sinq1,cosq1);
    q2 = atan2(sinq2,cosq2);
    solutions(:,i) = [q1;q2];
end

if solutions(1,1) > solutions(1,2)
    solutions(1,3) = solutions(1,1) - 2*pi;
    solutions(1,4) = solutions(1,2) + 2*pi;
else
    solutions(1,3) = solutions(1,1) + 2*pi;
    solutions(1,4) = solutions(1,2) - 2*pi;
end

if solutions(2,1) > solutions(2,2)
    solutions(2,3) = solutions(2,1) - 2*pi;
    solutions(2,4) = solutions(2,2) + 2*pi;
else
    solutions(2,3) = solutions(2,1) + 2*pi;
    solutions(2,4) = solutions(2,2) - 2*pi;
end

delta = 10000;
index = 0;
for i=1:4
    if norm(q - solutions(:,i))<delta
        delta = norm(q - solutions(:,i));
        index = i;
    end
end
sol = solutions(:,index);
end

author: xuhuairuogu

repo: V-REP-Simulation-Projects


text:
function P = subsasgn(P, index, val)
%
%  SUBSASGN  Define index assignment for psf object.
%
%               P = subsref(P, index, val)
%
%          This is called whenever an assignment statement of a
%          psf object is made, such as:
%              P(i) = val, i = 1, 2
%              P.fieldname = val, fieldname = image, center
%

%  J. Nagy 5/1/01

switch index.type
case '()'
  switch index.subs{:}
  case 1
    P.image = val;
  case 2
    P.center = val;
  otherwise
    error('Index out of range.')
  end

case '.'
  switch index.subs
  case 'image'
    P.image = val;
  case 'center'
    P.center = val;
  otherwise
    error('Invalid field names.');
  end

case '{}'
  error('Cell array indexing not supported for psf object.')
end

repo: IRtools


text:
function Population = EnvironmentalSelection(Population,N,objective_set)
% The environmental selection of SIBEA-kEMOSS

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

% This function is written by Liangli Zhen

    PopObj = Population.objs;

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

    %% Select the solutions in the last front based on their HV loss
    Last = find(FrontNo==MaxFNo);
    % Calculate the WHV loss of each solution in the last front
    HVLoss   = CalHVLoss(PopObj(Last,objective_set),FrontNo(Last));
    [~,Rank] = sort(HVLoss,'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;

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

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%kcmindist_classify 'Classify an object using the Minimum Distance algorithm'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros cmindist_classify.pane file
%
% Parameters:
% InputFile: iimage 'Input Image', required: 'input image file'
% InputFile: isigs 'Input Signatures', required: 'input signatures file'
% Integer: distancerank 'Use distance rank', default: 1: 'distance rank integer'
% OutputFile: oclass 'Output Classified', required: 'output classified image'
% OutputFile: oprob 'Output Probabilities', optional: 'output a posteriori probabilities'
% OutputFile: oinfo 'Output Information', optional: 'output information about classification'
%
% Example: [oclass, oprob, oinfo] = kcmindist_classify({iimage, isigs}, {'iimage','';'isigs','';'distancerank',1;'oclass','';'oprob','';'oinfo',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% cmindist_classify - Classify an object using the Minimum Distance algorithm
%
%  DESCRIPTION
% This routine classifies an object using the Minimum Distance classification algorithm (or, more precisely, the Minimum Euclidean Distance to Class Means algorithm). The Minimum Distance algorithm is a simple algorithm that assigns a class C to a pixel X if the distance of a prototype of C to the vector X is the smallest between all known classes. More details about the Minimum Distance classifier are on the Classify Toolbox Manual.
% This routine requires an input object to be classified (specified by the parameter [-iimage]) and a set of classes's signatures object created by "cmindist_signature" and appended by \fIkappend\fP (specified by the parameter [-isigs]).
% This routine will create the output classification result in the file specified by [-oclass]. Optionally the final distances for each point and each class can be created if the parameter [-oprob] is used. If a file is specified in the  parameter [-oinfo] the classification results will be written to that file in ASCII (can get large for large images).
% The expected dimensions of the input and output objects are shown below:
% The input object which will be classified will have dimensions WxHxDxTxF, where F is the number of features. If it has mask, the masked points won't be classified, and the corresponding points in the output will have value 0 and mask 0.
% The input signatures must have dimensions Fx2x1xNx1, where N is the number of classes. It must be created by using "cmindist_signature" and \fIkappend\fP.
% The output object (specified with the parameter [-oclass]) will have dimensions WxHxDxTx1, and its value segment will have values on the range 0..N, where 0 means that the pixel was rejected (see below). It will also have a corresponding mask segment with values 0 for the rejected pixels and 1 for the non-rejected pixels. The values for the pixels will be associated accordingly to the order the signatures were appended with "kappend".
% If the parameter [-oprob] is used, it will have dimensions WxHxDxTxN.
% Alternatively to the Minimum Distance classification, the second (or third, or N-th) distance result can be obtained by specifying an index in the [-distancerank] parameter.
% Points can be optionally rejected (with value and mask 0 in the output) if their distance to the center (or mean) in the signature is larger than an absolute value (specified with the [-reject] parameter) or a number of standard deviations (specified with the [-reject] parameter and the [-usestddev] flag).
% Please refer to the Classify toolbox manual or for one of the example workspaces for usage examples and details.
%
%
%
%  EXAMPLES
% All examples for the Classify toolbox are listed on the Classify Toolbox Manual. For examples of this program, please see the Classify:workspaces:MINDIST and Classify:workspaces:MINDIST-Classify example workspaces.
%
%  "SEE ALSO"
% cmindist_signature, kappend.
%
%  RESTRICTIONS
% Expects the signatures to be valid signatures for the Minimum Distance algorithm.
%
%  REFERENCES
% All references for the Classify toolbox are listed on the Classify Toolbox Manual.
%
%  COPYRIGHT
% Copyright (C) 1997 Rafael Santos. Khoros (C) Khoral Research, Inc.
%

function varargout = kcmindist_classify(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kcmindist_classify(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'iimage', '__input';'isigs', '__input';'distancerank', 1;'oclass', '__output';'oprob', '__output';'oinfo', '__output'};
maxval={0,0,2,0,1,1};
minval={0,0,2,0,1,1};
istoggle=[0,0,1,0,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','Integer','OutputFile','OutputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; 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(1,1);
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
[varargout{:}]=callKhoros([w 'cmindist_classify"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function p00_tan_test ( problem_num )

%*****************************************************************************80
%
%% P00_TAN_TEST computes and tests the tangent vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM_NUM, the number of problems.
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'P00_TAN_TEST\n' );
  fprintf ( 1, '  Compute the tangent vector TAN(X) at the starting point.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Verify that JAC(X) * TAN(X) = 0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Verify that det ( JAC ) > 0\n' );
  fprintf ( 1, '                  ( TAN )    \n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   Problem    Option    ||Jac*Tan||     det(Jac|Tan)\n' );

  for problem = 1 : problem_num

    option_num = p00_option_num ( problem );

    fprintf ( 1, '\n' );

    for option = 1 : option_num

      nvar = p00_nvar ( problem, option );

      x0 = p00_start ( problem, option, nvar );

      jac = p00_jac ( problem, option, nvar, x0 );

      tan = p00_tan ( problem, option, nvar, x0 );

      jt = jac * tan;
      jtn = norm ( jt );

      jac(nvar,1:nvar) = tan(1:nvar);
      jtd = det ( jac );

      fprintf ( 1, '  %8d  %8d  %14e  %14e\n', problem, option, jtn, jtd );

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function IEEE80211n_open_graphics(blk, forceClose);
% forceClose: 1 = force figure to close; 0 = toggle open/close
gfig = findobj('type', 'figure', 'tag', blk);
if isempty(gfig)
    if ~forceClose
        gfig = openfig(get_param(blk, 'graphicsName'));
        set(gfig, 'tag', blk);
    end
else
    close(gfig);
end


text:
function [ c, seed ] = i4mat_uniform_ab ( m, n, a, b, seed )

%*****************************************************************************80
%
%% I4MAT_UNIFORM_AB returns a scaled pseudorandom I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley Interscience, page 95, 1998.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer M, N, the row and column dimensions of the matrix.
%
%    Input, integer A, B, the minimum and maximum acceptable values.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, integer C(M,N), the randomly chosen integer vector.
%
%    Output, integer SEED, the updated seed.
%
  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_UNIFORM_AB - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'I4MAT_UNIFORM_AB - Fatal error!' );
  end

  seed = floor ( seed );
  a = round ( a );
  b = round ( b );

  for j = 1 : n

    for i = 1 : m

      seed = mod ( seed, 2147483647 );

      if ( seed < 0 )
        seed = seed + 2147483647;
      end

      k = floor ( seed / 127773 );

      seed = 16807 * ( seed - k * 127773 ) - k * 2836;

      if ( seed < 0 )
        seed = seed + 2147483647;
      end

      r = seed * 4.656612875E-10;
%
%  Scale R to lie between A-0.5 and B+0.5.
%
      r = ( 1.0 - r ) * ( min ( a, b ) - 0.5 ) ...
        +         r   * ( max ( a, b ) + 0.5 );
%
%  Use rounding to convert R to an integer between A and B.
%
      value = round ( r );

      value = max ( value, min ( a, b ) );
      value = min ( value, max ( a, b ) );

      c(i,j) = value;

    end
  end

  return
end


text:
function splitNSx(splitCount)

% splitNSx
%
% Opens and splits an NSx file in smaller pieces, timewise.
%
% Use splitNSx(splitCount)
%
% All input arguments are optional. Input arguments can be in any order.
%
%   splitCount:   Defines the number of splits.
%                 DEFAULT: Splits the file in 2 pieces.
%
%   Example 1:
%   splitNSx(4);
%
%   In the example above, the user will be prompted to select a file. The
%   loaded file will be split in 4 samller files. For example, if the file
%   is 1 hour long then it will be split into four 15-minute files.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Kian Torab
%   support@blackrockmicro.com
%   Blackrock Microsystems
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Version History
%
% 1.0.0.0:
%   - Initial release.
%
% 1.1.0.0:
%   - Fixed a bug related to a case where initial timestamp of the first
%     data segment was not 0.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

% Validating input parameter
if ~exist('splitCount', 'var')
    splitCount = 2;
end

% Getting the file name
[fname, path] = getFile('*.nev', 'Choose an NSx file...');
if fname == 0
    disp('No file was selected.');
    if nargout
        clear variables;
    end
    return;
end
fext = fname(end-3:end);

% Loading the file
%% Reading Basic Header from file into NSx structure.
FID                       = fopen([path fname], 'r', 'ieee-le');
BasicHeader               = fread(FID, 336, '*uint8');
Trackers.fExtendedHeader  = double(typecast(BasicHeader(13:16), 'uint32'));
Trackers.countPacketBytes = double(typecast(BasicHeader(17:20), 'uint32'));

%% Doing Trackers
fseek(FID, 0, 'eof');
Trackers.fData = ftell(FID);
Trackers.countDataPacket = (Trackers.fData - Trackers.fExtendedHeader)/Trackers.countPacketBytes;

fseek(FID, Trackers.fExtendedHeader, 'bof');
tRawData  = fread(FID, [10 Trackers.countDataPacket], '10*uint8=>uint8', Trackers.countPacketBytes - 10);
Timestamp = tRawData(1:4,:);
Timestamp = typecast(Timestamp(:), 'uint32').';

splitPacketStarts = find(diff(Timestamp)<0);
splitPacketBytes = Trackers.countPacketBytes * splitPacketStarts;

% Reading headers and seeking to beginning of data
fseek(FID, 0, 'bof');
fileHeader = fread(FID, Trackers.fExtendedHeader, '*uint8');

for idx = 1:length(splitPacketStarts)
    % Opening a file for saving
    FIDw = fopen([path fname(1:end-4) '-s' sprintf('%03d', idx) fname(end-3:end)], 'w+', 'ieee-le');
    fprintf('\nReading segment %d... ', idx);
    % Reading the segment
    dataSegment = fread(FID, splitPacketBytes(idx), 'char');
    fprintf('Writing segment %d... ', idx);
    % Writing the segmented data into file
    fwrite(FIDw, fileHeader, 'char');
    fwrite(FIDw, dataSegment, 'char');
    % Clearing variables and closing file
    clear dataSegment;
    fclose(FIDw);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function egrad = getEuclideanGradient(problem, x, storedb, key)
% Computes the Euclidean gradient of the cost function at x.
%
% function egrad = getEuclideanGradient(problem, x)
% function egrad = getEuclideanGradient(problem, x, storedb)
% function egrad = getEuclideanGradient(problem, x, storedb, key)
%
% Returns the Euclidean gradient at x of the cost function described in the
% problem structure.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% Because computing the Hessian based on the Euclidean Hessian will require
% the Euclidean gradient every time, to avoid overly redundant
% computations, if the egrad function does not use the store caching
% capabilites, this implements an automatic caching functionality. Writing
% egrad to accept the optional store or storedb parameter will disable
% automatic caching, but allow user controlled caching.
%
% See also: getGradient canGetGradient canGetEuclideanGradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, July 9, 2013.
% Contributors:
% Change log:
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.
%
%   June 28, 2016 (NB):
%       Added support for getPartialEuclideanGradient
%
%   July 26, 2018 (NB):
%       The Euclidean gradient is now automatically cached if the Euclidean
%       Hessian is also computable from the problem description. This
%       differs from previous behavior where it would only be cached if
%       problem.egrad did not accept store or storedb as input; the
%       converse was taken as a sign that the user wants to deal with
%       caching on their own, but in reality it proved more confusing than
%       helpful.
%
%   Sep.  6, 2018 (NB):
%       The Euclidean gradient is now always cached. Caching conditioned on
%       the Euclidean Hessian being provided was problematic, because if no
%       Hessian is provided at all, then the default fall-back of finite
%       differences would be slowed down significantly without cache. Since
%       the new storedb functionalities around storedb.remove() much reduce
%       the number of stored stores, this should not be an issue.

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    % Contrary to most similar functions, here, we get the store by
    % default. This is for the special caching functionality described
    % below.
    store = storedb.getWithShared(key);
    store_is_stale = false;

    % We force caching of the Euclidean gradient. Look up first.
    force_egrad_caching = true;
    if force_egrad_caching && isfield(store, 'egrad__')
        egrad = store.egrad__;
        return;
    end

    if isfield(problem, 'egrad')
    %% Compute the Euclidean gradient using egrad.

        % Check whether this function wants to deal with storedb or not.
        switch nargin(problem.egrad)
            case 1
                egrad = problem.egrad(x);
            case 2
                [egrad, store] = problem.egrad(x, store);
            case 3
                egrad = problem.egrad(x, storedb, key);
                % The store structure in storedb might have been modified
                % (since it is passed by reference), so before caching
                % we'll have to update (see below).
                store_is_stale = true;
            otherwise
                up = MException('manopt:getEuclideanGradient:badegrad', ...
                    'egrad should accept 1, 2 or 3 inputs.');
                throw(up);
        end

    elseif canGetPartialEuclideanGradient(problem)
    %% Compute the Euclidean gradient using a full partial Euclidean gradient.

        d = problem.ncostterms;
        egrad = getPartialEuclideanGradient(problem, x, 1:d, storedb, key);
        store_is_stale = true;

    else
    %% Abandon computing the Euclidean gradient

        up = MException('manopt:getEuclideanGradient:fail', ...
            ['The problem description is not explicit enough to ' ...
             'compute the Euclidean gradient of the cost.']);
        throw(up);

    end

    % If we are not sure that the store structure is up to date, update.
    if store_is_stale
        store = storedb.getWithShared(key);
    end

    % Cache here.
    if force_egrad_caching
        store.egrad__ = egrad;
    end

    storedb.setWithShared(store, key);

end

author: NicolasBoumal

repo: manopt


text:
function varargout=corr(this,varargin)
% Computes correlation for time series: It is an interface to MATLAB implementation of corr
%
% ::
%
%    varargout = corr(db,varargin);
%
% Args:
%
%    db (ts object): times series object
%
%    varargin: varargin for corr function in MATLAB
%
% Returns:
%    :
%
%    varargout: output from corr function
%

if ~isempty(varargin) && isa(varargin{1},'ts')

    this=this & varargin{1};

    varargin=varargin(2:end);

end

data=this.data;

if size(data,3)>1

    error([mfilename,':: this operation is only defined for databases with one page'])

end

nout=nargout;

bad=any(isnan(data),2);

if any(bad)

    data(bad,:)=[];

    disp([int2str(sum(bad)),' nan rows above discarded'])

end

[varargout{1:nout}]=utils.stat.corr(data,varargin{:});

end

author: jmaih


text:
function [x, infos] = incremental_mu_nmf(V, rank, in_options)
% Incremental non-negative matrix factorization (NMF) with outliers (INCNMF) algorithm.
%
% Inputs:
%       matrix      V
%       rank        rank
%       in_options  options
% Output:
%       w           solution of w
%       infos       information
%
% References:
%       S. S. Bucak, B. Gunsel,
%       "Incremental Subspace Learning via Non-negative Matrix Factorization,"
%       Pattern Recognition, 2009.
%
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on Feb. 12, 2017
%
% Change log:
%
%       Oct. 27, 2017 (Hiroyuki Kasai): Fixed algorithm.
%
%       May. 20, 2019 (Hiroyuki Kasai): Added initialization module.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.max_inneriter     = 1;
    local_options.online            = 0;
    local_options.tolcostdegrease   = 1e-8;
    local_options.alpha             = 0.5;
    local_options.beta              = 0.5;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    W = init_factors.W;
    H = init_factors.H;

    % initialize
    method_name = 'Incremental-MU';
    epoch = 0;
    grad_calc_count = 0;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    if options.online
        ht = H(:, end); %why
    else
        % Do nothing
    end
    A = V * H';
    B = H * H';

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, W, H, [], options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s: Epoch = 0000, cost = %.16e, optgap = %.4e\n', method_name, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main outer loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % main inner loop
        for t = 1 : options.batch_size : n - 1

            vt = V(:, t:t+options.batch_size -1);

            % Need to be considered more carefully
            if options.online
                % Do nothing (??)
            else
                ht = H(:, t:t+options.batch_size-1);
            end

            for j = 1 : options.max_inneriter

                % Compute new ht (column vector at t)
                ht = ht .* (W' * vt) ./ (W' * (W * ht) + 1e-9);
                ht = ht + (ht<eps) .* eps;

                % Compute new W
                vht = vt * ht';
                hht = ht * ht';
                A_tmp = options.beta * A + options.alpha * vht;
                B_tmp = options.beta * B + options.alpha * hht;
                W = W .* ( A_tmp ./ (W * B_tmp + 1e-9) ); % Add 1e-9 to avoid 0 in the denom.
                W = W + (W<eps) .* eps;

               if j > 1
                    oldobj = newobj;
                end
                newobj = ((sum(sum((vt-W*ht).^2)))/m);

                if options.verbose > 2
                    fprintf('\t[%d-%d-%d] %e\n', epoch, t, j, newobj);
                end

                if j > 1 && (oldobj-newobj)/newobj < options.tolcostdegrease
                    break;
                end

            end

            grad_calc_count = grad_calc_count + m * options.batch_size;

            % update A and B
            A = options.beta * A + options.alpha * vht;
            B = options.beta * B + options.alpha * hht;

            % update ht (columb of H at t)
            H(:,t:t+options.batch_size-1) = ht;

            if options.verbose > 2
                % measure cost
                f_val = nmf_cost(V, W, H, []);
                fprintf('%s: inner [%03d-%03d] %e\n', method_name, epoch, t, f_val);
            end
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, W, H, [], options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = W;
    x.H = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function t = Xaxis_string(Xt,bnames,R,fs);

ax = axis;
set(gca,'XLimMode','manual')
Yl = ax(3:4);
t = text(Xt,Yl(1)*ones(1,length(Xt)),bnames);
if R == 0
  set(t,'HorizontalAlignment','center','VerticalAlignment','top', 'Rotation',R,'fontsize',fs);
elseif R == 90
  set(t,'HorizontalAlignment','right','VerticalAlignment','middle', 'Rotation',R,'fontsize',fs);
else
  set(t,'HorizontalAlignment','right','VerticalAlignment','top', 'Rotation',R,'fontsize',fs);
end
set(gca,'XTickLabel','')

author: andrewssobral


text:
function [h1, h2] = ivmPrintPlot(model, plotType, iter, X, ...
                                 y, capName, experimentNo)

% IVMPRINTPLOT Make a 3-D or contour plot of the IVM.
% FORMAT
% DESC makes a 3-D or contour plot from an IVM model to show
% results and prints it out to various directories (in eps and png form).
% ARG model : the model from which the plot is generated.
% ARG plotType : type of plot to be used, options include
% 'ncnmContour', which gives the contours at the edge of the null
% category region for the null category noise model, and
% 'ivmContour'.
% ARG iter : iteration number, if give it is used at the title of
% the plot.
% ARG X : optional argument, if given it is the plotted input
% locations, otherwise model.X is used.
% ARG y : optional argument, if given it is the plotted target
% locations, otherwise model.y is used.
% ARG capName : the name of the saved plots.
% ARG experimentNo : the experiment number to assign to the files.
% RETURN h1 : a handle to the data in the plot.
% RETURN h2 : a handle to any contours, or functions plotted
% derived from the IVM.
%
% SEEALSO : noise3dPlot, noisePointPlot, ivmMeshVals, ivmCreate
%
% COPYRIGHT : Neil D. Lawrence, 2007

% IVM

ivm3dPlot(model, plotType, iter, X, y);
% display active points.
model = ivmOptimiseIvm(model, 2);
fileName = ['dem' capName 'Ivm' num2str(experimentNo)];

printPlot(fileName, '../tex/diagrams', '../html');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%% SATISFY_SCRIPT is a script to call SATISFY_FUN.
%
%  Discussion:
%
%    The BATCH command runs scripts, not functions.  So we have to write
%    this short script if we want to work with BATCH!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SATISFY_SCRIPT:\n' );
  fprintf ( 1, '  Call SATISFY_FUN for 0/1 variable assignments that satisfy a formula.\n' );

  tic
  solutions = satisfy_fun ( );
  toc


text:
function figbolden(pos,pv,m)
% embolden the current figure
% Inputs: pos = [xmin ymin width height] gives the lower left corner position and the window size in pixels
%               [width height] leaves the lower left corner alone
%               [width] has a standard aspect ratio of 4:3
%               [-width/height] leaves the area unchanged but fixes the aspect ratio
%         pv is a cell name containing attribute-value pairs.
%            default = {'FontName' 'Ariel'; 'FontSize' 16; 'LineWidth' 2; 'MarkerSize' 8}
%         m is a mode string:
%                'l' lists the changes made
%                'd' uses default pv settings
%                'x' suppresses changes

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

ps={'Title' 'XLabel' 'YLabel' 'Children'};
if nargin<3
    m='';
end
if nargin<2 || any(m=='d')
    pv={'FontName' 'Arial'; 'FontSize' 16; 'LineWidth' 2; 'MarkerSize' 8};
end
pp={'Symbol';'Wingdings'};      % protected fonts
if nargin<1
    pos=[];
end
mlist=any(m=='l');  % list changes
mnotx=~any(m=='x'); % do changes
scsz=get(0,'screensize');
if length(pos)
    po=get(gcf,'position');
    if length(pos)>2            % position is specified
        po(1:2)=pos(1:2);
        pos(1:2)=[];      % remove xmin,ymin
    end
    if length(pos)>1
        po(3:4)=pos(1:2);
    else
        if pos(1)>0
            po(3:4)=[1 0.75]*pos(1);
        else
            po(3:4)=[-pos(1) 1]*sqrt(-po(3)*po(4)/pos(1)); % preserve area
        end
    end
    set(gcf,'position',po);
end
hlist=get(gcf,'children');
while length(hlist)
    pl=get(hlist(1));
    %fprintf('list length = %d, handle = %f\n',length(hlist),hlist(1));
    for i=1:size(pv,1)
        if isfield(pl,pv{i,1})
            if i>1 || all(~strcmpi(get(hlist(1),pv{i,1}),pp))
                pval=get(hlist(1),pv{i,1});
                if ~all(size(pval)==size(pv{i,2})) || ~all(pval(:) == pv{i,2}(:))
                    if mnotx
                        set(hlist(1),pv{i,1},pv{i,2})
                    end
                    if mlist
                        if ischar(pval)
                            fch='s';
                        else
                            fch='g';
                        end
                        fprintf(['change %f %s: %' fch ' -> %' fch '\n'],hlist(1),pv{i,1},pval,pv{i,2});
                    end
                end
            end
        end
    end
    for i=1:length(ps)
        if isfield(pl,ps{i})
            hlist=[hlist; get(hlist(1),ps{i})];
            %fprintf('add %f:%s\n',hlist(1),ps{i});
        end
    end
    hlist(1)=[];
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
%IM_MEASURE Fixed mapping computating features by DIP_Image
%
%		F = IM_MEASURE(A,GRAY,FEATURES)
%
% INPUT
%   A        Dataset with binary object images dataset (possibly multi-band)
%   GRAY     Gray-valued images (matched with A, optional)
%   FEATURES Features to be computed
%
% OUTPUT
%   F        Dataset with computed features
%
% In each image of the measurement set GRAY the features given in FEATURES
% are measured. In A a segmented version of GRAY has to be supplied.
% When no GRAY is supplied, the binary images in A are used. Only
% the largest object in each image is considered.
%
% The following features may be computed:
%    'dimension','mean','stddev','gravity','size','center','max','min',
%    'maxval','minval','feret'','inertia','ccbendingenergy'.
% Note that some features like 'mean' (mean image intensity) and 'stddev'
% (standard deviation of image intensity) are not useful for binary images.
% Run MEASUREHELP to get some information on these measures.
%
% Use FEATURES = 'all' for computing all features.
% Use MEASUREHELP for some description of the features.
% Use IM_FEATURES for a set of features not based on DIP_Image
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES, MEASURE, MEASUREHELP

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

function b = im_measure(a,gray,features)

	checktoolbox('diplib');
	if nargin < 3  features = []; end
	if nargin < 2  gray = []; end
	%if nargin < 2 | isempty(gray), gray = a; end

	if nargin < 1 | isempty(a)
		b = prmapping(mfilename,'fixed',{gray,features});
		b = setname(b,'DIP measurements');
	elseif isdataset(a)
		if (nargin < 3 | isempty(features)) & ~isdataset(gray)
			features = gray;
			gray = a;
		end
		if ~isdataset(gray)
			error('Binary and gray images should be both datasets')
		end
		fsize = getfeatsize(a);
		if any(getfeatsize(gray) ~= fsize)
			error('Image structures of binary and gray images should be identical')
		end
		if length(fsize) == 2, fsize = [fsize 1]; end
		if size(a,1) ~= size(gray,1)
			error('Same number of binary and gray images expected')
		end
		out = [];
		binim = data2im(a);
		grim = data2im(gray);
		nim = size(a,1)*fsize(3);
		s = sprintf('Measuring %i images',nim);
		prwaitbar(nim,s);
		for i=1:size(a,1)
			for j=1:fsize(3)
				prwaitbar(nim,(i-1)*fsize(3)+j);
				f = feval(mfilename,binim(:,:,j,i),grim(:,:,j,i),features);
				if isempty(out)
					out = repmat(f(:)',[size(a,1),1,fsize(3)]);
					%out = reshape(f(:)',[size(a,1),1,fsize(3)]);
				else
					out(i,:,j) = f;
				end
			end
		end
		prwaitbar(0);
		b = setdat(a,out);
		b = setfeatsize(b,[length(f),fsize(3)]);
		b = setfeatlab(b,getfeaturelabels(features));
	elseif isdatafile(a)
		if nargin < 3 | isempty(features) & ~isdatafile(gray)
			features = gray;
			gray = a;
		end
		if ~isdatafile(gray)
			error('Binary and gray images should be both datafiles')
		end
		b = dyadic(a,mfilename,gray,{features});
		%b = setfeatlab(b,getfeaturelabels(features));
  elseif isa(a,'double') | isa(a,'dip_image') % here we have a single image
		if isempty(features), features = 'dimension'; end
		gray = 1.0*dip_image(gray);
		%labim = label(dip_image(im_select_blob(a),'bin'));
		labim = label(dip_image(a,'bin'));
    c = measure(labim,gray,'size',[],2);
    labid = c.id;
    sz = c.size;
    [bb,mm] = max(sz);
    labid = labid(mm);
		if strcmp(features,'all')
			features = {'dimension','mean','stddev','gravity',...
			'size','center','max','min', 'maxval','minval',...
			'feret','inertia', 'ccbendingenergy'};
		end
		b = measure(labim,gray,features,labid,2);
    b = double(b);
	else
		error('Wrong input')
	end

return

function names = getfeaturelabels(features)

names = {};
for i=1:length(features)
	switch features{i}
	case 'dimension'
		names{end+1} = 'imagewidth';
		names{end+1} = 'imageheight';
	case 'mean'
		names{end+1} = 'mean int';
	case {'mean', 'sum'}
		names{end+1} = 'mass';
	case 'stddev'
		names{end+1} = 'standard dev.';
	case 'gravity'
		names{end+1} = 'gravity x';
		names{end+1} = 'gravity y';
	case 'size'
		names{end+1} = 'size';
	case 'center'
		names{end+1} = 'center x';
		names{end+1} = 'center y';
	case 'max'
		names{end+1} = 'max x coord';
		names{end+1} = 'max y coord';
	case 'min'
		names{end+1} = 'min x coord';
		names{end+1} = 'min y coord';
	case 'maxval'
		names{end+1} = 'max int';
	case 'minval'
		names{end+1} = 'min int';
	case 'perimeter'
		names{end+1} = 'perimeter';
	case 'feret'
		names{end+1} = 'max diameter';
		names{end+1} = 'min diameter';
		names{end+1} = 'max perp. diameter';
	case 'inertia'
		names{end+1} = 'inertia moment 1';
		names{end+1} = 'inertia moment 2';
	case 'ccbendingenergy'
		names{end+1} = 'bending energy perimeter';
	otherwise
		error('I do not know feature %s.',features{i});
	end
end

return

author: marianux

repo: ecg-kit

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


text:
function reconOffline(pfileDir,writeEfileFlag,bigEndian,seqType);
% reconOffline([pfileDir,writeEfileFlag,bigEndian,seqType]);
%
% (For linux and Mac OSX only):
% Use gary's recon code (in /biac2/kgs/dataTools/) to
% make Pmags and E-file headers from Pfiles without
% logging into lucas.
%
% pfileDir: where to find the Pfiles. Defaults to the
% current directory, and also checks if there's a Raw/Pfiles
% directory relative to that directory. So you can run this
% straight from the session directory.
%
% writeEfileFlag: if 0, won't re-write the E-file headers. By
% default, it does (1).
%
% bigEndian: flag on whether the Pfiles are big-endian or
% little-endian format. Short answer: if the data were collected
% (at stanford lucas center) prior to late February 2005, set
% this flag to 1; otherwise, ignore it (it defaults to 0).
%
% seqType: flag for sequence type, when the sequence has special needs.
% defaults to 0 (standard sequence), 1 = 3D sequence (needs to run 3dcode),
% other values currently ignored...
%
% This uses grecons.lnx and the -O (one file) and -b (byte swap, *Endian dependent)
% flags -- otherwise the Pmags come out messy.
%
% 09/04 ras.
% 03/05 ras: updated to deal w/ little-endian files, after
% the 3T at Lucas switched to the Excite system.
% 01/07 dar: included seqType flag, updated call to target
% "grecons.lnx" (not "grecons")
% 10/10 kw/al: updated to version 95 of grecon(grecon15) and writeihdr
% To do:
% +could use a good way to auto-detect if 3dcode should be run.

if ~isunix
    fprintf('Sorry, unix systems only (linux and Mac OS X).\n');
    return
end

if ~exist('pfileDir','var') | isempty(pfileDir)
    pfileDir = pwd;
end

if ~exist('writeEfileFlag','var') | isempty(writeEfileFlag)
    writeEfileFlag = 1;
end

if ~exist('bigEndian','var') | isempty(bigEndian)
    bigEndian = 0;
end

if ~exist('seqType','var') | isempty(seqType)
    seqType = 0;
end

pattern = fullfile(pfileDir,'P*.7');
w = dir(pattern);

% may not be pfiles in the current directory, but
% in Raw/Pfiles. Check this as a backup:
if isempty(w)
    pfileDir = fullfile(pfileDir,'Raw','Pfiles');
    pattern = fullfile(pfileDir,'P*.7');
    w = dir(pattern);
end

% if still no pfiles, you're outta luck
if isempty(w)
    error(sprintf('No Pfiles found in %s.',pfileDir));
end

callingDir = pwd;
cd(pfileDir);

bfilepattern = fullfile(pfileDir,'B0*'); %check for preexisting B0* files.  If present, will not delete later
b = dir(bfilepattern);

for i = 1:length(w)
    fname = w(i).name;

    if writeEfileFlag==1
        cmd = sprintf('/biac2/kgs/dataTools/writeihdr15 %s',fname);
        unix(cmd);
        fprintf('Wrote E-file header for %s.\n',fname);
    end

    if seqType == 1
        %this block MUST come before recon call...creates files to be used for
        %3d recon - Original 'P*****.7' is converted to '3d_P*****.7' which becomes
        %the active file
        cmd = sprintf('/biac2/kgs/dataTools/3dcode %s 3d_%s',fname,fname);
        unix(cmd);
        fprintf('3dcode has converted %s to 3d_%s - ready for recon\n',fname,fname);
        oldfname = fname;
        fname = sprintf('3d_%s',fname);
    end

    if bigEndian
        % big endian: for scans before the upgrade in Feb 2005
        cmd = sprintf('/biac2/kgs/dataTools/grecons.lnx -O -b %s',fname);
    else
        % little endian: after the upgrade, Pfiles now don't
        % need to be byte-swapped
        cmd = sprintf('/biac2/kgs/dataTools/grecons15.lnx -O %s',fname);
    end
    unix(cmd);

    if seqType == 1
        cmd = sprintf('mv %s.mag %s.mag', fname, oldfname);
        unix(cmd);
        fprintf('\n\n\t*****Reconned %s to %s.mag*****\n\n',fname,oldfname);
        fname = oldfname;
    else
        fprintf('\n\n\t*****Reconned %s*****\n\n',fname);
    end

    if isempty(b)
        % move the B0 estimate files into their own subdirectory,
        % and rename the B0.* files to I.*, so code like loadVolume
        % which looks for the pattern I* can find them. (ras, 02/06)
        % modified to create subdirectory for each scan (dar 01/07)
        try
            if ~isempty(dir('B0.*'))
                dirname = sprintf('B0_%s',fname(1:6));
                cmd = sprintf('mkdir %s',dirname);
                unix(cmd)
                nfiles = length(dir('B0.*'));
                for i = 1:nfiles
                    src = sprintf('B0.%03.0f', i);
                    tgt = sprintf('%s/I.%03.0f',dirname, i);
                    movefile(src, tgt);
                end
            end
        end
    end
end
if ~isempty(dir('B0_*'))
    mkdir B0
    cmd = sprintf('mv B0_*/ B0/')
    unix(cmd)
end
cd(callingDir);
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% PERTURBK - Perturbs the current solution to the solution valid for the
%            given kernel parameter.
%
% Syntax: [a,b,g,ind,X_mer,y_mer,Rs,Q] = perturbk(scale)
%
%      a: alpha coefficients
%      b: bias
%      g: partial derivatives of cost function w.r.t. alpha coefficients
%    ind: cell array containing indices of margin, error and reserve vectors
%         ind{1}: indices of margin vectors
%         ind{2}: indices of error vectors
%         ind{3}: indices of reserve vectors
%  X_mer: matrix of margin, error and reserve vectors stored columnwise
%  y_mer: column vector of class labels (-1/+1) for margin, error and reserve vectors
%     Rs: inverse of extended kernel matrix for margin vectors
%      Q: extended kernel matrix for all vectors
%  scale: kernel scale
%
% Version 3.22e -- Comments to diehl@alumni.cmu.edu
%

function [a,b,g,ind,X,y,Rs,Q] = perturbk(new_scale)

% flags for example state
MARGIN    = 1;
ERROR     = 2;
RESERVE   = 3;
UNLEARNED = 4;

% define global variables
global a;             % alpha coefficients
global b;             % bias
global C;             % regularization parameters
global deps;          % jitter factor in kernel matrix
global g;             % partial derivatives of cost function w.r.t. alpha coefficients
global ind;           % cell array containing indices of margin, error, reserve and unlearned vectors
global num_unlearned; % number of unlearned vectors initially
global perturbations; % number of perturbations
global Q;             % extended kernel matrix for all vectors
global Rs;            % inverse of extended kernel matrix for margin vectors
global scale;			 % kernel scale
global type;          % kernel type
global X;             % matrix of margin, error, reserve and unlearned vectors stored columnwise
global y;             % column vector of class labels (-1/+1) for margin, error, reserve and unlearned vectors

kernel_evals_begin = kevals;
num_examples = size(X,2);

% sum{k in U} Qik lambda k and sum{k in U} yk lambda k
SQl = zeros(num_examples,1);
Syl = 0;

% adjust kernel scale
scale = new_scale;

% recompute g for the margin, error and reserve vectors
inda = [ind{MARGIN} ind{ERROR} ind{RESERVE}];
[f,K] = svmeval(X(:,inda));
g(inda) = y(inda).*f + a(inda)*deps - 1;

% identify the unlearned vectors and compute their coefficient sensitivities
lambda = zeros(num_examples,1);

% find all error vectors with g >= 0
flag = (g(ind{ERROR}) >= 0);
i = find(flag);

if (length(i) > 0)

   % relabel error vectors with g >= 0 as unlearned
   ind{UNLEARNED} = [ind{UNLEARNED} ind{ERROR}(i)];

   % coefficient sensitivities
   lambda(ind{ERROR}(i)) = -a(ind{ERROR}(i));

   % update sums
   SQl(inda) = SQl(inda) + ((y(ind{ERROR}(i))*y(inda)').*K(length(ind{MARGIN})+i,:))'*lambda(ind{ERROR}(i));
   Syl = Syl + y(ind{ERROR}(i))'*lambda(ind{ERROR}(i));

   % keep remaining error vectors labeled as such
   ind{ERROR}(i) = [];

end;

% find all reserve vectors with g <= 0
flag = (g(ind{RESERVE}) <= 0);
i = find(flag);

if (length(i) > 0)

   % relabel reserve vectors with g <= 0 as unlabeled
   ind{UNLEARNED} = [ind{UNLEARNED} ind{RESERVE}(i)];

   % coefficient sensitivities
   lambda(ind{RESERVE}(i)) = C(ind{RESERVE}(i));

   % update sums
   SQl(inda) = SQl(inda) + ((y(inda)*y(ind{RESERVE}(i))').*kernel(X(:,inda),X(:,ind{RESERVE}(i)),type,scale))*lambda(ind{RESERVE}(i));
   Syl = Syl + y(ind{RESERVE}(i))'*lambda(ind{RESERVE}(i));

   % keep remaining reserve vectors labeled as such
   ind{RESERVE}(i) = [];

end;

% find all margin vectors with g > 0
flag = (g(ind{MARGIN}) > 0);
i = find(flag);

if (length(i) > 0)

   % coefficient sensitivities
   lambda(ind{MARGIN}(i)) = -a(ind{MARGIN}(i));

   % update sums
   SQl(inda) = SQl(inda) + ((y(ind{MARGIN}(i))*y(inda)').*K(i,:))'*lambda(ind{MARGIN}(i));
   Syl = Syl + y(ind{MARGIN}(i))'*lambda(ind{MARGIN}(i));

end;

% find all margin vectors with g <= 0
flag = ~flag;
i = find(flag);

if (length(i) > 0)

   % coefficient sensitivities
   lambda(ind{MARGIN}(i)) = C(ind{MARGIN}(i))-a(ind{MARGIN}(i));

   % update sums
   SQl(inda) = SQl(inda) + ((y(ind{MARGIN}(i))*y(inda)').*K(i,:))'*lambda(ind{MARGIN}(i));
   Syl = Syl + y(ind{MARGIN}(i))'*lambda(ind{MARGIN}(i));

end;

% relabel margin vectors as unlearned
ind{UNLEARNED} = [ind{UNLEARNED} ind{MARGIN}];
ind{MARGIN} = [];

% add jitter factor
SQl(ind{UNLEARNED}) = SQl(ind{UNLEARNED}) + deps*lambda(ind{UNLEARNED});

% number of unlearned vectors initially
num_unlearned = length(ind{UNLEARNED});
s = sprintf('Number of unlearned vectors: %d',num_unlearned);
disp(s);

% reset Q and Rs
Q = Q(1,:);
Rs = Inf;

p_s = 0;
num_MVs = length(ind{MARGIN});
num_learned = 0;
perturbations = 0;
while ((length(ind{UNLEARNED}) > 0) | ((p_s < 1) & (length(ind{UNLEARNED}) == 0)))

   perturbations = perturbations + 1;

   % compute beta and gamma
   if (num_MVs > 0)

      v = zeros(num_MVs+1,1);
      if (p_s < 1-eps)
         v(1) = -Syl - sum(y.*a)/(1-p_s);
      else
         v(1) = -Syl;
      end;
      v(2:num_MVs+1) = -SQl(ind{MARGIN});
      beta = Rs*v;
      gamma = zeros(size(Q,2),1);
      ind_temp = [ind{ERROR} ind{RESERVE} ind{UNLEARNED}];
      gamma(ind_temp) = Q(:,ind_temp)'*beta + SQl(ind_temp);

   else

      beta = 0;
      gamma = SQl;

   end;

   % minimum acceptable parameter change
   [min_dps,indss,cstatus,nstatus] = min_delta_p_s(p_s,gamma,beta,lambda);

   % update a, b, g and p_s
   if (length(ind{UNLEARNED}) > 0)
      a(ind{UNLEARNED}) = a(ind{UNLEARNED}) + lambda(ind{UNLEARNED})*min_dps;
   end;
   if (num_MVs > 0)
      a(ind{MARGIN}) = a(ind{MARGIN}) + beta(2:num_MVs+1)*min_dps;
   end;
   b = b + beta(1)*min_dps;
   g = g + gamma*min_dps;
   p_s = p_s + min_dps;

   % perform bookkeeping
   indco = bookkeeping(indss,cstatus,nstatus);

   % set g(ind{MARGIN}) to zero
   g(ind{MARGIN}) = 0;

   % update Rs and Q if necessary
   if (nstatus == MARGIN)

      num_MVs = num_MVs + 1;
      if (num_MVs > 1)

         % compute beta and gamma for indss
         beta = -Rs*Q(:,indss);
         gamma = kernel(X(:,indss),X(:,indss),type,scale) + deps + Q(:,indss)'*beta;

      end;

      % expand Rs and Q
      updateRQ(beta,gamma,indss);

   elseif (cstatus == MARGIN)

      % compress Rs and Q
      num_MVs = num_MVs - 1;
      updateRQ(indco);

   end;

   % update SQl and Syl when the current status of indss is UNLEARNED
   if (cstatus == UNLEARNED)
      num_learned = num_learned + 1;
      if (nstatus == MARGIN)
         SQl = SQl - Q(num_MVs+1,:)'*lambda(indss);
      else
         SQl = SQl - ((y*y(indss)).*kernel(X,X(:,indss),type,scale))*lambda(indss);
         SQl(indss) = SQl(indss) - deps*lambda(indss);
      end;
      Syl = Syl - y(indss)*lambda(indss);

      if (mod(num_learned,50) == 0)
         s = sprintf('Learned %d examples.',num_learned);
         disp(s);
      end;
   end;

end;
disp('Perturbation complete!');

% summary statistics
s = sprintf('\nMargin vectors:\t\t%d',length(ind{MARGIN}));
disp(s);
s = sprintf('Error vectors:\t\t%d',length(ind{ERROR}));
disp(s);
s = sprintf('Reserve vectors:\t%d',length(ind{RESERVE}));
disp(s);
s = sprintf('Kernel evaluations:\t%d\n',-kernel_evals_begin+kevals);
disp(s);


text:
function test_suite = test_mirdwt
  disp("mrdwt")
  test_mirdwt_1
  test_mirdwt_2

function test_mirdwt_1
       xin = makesig('Leopold',8);
       h = daubcqf(4,'min');
       Lin = 1;
       [yl,yh,L] = mrdwt(xin,h,Lin);
       [x,L] = mirdwt(yl,yh,h,L);

assertEqual(L,Lin);
assertVectorsAlmostEqual(x, xin,'relative',0.0001);

function test_mirdwt_2
       load ../lena512;
       x = lena512;
       h = daubcqf(6);
       [yl,yh,L] = mrdwt(x,h);
assertEqual(L,9);
       [x_new,L] = mirdwt(yl,yh,h);
assertEqual(L,9);
assertVectorsAlmostEqual(x, x_new,'relative',0.0001);

author: ricedsp


text:
function [] = mci_plot_dist_multi (dist,name,P)
% Plot (multiple) densities
% FORMAT [] = mci_plot_dist_multi (dist,name,P)
%
% dist{i}   ith distribution
%
% .Ep       mean
% .P        [Np x Ns] sample matrix
% .ind      indices of samples (eg. post burn-in)
% .names    names of variables
% .color    eg 'r','k','b'
% .order    eg. [1,3,4,2] to plot only variables 1,3,4 and 2 in that order
%
% name      name of parameters
% P         true parameters (optional)
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: mci_plot_dist_multi.m 6697 2016-01-27 14:57:28Z spm $

Ndist=length(dist);
for i=1:Ndist
    if strcmpi(dist{i}.type,'sample')
        figure;
        plot(dist{i}.P');
        xlabel('Sample');
        ylabel('Parameters');
        title(sprintf('%s trajectories',name));
    end
end

Np=length(dist{1}.Ep);
try
    order=dist{1}.order;
    Np=length(order);
catch
    order=[1:Np];
end

% Plot univariate densities
rNp=ceil(sqrt(Np));
h=figure;
set(h,'Name',name);
for j=1:Np,
    subplot(rNp,rNp,j);
    jp=order(j);
    for i=1:Ndist,
        mci_plot_dist(dist{i},jp);
        hold on
    end
    yl=get(gca,'YLim');
    if nargin > 2
        plot([P(jp),P(jp)],[0,yl(2)],'r','LineWidth',2);
    end
end
if nargin > 2
    disp('True parameters shown in red');
end

plot_bivariate=0;
if plot_bivariate
    % Plot bivariate densities
    rNp=ceil(sqrt(Np));
    h=figure;
    set(h,'Name',name);
    k=1;
    q=dist{1}.P(:,dist{1}.ind);
    for i=1:Np,
        for j=1:Np,
            if j > i
                subplot(Np,Np,k);
                plot(q(i,:),q(j,:),'k.');
                xlabel(dist{1}.names{i});
                ylabel(dist{1}.names{j});
                hold on
                plot(P(i),P(j),'ro');
            end
            k=k+1;
        end
    end
    if nargin > 2
        disp('True parameters shown in red');
    end

    % Posterior Correlation Matrix
    for i=1:Ndist,
        if strcmp(lower(dist{i}.type),'sample')
            disp(' ');
            disp(sprintf('Distribution %d',i));
            disp('Posterior Correlation Matrix:');
            corrcoef(dist{i}.P(:,dist{i}.ind)')
        end
    end
end

author: spm


text:
function sn = GetSn(Y, range_ff, method)
%% Estimate noise standard deviation

%% inputs:
%   Y: N X T matrix, fluorescence trace
%   range_ff : 1 x 2 vector, nonnegative, max value <= 0.5, range of frequency (x Nyquist rate) over which the spectrum is averaged
%   method: string, method of averaging: Mean, median, exponentiated mean of logvalues (default)

%% outputs:
%   sn: scalar, std of the noise

%% Authors: Pengcheng Zhou, Carnegie Mellon University, 2016
% adapted from the MATLAB implemention by Eftychios Pnevmatikakis and the
% Python implementation from Johannes Friedrich

%% References
% Pnevmatikakis E. et.al., Neuron 2016, Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data

%% input arguments
if ~exist('range_ff', 'var') || isempty(range_ff)
    range_ff = [.25, .5];
end
if ~exist('method', 'var') || isempty(method)
    method = 'logmexp';
end
if any(size(Y)==1)
    Y = reshape(Y, [], 1);
else
    Y = Y';
end

%% estimate the noise
[psdx, ff] = pwelch(double(Y), [],[],[], 1);
indf = and(ff>=range_ff(1), ff<=range_ff(2));
switch method
    case 'mean'
        sn=sqrt(mean(psdx(indf, :)/2));
    case 'median'
        sn=sqrt(median(psdx(indf,:)/2));
    case 'logmexp'
        sn = sqrt(exp(mean(log(psdx(indf,:)/2))));
    otherwise
        fprintf('wrong method! use logmexp instead.\n');
        sn = sqrt(exp(mean(log(psdx(indf,:)/2))));
end
sn = sn';

author: zhoupc


text:
function O = PetscOpenSocket(socketnumber)
%
%   O = PetscOpenSocket(socketnumber) - waits for a socket connection (from PETSc socket viewer)
%
%  This provides an object oriented interface to the PETSc provided Matlab routines sopen(), sread() and sclose()
%  allowing PETSc  MATLAB utilities like PetscBinaryRead.m to work cleanly with either binary
%  files or sockets
%
if nargin == 0
  S = struct('fd', sopen());
else
  S = struct('fd', sopen(socketnumber));
end
O = class(S,'PetscOpenSocket');

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function y = log_sum_exp( x, dim )

%LOG_SUM_EXP    log(sum(exp(x))).
%   LOG_SUM_EXP(X) = LOG(SUM(EXP(X)).
%
%   When used in a CVX model, LOG_SUM_EXP(X) causes CVX's successive
%   approximation method to be invoked, producing results exact to within
%   the tolerance of the solver. This is in contrast to LOGSUMEXP_SDP,
%   which uses a single SDP-representable global approximation.
%
%   If X is a matrix, LOGSUMEXP_SDP(X) will perform its computations
%   along each column of X. If X is an N-D array, LOGSUMEXP_SDP(X)
%   will perform its computations along the first dimension of size
%   other than 1. LOGSUMEXP_SDP(X,DIM) will perform its computations
%   along dimension DIM.
%
%   Disciplined convex programming information:
%       LOG_SUM_EXP(X) is convex and nondecreasing in X; therefore, X
%       must be convex (or affine).

error( nargchk( 1, 2, nargin ) );
if ~isreal( x ),
    error( 'Argument must be real.' );
end
y = exp( x );
if nargin == 1,
    y = sum( y );
else
    y = sum( y, dim );
end
y = log( y );

% 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:
%  Figure 10.78      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% fig10_78.m is a script for the LQR(LQG) design for the disk
% drive case study

clf;
wm= 5*pi;
z=.05;
GM = 4;
a=0.1;
numG= [1/(50*pi) 1];
denG=[1/(25*pi^2) 1/(50*pi) 1  0 0];
sysG=tf(numG,denG);
% form the system with y + a*ydot output
numGv=[0 numG]+.09*[numG 0];
sysGv=tf(numGv,denG);
% convert the design system to state form
[f,g,h,j]=ssdata(sysGv);
sysGvs=ss(f,g,h,j);
K=lqry(sysGvs,10000,1)
% form the state feedback
fc=f-g*K;
sysR=ss(fc,g,h,j);
[num,den]=tfdata(sysR,'v');
% remove the zero used to weight ydot
sysCLR=tf(numG,den);
[ac,b,c,d]=ssdata(sysCLR);
% normalize the gain for unity dc
k=c*inv(-ac)*b
b1=b/k;
sysCLR=ss(ac,b1,c,d);
t=0:.01:1.5;
[y,t]=step(sysCLR,t);
plot(t,y);
xlabel('Time (msec)');
ylabel('Amplitude');
grid;
title('Fig. 10.61 Step response of the LQR design for the disk');
%pause;
pc=eig(sysCLR);
pe=10*pc;
[f,g,h,j]=ssdata(sysG);
L=place(f',h',pe)';
nbar=1
a=[f-g*K-L*h];
b=[L -L];
c=-K;
d=[0 0];
sysF=ss(a,b,c,d);
sysOL=sysF*sysG;
sysH=tf(1,1);
sysCL=feedback(sysOL,sysH,2,1);
%
%grid
nicegrid

author: Sable


text:
function [MinPixVal, MaxPixVal]=max_min_pix_value(im_array);
% function [MinPixVal, MaxPixVal]=max_min_pix_value(im_array);
%
% function for calculating most often occuring Min and Max intensity values
% in images
% Perfromed by calculating means of pixels within the moving boxes of 5x5 pixel size throughout the
% image
% three most minimum values are averaged to get MinPixVal
% three most maximum values are averaged to get MaxPixVal
% input im_array:
% can be one image or array of images
%
% written by K.Artyushkova
% 102003

% Kateryna Artyushkova
% Postdoctoral Scientist
% Department of Chemical and Nuclear Engineering
% The University of New Mexico
% (505) 277-0750
% kartyush@unm.edu

[n,m,p]=size(im_array);
for k=1:p
    a=im_array(:,:,k);
    i=1:5:n-4;
    [e,f]=size(i);
    j=1:5:m-4;
    [g,h]=size(j);
    for I=1:f
        for J=1:h
            b(I,J)=mean(mean(a(i(I):i(I)+4,j(J):j(J)+4)));
        end
    end
    [r,p]=size(b);
    c=reshape(b,[r*p 1]);
    d=sort(c);
    MinPixVal(k)=mean(d(1:3));
    MaxPixVal(k)=mean(d((r*p-2):r*p));
end


text:
% Copyright 2018 Marc René Schädler
%
% This file is part of the mobile hearing aid prototype project
% The the mobile hearing aid prototype project 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.
%
% The mobile hearing aid prototype project 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 the mobile hearing aid prototype project. If not, see http://www.gnu.org/licenses/.

function signal = sinesweepphase(samples, frequency, phase)
%usage: signal = sinesweepphase(samples, frequency, phase)
%
% samples          - number of samples
% frequency        - normalized frequencies [0..1]
% phase            - phase [0..2*pi]

phase_diff_start = 2.*pi.*frequency(1);
phase_diff_stop = 2.*pi.*frequency(2);

phase_diff = logspace(log10(phase_diff_start),log10(phase_diff_stop),samples);

if nargin >= 4
  phase_diff(1) = phase_diff(1) + phase;
end

signal = sin(cumsum(phase_diff)).';
end

author: m-r-s

repo: hearingaid-prototype


text:
function matlab_install(SuiteSparse_path)
%  Matlab function to compile all the c-files to mex in the GPstuff toolbox.
%
%  Some of the sparse GP functionalities in the toolbox require
%  SuiteSparse toolbox by Tim Davis:
%    http://www.cise.ufl.edu/research/sparse/SuiteSparse/current/SuiteSparse/
%
%  This package includes the SuiteSparse version 4.4.2.
%
%  * To install without SuiteSparse run matlab_install
%  * To install with SuiteSparse run matlab_install('SuiteSparseOn')
%
%   The function matlab_install compiles the mex-files and prints on
%   the screen, which directories should be added to Matlab paths.

% Copyright (c) 2008-2012, 2016 Jarno Vanhatalo
% Copyright (c) 2016            Eero Siivola

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

    if nargin < 1
        SuiteSparse_path = [];
        fprintf('\n The path to the SuiteSparse package is not provided. \n')
        fprintf('\n Installing GPstuff without compactly supported covariance functions. \n')
        fprintf(' You are not able to use the following covariance functions:  \n')
        fprintf(' gpcf_ppcs0  \n gpcf_ppcs1 \n gpcf_ppcs2 \n gpcf_ppcs3 \n\n\n')
    elseif strcmp(SuiteSparse_path, 'SuiteSparseOn')
        cdir = pwd;
        cd SuiteSparse % Change this to correspond your own settings!
        SuiteSparse_path = path_spaces([pwd '/']);

        % Compile SuiteSparse
        fprintf('Compiling SuiteSparse. This may take a while \n \n')
        SuiteSparse_install(false);

        paths{1} = SuiteSparse_path;
        paths{2} = [SuiteSparse_path 'UMFPACK/MATLAB'];
        paths{3} = [SuiteSparse_path 'CHOLMOD/MATLAB'];
        paths{4} = [SuiteSparse_path 'AMD/MATLAB'];
        paths{5} = [SuiteSparse_path 'COLAMD/MATLAB'];
        paths{6} = [SuiteSparse_path 'CCOLAMD/MATLAB'];
        paths{7} = [SuiteSparse_path 'CAMD/MATLAB'];
        paths{8} = [SuiteSparse_path 'CXSparse/MATLAB/Demo'];
        paths{9} = [SuiteSparse_path 'CXSparse/MATLAB/CSparse'];
        paths{10} = [SuiteSparse_path 'LDL/MATLAB'];
        paths{11} = [SuiteSparse_path 'BTF/MATLAB'];
        paths{12} = [SuiteSparse_path 'KLU/MATLAB'];
        paths{13} = [SuiteSparse_path 'SPQR/MATLAB'];
        paths{14} = [SuiteSparse_path 'RBio/RBio'];
        paths{15} = [SuiteSparse_path 'MATLAB_Tools'];
        paths{16} = [SuiteSparse_path 'MATLAB_Tools/Factorize'];
        paths{17} = [SuiteSparse_path 'MATLAB_Tools/MESHND'];
        paths{18} = [SuiteSparse_path 'MATLAB_Tools/LINFACTOR'];
        paths{19} = [SuiteSparse_path 'MATLAB_Tools/find_components'];
        paths{20} = [SuiteSparse_path 'MATLAB_Tools/GEE'];
        paths{21} = [SuiteSparse_path 'MATLAB_Tools/shellgui'];
        paths{22} = [SuiteSparse_path 'MATLAB_Tools/waitmex'];
        paths{23} = [SuiteSparse_path 'MATLAB_Tools/spqr_rank'];
        paths{24} = [SuiteSparse_path 'MATLAB_Tools/spqr_rank/SJget'];
        paths{25} = [SuiteSparse_path 'MATLAB_Tools/UFcollection'];
        paths{26} = [SuiteSparse_path 'MATLAB_Tools/SSMULT'];
        paths{27} = [SuiteSparse_path 'MATLAB_Tools/dimacs10'];
        paths{28} = [SuiteSparse_path 'MATLAB_Tools/spok'];
        paths{29} = [SuiteSparse_path 'MATLAB_Tools/sparseinv'];

        cd(cdir)
        fprintf('Compiling GPstuff. This may take a while \n \n')
    else
        error('Unknown input argument. See help matlab_install for usage.')
    end

    % Go to diag/ and compile the mex-functions
    fprintf('\n Compiling files in diag.\n \n')
    cd('diag')
    diag_install
    cd('..')

    % Go to dist/ and compile the mex-functions
    fprintf('\n Compiling files in dist.\n \n')
    cd('dist')
    dist_install
    cd('..')

    % Go to gp/ and compile the mex-functions
    fprintf('\n Compiling files in gp.\n \n')
    cd('gp')
    gp_install(SuiteSparse_path)
    cd('..')

    % Go to mc/ and compile the mex-functions
    fprintf('\n Compiling files in mc. \n \n')
    cd('mc')
    mc_install
    cd('..')

    PP = pwd;
    S{1} = [PP '/diag'];
    S{2} = [PP '/dist'];
    S{3} = [PP '/gp'];
    S{4} = [PP '/mc'];
    S{5} = [PP '/misc'];
    S{6} = [PP '/optim'];
    S{7} = [PP '/tests'];

    fprintf ('\n The following paths have been added.  You may wish to add them\n') ;
    fprintf ('permanently, using the MATLAB pathtool command or copying the below\n') ;
    fprintf ('lines to your startup.m file. \n\n');
    for i = 1:length(S)
       addpath(S{i});
       fprintf ('addpath %s\n', S{i}) ;
    end

    if nargin==1
        fprintf ('\n')
        for k = 1:length (paths)
            fprintf ('addpath %s\n', paths {k}) ;
        end
    end
end

function path = path_spaces(path)
% Build correct path if path includes spaces

space_ind=strfind(path, ' ');
path=strrep(path, ' ', ''' ');
space_ind=space_ind+length(space_ind);
for i=1:length(space_ind)

  indd=strfind(path(space_ind(i):end), '/');
  path=[path(1:space_ind(i)+indd(1)-2) '''/' path(space_ind(i)+indd(1):end)];
  space_ind=space_ind+1;

end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function value = i4_bit_reverse ( i, n )

%*****************************************************************************80
%
%% I4_BIT_REVERSE reverses the bits in an I4.
%
%  Discussion:
%
%    An I4 is an integer value.
%
%  Example:
%
%       I      N  2^N     I4_BIT_REVERSE ( I, N )
%    ----    --------  -----------------------
%       0      0    1     0
%       1      0    1     1
%
%       0      3    8     0
%       1      3    8     4
%       2      3    8     2
%       3      3    8     6
%       4      3    8     1
%       5      3    8     5
%       6      3    8     3
%       7      3    8     7
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the integer to be bit reversed.
%    I should be nonnegative.  Normally I < 2^N.
%
%    Input, integer N, indicates the number of bits to
%    be reverse (N+1) or the base with respect to which the integer is to
%    be reversed (2^N).  N should be nonnegative.
%
%    Output, integer VALUE, the bit reversed value.
%
  i = round ( i );
  n = round ( n );

  if ( i < 0 )

    value = -1;

  elseif ( n < 0 )

    value = -1;

  else

    b = 2^n;
    j = mod ( i, b );

    value = 0;

    while ( 1 )

      if ( b == 1 )

        value = value + j;
        j = 0;
        break

      else

        if ( mod ( j, 2 ) == 1 )
          value = value + floor ( b / 2 );
          j = j - 1;
        end

        j = floor ( j / 2 );
        b = floor ( b / 2 );

      end

    end

  end

  return
end

author: johannesgerer


text:
function C = spm_mesh_isoline(M, T, t)
% Compute isolines on a triangular mesh
% FORMAT C = spm_mesh_isoline(M, T, t)
% M   - a GIfTI object or patch structure
% T   - [vx1] data vector
% t   - isovalue [Default: 0]
%
% C   - struct array of isolines with fields 'xdata', 'ydata', 'zdata' and
%       'isopen'
%__________________________________________________________________________
%
% M = gifti(fullfile(spm('Dir'),'canonical','cortex_20484.surf.gii'));
% M = export(M,'patch');
% M = spm_mesh_inflate(M);
% T = randn(size(M.vertices,1),1);
% T = spm_mesh_smooth(M,T,100);
% H = spm_mesh_render('Disp',M);
% H = spm_mesh_render('Overlay',H,T);
% hold on
% t = linspace(min(T),max(T),20);
% for i=1:numel(t)
%   C = spm_mesh_isoline(M,T,t(i));
%   for j=1:numel(C)
%     plot3(C(j).xdata,C(j).ydata,C(j).zdata,'k-');
%   end
% end
%__________________________________________________________________________
% Copyright (C) 2019 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_isoline.m 7618 2019-06-17 12:29:46Z guillaume $

if nargin < 3, t = 0; end

C = spm_mesh_contour(M,struct('T',T,'t',t));


text:
function output = callsparsecolo(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
x0      = interfacedata.x0;
ub      = interfacedata.ub;
lb      = interfacedata.lb;

% Bounded variables converted to constraints
if ~isempty(ub)
    [F_struc,K] = addStructureBounds(F_struc,K,ub,lb);
end

A = -F_struc(:,2:end)';
C = F_struc(:,1);
b = -c;

if options.savedebug
    ops = options.sparsecolo;
    save sparsecolodebug K A C b
end

ops = options.sparsecolo;
ops.SDPsolver = lower(interfacedata.solver.sdpsolver.tag);
ops.SDPAoptions = options.sdpa;
ops.sedumipar = options.sedumi;
ops.sdpt3OPTIONS = options.sdpt3;
ops.sdpt3OPTIONS.printlevel = options.verbose;
ops.printlevel = options.verbose;
if options.verbose==0
    ops.SDPAoptions.print = 'no';
else
    ops.SDPAoptions.print = 'display';
end

if options.savedebug
    save sparsecolodebug K A C b ops
end

if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end
solvertime = tic;
if options.verbose==0 % Sparsecolo does not run silent
    evalc('[x,y,infoCoLO,cliqueDomain,cliqueRange,LOP] = sparseCoLO(A,b,C,K,[],ops);');
else
    [x,y,infoCoLO,cliqueDomain,cliqueRange,LOP] = sparseCoLO(A,b,C,K,[],ops);
end
solvertime = toc(solvertime);

% Create YALMIP dual variable and slack
Dual = x;
Primal = y;

switch  lower(interfacedata.solver.sdpsolver.tag)
    case 'sdpt3'
        switch infoCoLO.SDPsolver.termcode
            case 0
                problem = 0; % No problems detected
            case {-1,-5}
                problem = 5; % Lack of progress
            case {-2,-3,-4,-7}
                problem = 4; % Numerical problems
            case -6
                problem = 3; % Maximum iterations exceeded
            case -10
                problem = 7; % YALMIP sent incorrect input to solver
            case 1
                problem = 2; % Dual feasibility
            case 2
                problem = 1; % Primal infeasibility
            otherwise
                problem = -1; % Unknown error
        end
    case 'sedumi'
        problem = sedumicode(infoCoLO.SDPsolver,options);

    case 'sdpa'
        switch (infoCoLO.SDPsolver.phasevalue)
            case 'pdOPT'
                problem = 0;
            case {'noINFO','pFEAS','dFEAS','pdFEAS'}
                problem = 3;
            case 'pFEAS_dINF'
                problem = 2;
            case 'pINF_dFEAS'
                problem = 1;
            case 'pUNBD'
                problem = 2;
            case 'dUNBD'
                problem = 1;
            case 'pdINF'
                problem = 12;
            otherwise
                problem = -1;
        end

end

if options.savesolveroutput
    solveroutput.obj = obj;
    solveroutput.X = X;
    solveroutput.y = y;
    solveroutput.Z = Z;
    solveroutput.info = info;
    solveroutput.runhist = runhist;
else
    solveroutput = [];
end

if options.savesolverinput
    solverinput.blk = blk;
    solverinput.A   = A;
    solverinput.C   = C;
    solverinput.b   = b;
    solverinput.X0   = [];
    solverinput.y0   = x0;
    solverinput.Z0   = [];
    solverinput.options   = options.sdpt3;
else
    solverinput = [];
end

% Standard interface
output = createOutputStructure(Primal,Dual,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

function problem = sedumicode(info,options)
temp = info.pinf;
pinf = info.dinf;
dinf = temp;

% Check for reported errors
if (pinf==0) &  (dinf==0)
    problem = 0; % No problems
end

% We can only report one error, use priorities
if (problem==0) & (pinf==1)
    problem = 1; % Primal infeasability
end

if (problem==0) & (dinf==1)
    problem = 2; % Dual infeasability
end

if (problem==0) & (info.numerr==1) | (info.numerr==2)
    problem = 4; %Numerical problems
end

if (problem==0) & (info.iter >= options.sedumi.maxiter)
    % Did we need exactly maxiter iterations to find optimum
    if (pinf==0) & (dinf==0) & (info.numerr==0)
        problem = 0; % Yes
    else
        problem = 3; % No, we are not optimal yet
    end
end

if (problem==0) & (info.feasratio<0.98)
    problem = 4;
end

% Fix for cases not really explained in documentation of sedumi?
if (abs(info.feasratio+1)<0.1) & (pinf==0) &  (dinf==0)
    problem = 1;
end
if (abs(info.feasratio+1)<0.1) & (pinf==0) &  (dinf==0) & (c'*y_s<-1e10)
    problem = 2;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP