text:
% load Human-GEM
model = importYaml('../../model/Human-GEM.yml');

% specify pairs of duplicate reactions, where reactions in first column are
% those that will be kept and those in the second column removed
rxns = {
    'MAR08971','MAR04306'
    'MAR00121','MAR00607'
    'MAR04411','MAR04911'
    'MAR00079','MAR06324'
    'MAR00101','MAR02199'
    'MAR02339','MAR00445'
    'MAR01728','MAR07441'
    'MAR01827','MAR07742'
    'MAR11237','MAR04455'
    'MAR11251','MAR04463'
    'MAR11252','MAR04478'
    'MAR11253','MAR04502'
    'MAR11254','MAR04517'
    'MAR11247','MAR04576'
    'MAR11249','MAR04609'
    'MAR11248','MAR04622'
    'MAR11255','MAR04626'
    'MAR11256','MAR04653'
    'MAR11257','MAR04659'
    'MAR11258','MAR04669'
    'MAR11260','MAR04713'
    'MAR11261','MAR04724'
    'MAR11262','MAR04747'
    'MAR11263','MAR04751'
    'MAR11264','MAR04753'
    'MAR11265','MAR04761'
    'MAR11238','MAR04794'
    'MAR11266','MAR04795'
    'MAR11267','MAR04798'
    'MAR11239','MAR04800'
    'MAR11268','MAR04801'
    'MAR11240','MAR06494'
    'MAR11241','MAR06497'
    'MAR11242','MAR06504'
    'MAR11243','MAR06541'
    'MAR11244','MAR06582'
    'MAR11246','MAR06590'
    'MAR11271','MAR06594'
    'MAR11245','MAR06706'
    'MAR11280','MAR06853'
    'MAR11278','MAR06856'
    'MAR11279','MAR06857'
    'MAR11259','MAR06858'
    'MAR09057','MAR07191'};

% verify that the reaction pairs have identical stoichiometry (either forward or reverse)
for i = 1:size(rxns, 1)
    indx = getIndexes(model, rxns(i,:), 'rxns');
    coeffs = full(model.S(any(model.S(:, indx) ~= 0, 2), indx));
    coeffs_norm = coeffs ./ max(abs(coeffs), [], 1);
    if ~isequal(coeffs_norm(:,1), coeffs_norm(:,2)) && ...
            ~isequal(-coeffs_norm(:,1), coeffs_norm(:,2))
        error('Reactions in row %u are not duplicates!', i);
    end
end

% normalize the stoichiometric coefficients of MAR08971
%  currently: 3 glucose-6-phosphate[c] + 3 NADP+[c] => 3 glucono-1,5-lactone-6-phosphate[c] + 3 H+[c] + 3 NADPH[c]
% normalized:   glucose-6-phosphate[c] +   NADP+[c] =>   glucono-1,5-lactone-6-phosphate[c] +   H+[c] +   NADPH[c]
rxn_indx = getIndexes(model, 'MAR08971', 'rxns');
model.S(:, rxn_indx) = model.S(:, rxn_indx) / 3;

% merge rxn annotations from the deleted duplicate rxn into the kept rxn
rxnAssocFile = '../../model/reactions.tsv';
rxnAssoc = importTsvFile(rxnAssocFile);
[~, rxn_indx] = ismember(rxns, rxnAssoc.rxns);
mergeFields = setdiff(fieldnames(rxnAssoc), {'rxns', 'rxnRetired', 'spontaneous'});
for i = 1:size(rxns, 1)
    for f = 1:numel(mergeFields)
        entry1 = strsplit(rxnAssoc.(mergeFields{f}){rxn_indx(i,1)}, ';');
        entry2 = strsplit(rxnAssoc.(mergeFields{f}){rxn_indx(i,2)}, ';');
        merged_entry = setdiff(union(entry1, entry2), {''});
        if ~isempty(merged_entry)
            rxnAssoc.(mergeFields{f}){rxn_indx(i,1)} = strjoin(merged_entry, ';');
        end
    end
end

% save annotations for to-be-deleted reactions to the deprecated rxns file
depRxnFile = '../../data/deprecatedIdentifiers/deprecatedReactions.tsv';
deprecRxns = importTsvFile(depRxnFile);
deprecRxnsTable = struct2table(deprecRxns);
rxnAssocTable = struct2table(rxnAssoc);
exportTsvFile([deprecRxnsTable; rxnAssocTable(rxn_indx(:,2), :)], depRxnFile);

% merge additional model fields from duplicate rxns with the kept rxns
mergeFields = {'eccodes'; 'rxnReferences'};
for i = 1:size(rxns, 1)
    for f = 1:numel(mergeFields)
        entry1 = strsplit(model.(mergeFields{f}){rxn_indx(i,1)}, ';');
        entry2 = strsplit(model.(mergeFields{f}){rxn_indx(i,2)}, ';');
        merged_entry = setdiff(union(entry1, entry2), {''});
        if ~isempty(merged_entry)
            model.(mergeFields{f}){rxn_indx(i,1)} = strjoin(merged_entry, ';');
        end
    end
end

% delete reactions from model and annotation file
model = removeReactions(model, rxns(:,2));
exportYaml(model, '../../model/Human-GEM.yml');
rxnAssocTable(rxn_indx(:,2), :) = [];
exportTsvFile(rxnAssocTable, rxnAssocFile);

author: SysBioChalmers

repo: Human-GEM


text:
function a = adaboost(C,hyper)

%  ADABOOST
%
%
% A=ADABOOST(H) implements basic boosting algorithm using max_margin linear hyperplanes.
%
%
% Hyperparameters, and their defaults
% kmax = 5                -- maximum number of weak learners
%
% Methods:
%  train, test
%
% Model:
%  child                  -- hyperplanes
% Example:
%
% % Use adaboost with 1-knn as weak learner and validate with 2 fold cross validation.
% c1=[2,0];
% c2=[-2,0];
% X1=randn(50,2)+repmat(c1,50,1);
% X2=randn(50,2)+repmat(c2,50,1);
%
% [r,a]=train(cv(adaboost(knn),'folds=2'),d);
%
% ========================================================================
% Reference : The boosting approach to machine learning: An overview
% Author    : Robert E. Schapire
% Link      : http://www.cs.princeton.edu/~schapire/uncompress-papers.cgi/msri.ps
% ========================================================================

  % model
  if (nargin <1)
      C=dualperceptron('margin=1');
  end
  a.alpha=[];
  a.child={};
  a.kmax=5;

  p=algorithm('adaboost');
  a= class(a,'adaboost',p);

  if nargin==2,
    eval_hyper;
  end;

  a.alpha=ones(a.kmax,1);

  for i=1:a.kmax
  	a.child{i}=C;
  end

  disp([num2str(a.kmax),' (', C.name, ') classifiers '])

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [cCombo,dataRecur]=getNextRucksackFill(param1,N)
%%GETNEXTRUCKSACKFILL Get the next combination of items with given positive
%               weights that will fit into a rucksack (knapsack, backpack)
%               with a capacity of N. This goes through all feasible
%               combinations (starting with putting nothing in the bag).
%               Going through all feasible combinations is not an efficient
%               approach to determining the set of items that can be placed
%               into the bag to maximize its weight.
%
%INPUTS: param1 To generate the first combination (which will be the empty
%               set - put nothing in the bag) this is w, an nX1 vector
%               holding the weights of each of the items that can be placed
%               in the bag. All weights>0, but need not be integer.
%               Otherwise, this is the dadaRecur output from the previous
%               call to this function.
%             N The capacity of the rucksack. N>0.
%
%OUTPUTS: cCombo The tX1 vector of indices of the items in w that form the
%                current combination of elements that fit into the
%                rucksack. An empty matrix is the first set (out nothing
%                in) and is also returned when the final value has been
%                passed.
%      dataRecur A data structure that can be passed to this function on
%                subsequent calls to get the next value in the sequence. If
%                an empty matrix is reurned, then the final value in the
%                sequence has been passed.
%
%This function implements Algorithm F of Section 7.2.1.3 of [1].
%
%EXAMPLE:
%This displays all of the combination indices for a simple problem.
% w=[1;1.5;5;3;1.5];
% N=6;
% [cCombo,dataRecur]=getNextRucksackFill(w,N);
% while(~isempty(dataRecur))
%     cCombo
%     [cCombo,dataRecur]=getNextRucksackFill(dataRecur);
% end
%
%REFERENCES:
%[1] D. E. Knuth, The Art of Computer Programming. Vol. 4A: Combinatorial
%    Algorithms, Part I, Boston: Addison-Wesley, 2011.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%If we are getting the first filling possibility (the empty sack).
if(nargin==2)
    w=param1;%The set of all-positive weights.
    [w,idx]=sort(w,'ascend');

    n=length(w);
    deltaJ=w(2:n)-w(1:(n-1));

    %The extra one at the beginning is for unassigned weights.
    c=zeros(n+1,1);

    %Step F1, Initialize
    t=0;
    c(0+1)=n;
    r=N;
    %The initial combination is the empty set.
    cCombo=[];
    dataRecur.w=w;
    dataRecur.idx=idx;
    dataRecur.deltaJ=deltaJ;
    dataRecur.c=c;
    dataRecur.r=r;
    dataRecur.t=t;
    return;
else
    dataRecur=param1;
    w=dataRecur.w;
    idx=dataRecur.idx;
    deltaJ=dataRecur.deltaJ;
    c=dataRecur.c;
    r=dataRecur.r;
    t=dataRecur.t;
end

%Step F3, try to add w(0).
if(c(t+1)>0&&r>=w(0+1))
    t=t+1;
    c(t+1)=0;
    r=r-w(0+1);

    %Step F2, visit the current combination.
    cCombo=idx(c(2:(t+1))+1);
    dataRecur.t=t;
    dataRecur.c=c;
    dataRecur.r=r;
    return;
end

while(1)
    if(t==0)%We have passed the last combination.
        cCombo=[];
        dataRecur=[];
        return;
    end

    if(c(t-1+1)>c(t+1)+1&&r>=deltaJ(c(t+1)+1))
        %Step F4
        c(t+1)=c(t+1)+1;
        r=r-deltaJ(c(t+1));
        %Step F2, visit the current combination.
        cCombo=idx(c(2:(t+1))+1);
        dataRecur.t=t;
        dataRecur.c=c;
        dataRecur.r=r;
        return
    else
        %Step F5
        r=r+w(c(t+1)+1);
        t=t-1;
    end
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function display(sAF,varargin)
% standard output

displayClass(sAF,inputname(1),varargin{:});

disp([' bandwidth: ' num2str(sAF.sF.bandwidth)]);

if sAF.sF.antipodal, disp(' antipodal: true'); end

end

author: mtex-toolbox

repo: mtex

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


text:
function Line=findlyap(MainHandles)
%FINDLYAP  determines the Lyapunov exponents and dimension
%
%    The alogrithm employed in this toolbox for determining Lyapunov
%    exponents is according to the algorithms proposed in
%
%    [1] A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano,
%        "Determining Lyapunov Exponents from a Time Series," Physica D,
%        Vol. 16, pp. 285-317, 1985.
%
%    [2] J. P. Eckmann and D. Ruelle, "Ergodic Theory of Chaos and Strange
%        Attractors," Rev. Mod. Phys., Vol. 57, pp. 617-656, 1985.
%
%    The algorithm given in [1] is used for first-order systems while
%    the QR-based algorithm proposed in [2] is applied for higher order
%    systems.

%   by Steve W. K. SIU, July 5, 1998.

%Print the results to the file every iteration
PrintStep=10;
%Display the results on the screen every iteration
DisplayStep=1;

%Clear the current axis objects
cla;
v=axis;
xc=(v(2)+v(1))/2;	%Center of the axis box [xc,y]
y=(v(4)+v(3))/2;
aW=v(2)-v(1);		%Axis width
x=xc-1/4*aW;

%Display the text "Loading... Please wait!"
th=text(x,y,'Loading... Please wait!','Color','r','FontSize',15,...
   	'FontAngle','italic');
drawnow;
%Clear the bottom text if any
set(MainHandles(1),'String','');

%Get the data stored in 'UserData' of the setting button
DATA=get(MainHandles(4),'UserData');
%Restore the data input by the user
output=DATA(1);		%Output checkbox: "on"=1; "off"=0
LEout=DATA(2);			%Output Lyapunov exponents to file: "yes"=1; "no"=0
LEprec=DATA(3);      %Precision of the output Lyapunov exponents
LDout=DATA(4);			%Output Lyapunov dimension to file: "yes"=1; "no"=0
LDprec=DATA(5);      %Precision of the Lyapunov dimension
IntMethod=DATA(6);	%Integration method: 1=Discrete map, 2=ODE45, 3=ODE23
							% 4=ODE113, 5=ODE23S, 6=ODE15S
InitialTime=DATA(7);	%Initial time
FinalTime=DATA(8);	%Final time
TimeStep=DATA(9);		%Time step
RelTol=DATA(10);		%Relative tolerance
AbsTol=DATA(11);		%Absolute tolerance
plot1=DATA(12);		%Plot imediately: 1="checked", 0="unchecked"
plot2=DATA(13);		%Plot according to specified iterations
ItrNum=DATA(14);		%No. of iteration for updating the plot
Color=DATA(15);		%Line Color option
							%Line color: 1="blue", 2="balck", 3="green", 4="red",
                     %5="yellow", 6="magenta", 7="cyan"
DiscardItr=DATA(16);	%Iterations to be discarded
UpdateStepNum=DATA(17);	%Lyapunov exponents updating steps
linODEnum=DATA(18);	%No. of linearized ODEs
ic=DATA(19:length(DATA));	%Initial conditions

%Get the output file and ODE function names stored in
%'UserData' of the "Start" button.
%Handles(2) is the handle of "Start" button
NAMES=get(MainHandles(2),'UserData');
OutputFile=rmspace(NAMES(1,:));
odefun=rmspace(NAMES(2,:));

%Construct a look-up table for the line colors
COLORS='bkgrymc';
%Map the "Line color" pop-up menu position to the look-up table
LineColor=COLORS(Color);

%Construct a look-up table for integration methods
Methods=char('Discrete map', 'ode45','ode23','ode113','ode23s','ode15s');
ODEsolver=strcat(Methods(IntMethod,:));

%Dimension of the linearized system (total: d x d ODEs)
d=sqrt(linODEnum);
%Initial conditions for the linearized ODEs
Q0=eye(d);
IC=[ic(:);Q0(:)];
ICnum=length(IC);		%Total no. of initial coniditions
%One iteration: Duration for updating the LEs
Iteration=UpdateStepNum*TimeStep;
DiscardTime=DiscardItr*Iteration+InitialTime;

%MATLAB's ODE functions will give the intermediate solutions if
%the duration between the initial time and the final time is only
%one time step, this will slow down the whole iteration process.
%To avoid this, reduce the time step by half.
if (UpdateStepNum==1 & IntMethod~=3)
   TimeStep=TimeStep/2;
end

T1=InitialTime;
T2=T1+Iteration;
TSpan=[T1:TimeStep:T2];
%Absolute tolerance of each components is set to the same value
options=odeset('RelTol',RelTol,'AbsTol',ones(1,ICnum)*AbsTol);

%Initialize variables
n=0;			%Iteration counter
k=0;			%Effective iteration counter
				% (discarded iterations are not counted)
h=1;			%No. of line handles sets (1 set = d line handles)
delLine=0;  %Indicator for deleting the drawn lines
Sum=zeros(1,d);
xData=[];
yData=[];
Line=[];
bufferSize=10000; %Max. no. of data can be stored in the buffer before creating a new line.
if ( output & (LEout | LDout) )
   %If the output file cannot be opened, warn the user
   msg=['Unable to open "' sprintf(OutputFile) '".'];
   Warn='errordlg(msg,''ERROR'',''replace''); problem=1;';
   eval('fid=fopen(OutputFile,''wt''); problem=0;',Warn)
   %Construct a look-up table for precision format
   Prec=char('%.4f','%.6f','%.8f','%.10f','%.12f');
   %Map the pop-up menu position to its corresponding format
   LEprecision=strcat(Prec(LEprec,:));
   LDprecision=strcat(Prec(LDprec,:));
   if ~problem
      fprintf(fid,'Time');
      if LEout
         for i=1:d
            Str1=['\tLE%d'];
            fprintf(fid,Str1,i);
         end
      end
      if LDout
         Str2=['\tLD'];
         fprintf(fid,Str2);
      end
      fprintf(fid,'\n');
   end
else
   problem=0;
end

%Start the stop watch
tic;
%Get the state of the "Stop" button
stop=get(MainHandles(3),'UserData');
if DiscardTime>0
   %Display the text "Discarding transient steps..."
   set(MainHandles(1),'String','Discarding transient steps...');
end

A=[];

%String that contains the integration command
IntegrationStr=['[t,X]=',ODEsolver,'(odefun,TSpan,IC,options);'];
%Main loop
while (~stop & ~problem)
   n=n+1;
   %Integration
   if IntMethod>1
      eval(IntegrationStr);
   else	%If it is a discrete map
      for i=1:UpdateStepNum
         X(i,:)=(feval(odefun,IC))';
      end
   end
   [rX,cX]=size(X);
   L=cX-linODEnum;      %No. of initial conditions for
                        %the original system
   for i=1:d
      m1=L+1+(i-1)*d;
      m2=m1+d-1;
      A(:,i)=(X(rX,m1:m2))';
   end
   %QR decomposition
   if d>1
      %The algorithm for 1st-order system doesn't require
      %QR decomposition
      [Q,R]=qr(A);
      if T2>DiscardTime
         Q0=Q;
      else
         Q0=eye(d);
      end
   else
      R=A;
   end

   %Delete the text "Loading...Please wait!"
   %before the first iteration
   if n==1
      delete(th);
      drawnow;
      %Display the final time
      set(MainHandles(11),'String',FinalTime);
   end

  %Any zero diagonal element will cause overflow
  %in the following calculation, so discard this step.
   permission=1;
   for i=1:d
      if R(i,i)==0
         permission=0;
         break;
      end
   end
  %To determine the Lyapunov exponents
   if (T2>DiscardTime & permission)
      k=k+1;
      T=k*Iteration;
      TT=n*Iteration+InitialTime;
      %There are d Lyapunov exponents
      Sum=Sum+log(abs(diag(R))');
      lambda=Sum/T;

      %Sort the Lyapunov exponents in descenting order
      Lambda=fliplr(sort(lambda));
      %To calculate the Lyapunov dimension (or Kaplan-Yorke dimension)
      LESum=Lambda(1);
      LD=0;
      if (d>1 & Lambda(1)>0)
         for N=1:d-1
            if Lambda(N+1)~=0
               LD=N+LESum/abs(Lambda(N+1));
               LESum=LESum+Lambda(N+1);
               if LESum<0
                  break;
               end
            end
         end
      end
      %Store the [x,y] data for plotting
      [rxD,cxD]=size(xData);
      [ryD,cyD]=size(yData);
      if rxD<=bufferSize
         xData=[xData;TT];
         yData=[yData;lambda];
      else
         %When the buffers are full, refresh them
         %Max. size of buffers = 10000 data
         xData=[xData(rxD);TT];
         yData=[yData(ryD,:);lambda];
         h=h+1;		%add one set of line handles
         delLine=0;	%After refreshing the buffers, the
                     % previous drawn lines must not be deleted
      end

      if ( output & ~problem & (LEout | LDout) & rem(k,PrintStep)==0)
         fprintf(fid,'%.2f',TT);
         if LEout
            for i=1:d
               Str1=['\t',LEprecision];
               fprintf(fid,Str1,Lambda(i));
            end
         end
         if LDout
            Str2=['\t',LDprecision];
            fprintf(fid,Str2,LD);
         end
         fprintf(fid,'\n');
      end

      %Draw lines immediately if "update the plot immediately" was chosen.
      if (plot1==1 | plot2==1 & rem(k,ItrNum)==0)
         if delLine
            %Clear the previous drawn line if any
            delete(Line(h,:));
         end
      	%Draw d lines
         for i=1:d
            %Set "Erase Mode" to "none" for increasing speed (less refresh)
            Line(h,i)=line('EraseMode','none','Color',LineColor,...
                      'xData',xData,'yData',yData(:,i));
            %Force MATLAB to draw immediately
            drawnow;
         end
         delLine=1;	%Set a flag to indicate that the lines now can be deleted
      end
      %Display the calculated Lyapunov exponents
      if rem(k,DisplayStep)==0
         set(MainHandles(1),'String',[num2str(Lambda),blanks(3),'( ',num2str(LD),' )']);
      end
   end

   %To see whether "Stop" is pressed
   stop=get(MainHandles(3),'UserData');

   %Display current time, and time used
   set(MainHandles(10),'String',num2str(round(T2)));		%Current time
   set(MainHandles(12),'String',num2str(round(toc)));		%Used time
   drawnow;

   %If calculation is finished or "stop" button is pressed, exit the loop.
   if (stop | T2>=FinalTime)
      %Reset the "Erase mode" to normal
      set(Line,'EraseMode','normal');
      %Show the final results (for making sure the final results being shown if DisplayStep>1)
      if (T2>DiscardTime & permission)
         set(MainHandles(1),'String',[num2str(Lambda),blanks(3),'( ',num2str(LD),' )']);
      end
      break;
   end
   %Update the initial conditions and time span for the next iteration
   if IntMethod>1
      ic=X(rX,1:L);
      T1=T1+Iteration;
      T2=T2+Iteration;
      TSpan=[T1:TimeStep:T2];
   else %For discrete map
      ic=X(UpdateStepNum,1:L);
      T2=T2+Iteration;
   end
   IC=[ic(:);Q0(:)];
end		%End of main loop

if T2>=FinalTime
   set(MainHandles([2,4,13,15]),'Enable','On');
   set(MainHandles(3),'Enable','Off');
end
if (output & ~problem)
   fclose(fid);
end

%----------------Subroutine-----------------------------------
function outStr=rmspace(inStr)
%RMSPACE		Function for removing the beginning and ending
%				spaces of a string

%Remove spaces at the end of the string
outStr=strcat(inStr);
%Delete spaces at the beginning of the string
if ~isempty(outStr)
   while isspace(outStr(1))
   	outStr=outStr(2:length(outStr));
   end
end

author: Sable

repo: mcbench-benchmarks


text:
function [V,D,Ann] = eign(A1, A2,A3, args)
% Find the principal eigenvalues and eigenvectors of a matrix with Nystr�m's low rank approximation method
%
% >>  D     = eign(A, nb)
% >> [V, D] = eign(A, nb)
%
% In the case of using this method for low rank approximation and
% decomposing the kernel matrix, one can call the function without
% explicit construction of the matrix A.
%
% >>  D     = eign(X, kernel, kernel_par, nb)
% >> [V, D] = eign(X, kernel, kernel_par, nb)
%
%
% Full syntax
%  (We denote the size of positive definite matrix A with a*a.)
%
%     1. Given the full matrix:
%
% >>  D    = eign(A,nb)
% >> [V,D] = eign(A,nb)
%
%
%       Outputs
%         V(*)  : a x nb matrix with estimated principal eigenvectors of A
%         D     : nb x 1 vector with principal estimated eigenvalues of A
%       Inputs
%         A     : a*a positive definite symmetric matrix
%         nb(*) : Number of approximated principal eigenvalues/eigenvectors
%
%
%     2. Given the function to calculate the matrix elements:
%
% >>  D = eign(X, kernel, kernel_par, nb)
% >> [V,D] = eign(X, kernel, kernel_par, nb)
%
%       Outputs
%         V(*)       : a x nb matrix with estimated principal eigenvectors of A
%         D          : nb x 1 vector with estimated principal eigenvalues of A
%       Inputs
%         X          : N x d matrix with the training data
%         kernel     : Kernel type (e.g. 'RBF_kernel')
%         kernel_par : Kernel parameter (bandwidth in the case of the 'RBF_kernel')
%         nb(*)      : Number of eigenvalues/eigenvectors used in the eigenvalue decomposition approximation
%
% See also:
%   eig, eigs, kpca, bay_lssvm

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

% AFUN?
 if nargin~=2
  X = A1;
  kernel = A2;
  kernel_par = A3;
  N = size(X,1);

  %eval(['if args<1, error(''strict positive number of eigenvalues required;'');else n = args; end;'],...
  %'n=min(6,ceil(N*.75));');
    if args<1
      error('Strict positive number of eigenvalues required');
  else
      n = args;
  end;

  % random sampling
  s = randperm(N); sr=s(n+1:end); s=s(1:n);
  %s = ceil(1:(N-1)/(n-1):N); s = s(1:n);

  ANn = zeros(n,N);
  ANn = kernel_matrix(X,kernel,kernel_par,X(s,:));
  Ann = ANn(s,:);

  % centering of matrix
  Zc = eye(n) - 1/n;
  %ZC = eye(N) - 1/N;
  Ann = Zc*Ann*Zc;
  %ANn = ZC*ANn*Zc;
  ANn = (Zc*(ANn*Zc)')';
else
  A = A1;
  N = size(A,1);

  %eval(['if args<1, error(''strict positive number of eigenvalues required;'');else n = args; end;'],...
  %'n=min(6,ceil(N*.75));');

  if A2<1
      error('Strict positive number of eigenvalues required');
  else
      n = A2;
  end;

  % random sampling
  s = randperm(N); sr=s(n+1:end); s=s(1:n);
  %s = ceil(1:(N-1)/(n-1):N); s = s(1:n);

  ANn = A(:,s);
  Ann = A(s,s);
end

%
% compute eigenvalues en vectors of low rank approximation
%
[Vn,Dn] = eig(Ann);Dn = diag(Dn);

%
% select only relevant eigenvalues and sort
% (only largest eigenvectors are orthogonal)
%
[Dn,peff] = sort(Dn(find(Dn>1000*eps)));
Dn = Dn(end:-1:1); peff = peff(end:-1:1);
Vn = Vn(:,peff);

%
% Nystrom correction
%
D = (N/n).*Dn;

%
% eigenvectoren correctie
%
if nargout>1,
  %V = zeros(n,length(peff));
  for i=1:length(peff),
    V(:,i) = (sqrt(n/N)/Dn(peff(i)))*ANn*Vn(:,peff(i));
  end

  %
  % correction of found eigenvectors:
  % orthogonal and unit length
  %

  % svd
  %[V,D2,ff] = svd(V*diag(D.^.5));   V = V(:,1:length(peff));

  % gram schmidt
  %[V,r] = gramschmidt2(V);
  %D = D.*r;

  %D = diag(D.^2);
else
  V=D;
end

author: peterhcharlton

repo: RRest


text:
function [Population,Archive,Fitness] = EnvironmentalSelection(Population,N)
% The environmental selection of SPEA2

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

    %% Calculate the fitness of each solution
    Fitness = CalFitness(Population.objs);

    %% Environmental selection
    Next    = Fitness < 1;
    Archive = Population(Next);
    if sum(Next) < N
        [~,Rank] = sort(Fitness);
        Next(Rank(1:N)) = true;
    elseif sum(Next) > N
        Del  = Truncation(Population(Next).objs,sum(Next)-N);
        Temp = find(Next);
        Next(Temp(Del)) = false;
        Archive = Population(Next);
    end
    % Population for Next generation
    Population = Population(Next);
    Fitness    = Fitness(Next);
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    %% Truncation
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,size(PopObj,1));
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function ang= getRelativeRigidOrientation(obj, frame, R, p)
    % Returns the symbolic representation of the Euler angles of a
    % rigid link.
    %
    % Parameters:
    % frame: the list of coordinate frame of the point
    % @type cell
    % R: ang is computed as the relative Euler angles from frame to R
    % @type matrix
    % p: the offset of the point from the origin of the frame
    % @type matrix
    %
    %
    % Return values:
    %  ang: the 3-dimensional Euler angles (roll,pitch,yaw) vector of the
    %  CoM of the system @type SymExpression
    %
    %
    % @note Syntax for ont point
    %
    % >> getEulerAngles(obj,pf,offset)
    %
    % @note Syntax for multiple points (offset should be np*3 matrix)
    %
    % >> getEulerAngles(obj,pfarray, offset)

    % the number of points (one less than the nargin)
    if nargin < 4
        p = [];
    end
    if nargin < 3
        R = eye(3);
    end
    c_str = getTwists(frame, p);
    c_str{1}.R = R;

    ang = eval_math_fun('ComputeRelativeRigidOrientation', c_str);

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
% test minq.m, minqdef.m, minqsep.m

disp('test of minq.m')

n=7;
m=10;
p=0.8;		% approx. fraction of activities and equalities
newdata=1;	% new data?

if newdata,
  % create random data satisfying the KKT conditions
  A=rand(m,n);
  c=rand(n,1);
  d=rand(n,1);
  ydes=rand(m,1)-p;
  act=( ydes>p );		% indices of nonactive inequalities
  ydes(act)=0*ydes(act);
  eq=( ydes<0 );		% indices of equations
  xdes=(A'*ydes-c)./d;
  res=rand(m,1);
  res(~act)=0*res(~act);
  b=A*xdes-res;
  save minq_data
 prt=0;
else
  load minq_data
  warning debug
  prt=1;
end;

for cas=1:2,
  if cas==1,
    disp('test of minqsep.m');
    [x,y,ier]=minqsep(c,d,A,b,eq,prt);
  else
    disp('test of minqdef.m');
    [x,y,ier]=minqdef(c,diag(d),A,b,eq,prt);
  end;

  compldes=[ydes,A*xdes-b];
  compl=[y,A*x-b];
  act_eq=[act,eq]'
  ydif=[ydes,y]
  disp('usually two equal columns');
  disp('but sometimes the dual is not unique')
  xdif=[xdes,x]
  disp('xdif should have two equal columns')
  ier
  if cas==2, break; end;
  cont=input('enter return (next test) or 0 (quit)>');
  if isempty(cont),	% continue
  elseif cont==0,      return;
  end;
end;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function H = hankel(c,r)

%Disciplined convex/geometric programming information for HANKEL:
%   HANKEL imposes no convexity restrictions on its arguments. Instead
%   of using the HANKEL function, however, consider creating a matrix
%   variable using the 'hankel' or 'upper_hankel' keyword; e.g.
%       variable X(5,5) hankel;
%       variable Y(4,4) upper_hankel;

%
% Check arguments
%

narginchk(1,2);
if nargin < 2,
    r = zeros(size(c));
else
    temp = cvx_subsref( r, 1 ) - cvx_subsref( c, numel(c) );
    if ~cvx_isnonzero( temp ),
        warning('MATLAB:hankel:AntiDiagonalConflict',['Last element of ' ...
               'input column does not match first element of input row. ' ...
               '\n         Column wins anti-diagonal conflict.'])
    end
end

%
% Compute indices and construct data vector
%

r  = vec( r );
c  = vec( c );
nc = length( c );
nr = length( r );
x  = [ c ; cvx_subsref( r, 2 : nr, 1 ) ];

%
% Construct matrix
%

cidx = ( 1 : nc )';
ridx = 0 : nr - 1;
H    = cidx(:,ones(nr,1)) + ridx(ones(nc,1),:);
H    = reshape( cvx_subsref( x, H( : ) ), size( H ) );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% Frequency response Magnet Driver circuit network using DS method
%  and G array.
% File:  c:\M_files\short_updates\UsingGarray1.m
% 02/10/07
% See Word file UsingGarray1.doc for background explanations.
clear;clc;
% unit suffixes
u=1e-6;K=1e3;m=1e-3;u=1e-6;p=1e-12;
% component values
R1=5*K;R2=1*K;R3=10*K;
C1=0.1*u;C2=400*p;
Ein=1; % Unity input for (normalized) transfer function
N=2; % Number of capacitors = order of circuit
%
% Form W, Q, S, and P arrays:
%
W=[1+R3/R1 1+R3/R1;R2 0];Q=[0 -1/R1;-1 1];S=[Ein/R1;0];P=diag([C1 C2]);
%
% Get A, B, D, & E arrays:
%
C=inv(W*P);A=C*Q;B=C*S;
% D & E from the first output equation:
D1=[0 R1/(R1+R3)];E1=R3/(R1+R3);
% D & E from the second (easier) output equation:
D2=[0 1];E2=0;
%
% Now we can form the G array = F*P as follows.
%
F=[R3 R3]; % R3 is the coefficient of both iC1 and iC2 in Vo2
% which are derivative terms (iC=dVc/dt, etc.)
G=F*P;
%
% * * * * * * * * * * * * Frequency response * * * * * * * * * * * *
%
I=eye(N); % identity matrix
%
% Log frequency sweep from BF to BF+ND
%
BF=0;ND=4;PD=25;NP=ND*PD+1;Fr=logspace(BF,BF+ND,NP);
for i=1:NP
   s=2*pi*Fr(i)*j; % j = sqrt(-1)
   stm=(s*I-A)\B;v1=abs(D1*stm+E1);Vo1(i)=20*log10(v1);
   v2=abs((D2+G*A)*stm+E2+G*B);Vo2(i)=20*log10(v2);
end
%
% Plot frequency (Vo1 & Vo2 separated by 0.2 dBV to ovoid overlay)
%
h=plot(log10(Fr),Vo1-0.2,'b',log10(Fr),Vo2,'r');
set(h,'LineWidth',2);
grid on;
axis auto
ylabel('dBV');title('Magnet Driver Output');
xlabel('Log Freq(Hz)');
legend('Vo1 - 0.2','Vo2 using G array');
figure(1);
%


text:
function plotCircle3D(center,normal,radius)

theta=0:0.01:2*pi;
v=null(normal);
points=repmat(center',1,size(theta,2))+radius*(v(:,1)*cos(theta)+v(:,2)*sin(theta));
plot3(points(1,:),points(2,:),points(3,:),'r-');

end


text:
function [HB, NL, level] = HBstructure3(elem,N0,HBmesh)
%% HBSTRUCTURE3 reconstructs a hierachical structure of a 3-D mesh.
%
% HB(:,2:3) are two parent nodes of the node HB(:,1):
%             HB(:,2) --- HB(:,1) --- HB(:,3)
% NL records the range of indices in each level. The indices of nodes in
% the k-th level is given by NL(k)+1:NL(k+1).
%
% See also HBstructure3
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('HBmesh','var'), HBmesh = []; end
N = max(elem(:));
HB = zeros(N,3);
% level = max(min(round(log2(N/N0)),16),2); % at least two level
level = 20;
NL(level+1) = N;
for k = level: -1 : 2
    [elem,HBmesh,newHB] = uniformcoarsen3red(elem,HBmesh);  % try coasen red refinement
    if isempty(newHB) && ~isempty(HBmesh)     % then try coarsen bisection
        [elem,HBmesh,newHB] = uniformcoarsen3(elem,HBmesh);
    end
    if (isempty(newHB)) || (size(elem,1)< 2*N0)
    % no nodes are removed or it reaches the coarsest level
        NL = NL(k:end);
        break;
    end
    NL(k) = NL(k+1) - size(newHB,1);
    HB(NL(k)+1:NL(k+1),1:3) = newHB(:,1:3);
end
level = length(NL)-1;

repo: ifem


text:
function img = ldrimread(filename, bDouble)
%
%       img = ldrimread(filename, bDouble)
%
%
%        Input:
%           -filename: the name of the file to open
%
%        Output:
%           -img: the opened image
%
%     Copyright (C) 2011-15  Francesco Banterle
%
%     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/>.
%

if(~exist('filename', 'var'))
    error('A filename with extension needs to be passed as input!');
end

if(~exist('bDouble', 'var'))
    bDouble = 1;
end

img = [];

try
    image_info = imfinfo(filename);

    if((image_info.BitDepth == 24) || (image_info.BitDepth == 8))
        if(bDouble)
            img = double(imread(filename)) / 255.0;
        else
            img = single(imread(filename)) / 255.0;
        end
    end

    if(image_info.BitDepth == 48)
        if(bDouble)
            img = double(imread(filename)) / 65535.0;
        else
            img = single(imread(filename)) / 65535.0;
        end
    end

catch err
    disp(['Error in reading ', filename]);
end

end

author: banterle


text:
%% Import from VPSC
%
%%
% <https://public.lanl.gov/lebenso/ VPSC> is a crystal plasticity code
% originally written by Ricardo Lebensohn and Carlos Tome from Los Alamos
% National Laboratory - USA.
%
% Original code can be requested to lebenso@lanl.gov
%
% https://public.lanl.gov/lebenso/
%
%% Import the orientations generated by VPSC
%
% Running a simulation in VPSC ussually results in an output file
% |TEX_PH1.OUT| which contains multiple sets of orientations for different
% strain levels. As these files does not contain any information on the
% crystal symmetry we need to specify it first

cs = crystalSymmetry('222', [4.762 10.225 5.994],'mineral', 'olivine');

%%
% In the next step the orientations are imported and converted into a list
% of <SO3Fun.SO3Fun.html ODFs> using the command <ODF.load.html |ODF.load|>.

% put in here the path to the VPSC output files
path2file = [mtexDataPath filesep 'VPSC'];

odf = ODF.load([path2file filesep 'TEX_PH1.OUT'],'halfwidth',10*degree,'CS',cs)

%%
% The individuel ODFs can be accessed by |odf{id}|

% lets plot the second ODF
plotSection(odf{2},'sigma','figSize','normal')

%%
% The information about the strain are stored as additional properties
% within each ODF variable

odf{1}.opt

%% Compare pole figures during deformation
%
% Next we examine the evaluation of the ODF during the deformation by
% plotting strain depended pole figures.

% define some crystal directions
h = Miller({1,0,0},{0,1,0},{0,0,1},cs,'uvw');

% generate some figure
fig = newMtexFigure('layout',[4,3],'figSize','huge');
subSet = 1:4;

% plot pole figures for different strain steps
for n = subSet
  nextAxis
  plotPDF(odf{n},h,'lower','contourf','doNotDraw');
  ylabel(fig.children(end-2),['\epsilon = ',xnum2str(odf{n}.opt.strain)]);
end
setColorRange('equal')
mtexColorbar

%% Visualize slip system activity
%
% Alongside with the orientation data VPSC also outputs a file
% |ACT_PH1.OUT| which contains the activity of the different slip systems
% during the deformation. Lets read this file as a table

ACT = readtable([path2file filesep 'ACT_PH1.OUT'],'FileType','text')

%%
% and plot the slip activity with respect to the strain for the different
% modes

% loop though the columns MODE1 ... MOD11
close all
for n = 3: size(ACT,2)

  % perform the plotting
  plot(ACT.STRAIN, table2array(ACT(:,n)),'linewidth',2,...
    'DisplayName',['Slip mode ',num2str(n-2)])
  hold on;
end
hold off

% some styling
xlabel('Strain');
ylabel('Slip activity');
legend('show','location','NorthEastOutside');

set(gca,'Ylim',[-0.005 1])
set(gcf,'MenuBar','none','units','normalized','position',[0.25 0.25 0.5 0.5]);

%for only one mode plot, e.g.,mode 3: cs = csapi(STRAIN,MODE{3});fnplt(cs,3,'color','b');hold off;


text:
function [ x, know ] = p05_sol ( m, know )

%*****************************************************************************80
%
%% P05_SOL returns known solutions for problem 5.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Harald Niederreiter, Kevin McCurley,
%    Optimization of functions by quasi-random search methods,
%    Computing,
%    Volume 22, Number 2, 1979, pages 119-123.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input/output, integer KNOW.
%    On input, KNOW is 0, or the index of the previously returned solution.
%    On output, KNOW is 0 if there are no more solutions, or it is the
%    index of the next solution.
%
%    Output, real X(M), the solution.
%
  if ( know == 0 )
    know = 1;
    x(1:m,1) = [ ...
       3.0 / 11.0, ...
       6.0 / 13.0, ...
      12.0 / 23.0, ...
       8.0 / 37.0 ]';
  else
    know = 0;
    x = zeros ( m, 1 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%
% Smooth point-set registration method using neighboring constraints
% -------------------------------------------------------------------
%
% Authors: Gerard Sanromà, René Alquézar and Francesc Serratosa
%
% Contact: gsanorma@gmail.com
% Date: 15/02/2012
%
% Centers and rescales input graphs according to the coefficients of the
% match matrix. It follows the ansatz by Softassign Procrustes.
%
% Input
%   g_M,g_D: 2D coordinates of the model and data graph nodes
%   S: matrix of matching coefficients
%
% Output
%   g_Mc,g_Dc: centered and rescaled coordinates of the model and data
%              graphs' nodes
%

function [g_Mc g_Dc] = center_rescale_weighted(g_M,g_D,S)
% Softassign Procrustes (Rangarajan)
% Update centroids and variances (eq.(5))
m_M = sum(g_M.*repmat(sum(S)',1,2))/(sum(S(:))+eps);
m_D = sum(g_D.*repmat(sum(S,2),1,2))/(sum(S(:))+eps);
r_M = g_M - repmat(m_M,size(g_M,1),1);
r_D = g_D - repmat(m_D,size(g_D,1),1);
v_M = sum(diag(r_M*r_M').*sum(S)');%/sum(S(:));
v_D = sum(diag(r_D*r_D').*sum(S,2));%/sum(S(:));
% Center and rescale
g_Mc = r_M/(sqrt(v_M)+eps);
g_Dc = r_D/(sqrt(v_D)+eps);


text:
function [recalls, allRecalls]= testCore(db, qFeat, dbFeat, varargin)
    opts= struct(...
        'nTestSample', inf, ...
        'recallNs', [1:5, 10:5:100], ...
        'printN', 10 ...
        );
    opts= vl_argparse(opts, varargin);

    searcherRAW_= @(iQuery, nTop) rawNnSearch(qFeat(:,iQuery), dbFeat, nTop);
    if ismethod(db, 'nnSearchPostprocess')
        searcherRAW= @(iQuery, nTop) db.nnSearchPostprocess(searcherRAW_, iQuery, nTop);
    else
        searcherRAW= searcherRAW_;
    end
    [res, recalls]= recallAtN( searcherRAW, db.numQueries,  @(iQuery, iDb) db.isPosQ(iQuery, iDb), opts.recallNs, opts.printN, opts.nTestSample );

    allRecalls= recalls;
    recalls= mean( allRecalls, 1 )';

end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
clc, clear
c=[3 8 2 10 3;8 7 2 9 7;6 4 2 7 5
   8 4 2 3 5;9 10 6 9 10];
c=c(:); a=zeros(10,25); intcon=1:25;
for i=1:5
   a(i,(i-1)*5+1:5*i)=1;
   a(5+i,i:5:25)=1;
end
b=ones(10,1); lb=zeros(25,1); ub=ones(25,1);
x=intlinprog(c,intcon,[],[],a,b,lb,ub);
x=reshape(x,[5,5])

author: Eurus-Holmes


text:
classdef PTKSaveCoronalAnalysisResults < PTKPlugin
    % PTKSaveCoronalAnalysisResults. Plugin for performing analysis of density using bins
    % along the anterior-posterior axis
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKSaveCoronalAnalysisResults divides the cranial-caudal axis into bins and
    %     performs analysis of the tissue density, air/tissue fraction and
    %     emphysema percentage in each bin.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Coronal analysis'
        ToolTip = 'Performs density analysis in bins along the anterior-posterior axis'
        Category = 'Slice-by-slice'
        Mode = 'Analysis'

        Context = PTKContextSet.LungROI
        AllowResultsToBeCached = true
        AlwaysRunPlugin = true
        PluginType = 'DoNothing'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '2'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = false
        Location = 31
        Icon = 'cor_slices.png'
    end

    methods (Static)
        function results = RunPlugin(dataset, context, reporting)

            % Generate the results over the lungs and lobes
            contexts = {PTKContextSet.Lungs, PTKContextSet.SingleLung, PTKContextSet.Lobe};
            results = dataset.GetResult('PTKCoronalAnalysis', contexts);

            % Convert the results into a MimResultsTable
            image_info = dataset.GetImageInfo;
            uid = image_info.ImageUid;
            patient_name = dataset.GetPatientName;
            table = PTKConvertMetricsToTable(results, patient_name, uid, reporting);

            % Save the results table as a series of CSV files
            dataset.SaveTableAsCSV('PTKSaveCoronalAnalysisResults', 'Coronal analysis', 'CoronalResults', 'Density analysis in bins along the anterior-posterior axis', table, MimResultsTable.ContextDim, MimResultsTable.SliceNumberDim, MimResultsTable.MetricDim, []);

            % Generate graphs of the results
            y_label = 'Distance along coronal axis (%)';

            context_list_both_lungs = [PTKContext.Lungs];
            PTKSaveCoronalAnalysisResults.DrawGraphAndSave(dataset, table, y_label, context_list_both_lungs, '_CombinedLungs', reporting);

            context_list_single_lungs = [PTKContext.LeftLung, PTKContext.RightLung];
            PTKSaveCoronalAnalysisResults.DrawGraphAndSave(dataset, table, y_label, context_list_single_lungs, '_Lungs', reporting);

            context_list_lobes = [PTKContext.LeftLowerLobe, PTKContext.LeftUpperLobe, PTKContext.RightLowerLobe, PTKContext.RightMiddleLobe, PTKContext.RightUpperLobe];
            PTKSaveCoronalAnalysisResults.DrawGraphAndSave(dataset, table, y_label, context_list_lobes, '_Lobes', reporting);

            results = [];
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded() && gui_app.IsCT();
        end
    end

    methods (Static, Access = private)
        function DrawGraphAndSave(dataset, table, y_label, context_list, file_suffix, reporting)
            figure_handle = PTKGraphMetricVsDistance(table, 'MeanDensityGml', 'StdDensityGml', context_list, [], y_label, reporting);
            dataset.SaveFigure(figure_handle, 'PTKSaveCoronalAnalysisResults', 'Coronal analysis', ['DensityVsCoronalDistance' file_suffix], 'Graph of density vs distance along the anterior-posterior axis');
            close(figure_handle);

            figure_handle = PTKGraphMetricVsDistance(table, 'EmphysemaPercentage', [], context_list, [], y_label, reporting);
            dataset.SaveFigure(figure_handle, 'PTKSaveCoronalAnalysisResults', 'Coronal analysis', ['EmphysemaVsCoronalDistance' file_suffix], 'Graph of emphysema vs distance along the anterior-posterior axis');
            close(figure_handle);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% Get the mask defining the points where the reconstructions should be
% computed.
function [ mask, lat, lon ] = metoffice_get_mask( dataset )

% filename from dataset
switch dataset
case {'hadsst2d1','hadsst2d2','mohsst5d1','mohsst5d2'}
    maskfile = 'mask';
case {'aari','aari8'}
    maskfile = 'mask_seaice';
otherwise
    error(sprintf('Unknown dataset %s',dataset));
end

% path from system
datapath = metoffice_get_path();

% load mask
if nargout == 3
    load( [ datapath, maskfile ], 'mask', 'lat', 'lon' )
elseif nargout == 1
    load( [ datapath, maskfile ], 'mask' )
end

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

function m = metoffice_get_mask_old(maskfile)

if nargin == 0
  warning('Give the maskfile');
  maskfile = 'mask';
end

os = system_dependent('getos');

if ~isempty(strfind(os,'Ubuntu'))
  % Finland
  maskfile = [ '/share/climate/data/UK_Met_Office/RecTest/DATASETS/' maskfile ];
elseif ~isempty(strfind(os,'Linux'))
  % ICS GridEngine
  maskfile = [ '/share/climate/data/UK_Met_Office/RecTest/DATASETS/' maskfile ];
else
  % MetOffice
  maskfile = [ '/net/home/h03/hadia/data/RecTest/' maskfile ];
end

S = load(maskfile);
m = S.mask(:)~=0;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
%% Copyright (C) 2015-2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym cbrt (@var{x})
%% Symbolic cbrt function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = cbrt (x)
%%   @result{} y = (sym)
%%       3 ___
%%       ╲╱ x
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = cbrt(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('cbrt', x);
end

%!error cbrt (sym(1), 2)
%!assert (isequaln (cbrt (sym(nan)), sym(nan)))

%!shared x, d
%! d = 2;
%! x = sym('2');

%!test
%! f1 = cbrt(x);
%! f2 = 1.2599210498948731647;
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = cbrt(A);
%! f2 = 1.2599210498948731647;
%! f2 = [f2 f2; f2 f2];
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! if (exist ('OCTAVE_VERSION', 'builtin'))
%!   A = cbrt (d);
%! else
%!   % Issue #742
%!   A = d^(1/3);
%! end
%! f = cbrt (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)

repo: octsympy


text:
% *************************************************************************
% Video Super-Resolution with Convolutional Neural Networks
%
% If you use this code, please cite the following publication:
%     A.Kappeler, S.Yoo, Q.Dai, A.K.Katsaggelos, "Video Super-Resolution
%     with Convolutional Neural Networks", to appear in IEEE Transactions
%     on Computational Imaging
%
% Installation Instructions:
% 1) install Caffe from "http://caffe.berkeleyvision.org/"
% 2) run the following command in Matlab:
%       cd external_functions/CLG-TV-matlab
%       mex applyBilateralFilterToDataTerms.cpp
%       cd ../..
% 3) specify the CAFFEPATH on line 53 in VSRnet_demo.m
% 4) set the experiment you want to execute to "true" (only one at the time)
%
%
% Version 1.0
%
% Created by:   Armin Kappeler
% Date:         02/19/2016
%
% http://ivpl.eecs.northwestern.edu/software
%
% *************************************************************************
%
% Because of the data size, we only provide a subset of our videos on the
% website. For additional testvideos and our training database, please
% contact us directly.
%
% *************************************************************************
% Copyright (C) 2016 Armin Kappeler
%
% 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.
%
% For a copy of the GNU General Public License,
% please see <http://www.gnu.org/licenses/>.
%
% *************************************************************************

clearvars
%% GENERAL PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CAFFEPATH = '/HOMES/baetz/develop/CAFFE/caffe/';   % path to caffe installation

LOW_MEMORY_MODE = 1;    % if Matlab crashes, try a higher number:
                        % 0 = high GPU memory usage -> fastest
                        % 1 = medium GPU memory usage -> fast
                        % 2 = low GPU memory usage -> slow

USE_GPU = false;         % set to false, if no GPU available -> slowest
GPU_ID = 0;             % GPU ID -> should normally be 0

%% RUN VSRnet on Video %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% specify TESTVIDEO_PATCH and UPSCALE_FACTOR and run
if true
    TESTVIDEO_PATH = 'data/lr_video/u2/foreman_LR.mat'; % path to input video frames
                                        % a variable "frames" with dimensions (width x height x nrFrames) is expected
    UPSCALE_FACTOR = 2;                 % upscale factor: 2,3,4 are available

    MOTIONCOMPENSATION = true;          % use Motion Compensation (much faster)
                                        % Note: MOTIONCOMPENSATION=false is only available for upscale factor 3
    ADAPTIVEMOTIONCOMPENSATION = false; % use Adaptive Motion Compensation
    TESTONLY1FRAME = 0;                 % 0: all frames will be tested, otherwise only one frame will be tested
    RESULTFILE_PATH = ['results/magi_quick_test'];% filename for results
    PREPROCESSED_INPUT = false;         % do not change
end

%% REPRODUCE RESULTS FROM PAPER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% these experiments reproduce the results from the paper
% the testvideos here are already preprocessed
% do not change the parameters below this point

% Myanmar upscale factor 2
if false
    TESTVIDEO_PATH = 'data/preprocessed/u2/matTestimages_Myanmar.mat';
    UPSCALE_FACTOR = 2;

    MOTIONCOMPENSATION = true;
    ADAPTIVEMOTIONCOMPENSATION = false;
    TESTONLY1FRAME = 0;
    RESULTFILE_PATH = ['results/dummy'];
    PREPROCESSED_INPUT = true;
end

% Myanmar upscale factor 3
if false
    TESTVIDEO_PATH = 'data/preprocessed/u3/matTestimages_Myanmar.mat';
    UPSCALE_FACTOR = 3;

    MOTIONCOMPENSATION = true;
    ADAPTIVEMOTIONCOMPENSATION = false;
    TESTONLY1FRAME = 0;
    RESULTFILE_PATH = ['results/dummy'];
    PREPROCESSED_INPUT = true;
end

% Myanmar upscale factor 4
if false
    TESTVIDEO_PATH = 'data/preprocessed/u4/matTestimages_Myanmar.mat';
    UPSCALE_FACTOR = 4;

    MOTIONCOMPENSATION = true;
    ADAPTIVEMOTIONCOMPENSATION = false;
    TESTONLY1FRAME = 0;
    RESULTFILE_PATH = ['results/dummy'];
    PREPROCESSED_INPUT = true;
end

% Myanmar upscale factor 3, no motion compensation
if false
    TESTVIDEO_PATH = 'data/preprocessed/u3/matTestimages_Myanmar_noMotionCompensation.mat';
    UPSCALE_FACTOR = 3;

    MOTIONCOMPENSATION = false;
    ADAPTIVEMOTIONCOMPENSATION = false;
    TESTONLY1FRAME = 0;
    RESULTFILE_PATH = ['results/dummy'];
    PREPROCESSED_INPUT = true;
end

% Foreman, upscale factor 3, normal motion compensation (MC)
if false
    TESTVIDEO_PATH = 'data/preprocessed/u3/matTestimages_foreman.mat';
    UPSCALE_FACTOR = 3;

    MOTIONCOMPENSATION = true;
    ADAPTIVEMOTIONCOMPENSATION = false;
    TESTONLY1FRAME = 7;
    RESULTFILE_PATH = ['results/dummy'];
    PREPROCESSED_INPUT = true;
end

% Foreman, upscale factor 3, adaptive motion compensation (AMC)
if false
    TESTVIDEO_PATH = 'data/preprocessed/u3/matTestimages_foreman_adaptiveMotionCompensation.mat';
    UPSCALE_FACTOR = 3;

    MOTIONCOMPENSATION = true;
    ADAPTIVEMOTIONCOMPENSATION = true;
    TESTONLY1FRAME = 7;
    RESULTFILE_PATH = ['results/dummy'];
    PREPROCESSED_INPUT = true;
end

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

%% start processing
addpath('functions','external_functions/CLG-TV-matlab',[CAFFEPATH '/matlab'])
run('main_VSRnet_test.m')
run('evaluate_result.m')

author: thomas-koehler

repo: SupER


text:
%FEACHECKMEX Check/test for valid mex file usage.
%
%   [ TF ] = FEACHECKMEX( MEX_FILE, ...ARGS ) Check/test for valid mex
%   file usage. MEX_FILE is an optional mex file to check/test against
%   (default chkmex/ftsparse), and ARGS optional input arguments in
%   test tex function call.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function sumll=LLtwoCIR(para,Y, tau, nrow, ncol)
% Initialize the parameters for CIR model
theta1=para(1);
kappa1=para(2);
sigma1=para(3);
lambda1=para(4);
theta2=para(5);
kappa2=para(6);
sigma2=para(7);
lambda2=para(8);
sigmai=para(9:end);

R=eye(ncol);
for i=1:ncol
    R(i,i)=sigmai(i)^2;
end
dt=1/12;

% System Matrices Initialization
C=[theta1*(1-exp(-kappa1*dt));theta2*(1-exp(-kappa2*dt))];   % eqn b.3
F=[exp(-kappa1*dt),0 ; 0,exp(-kappa2*dt)];                   % eqn b.3
A=zeros(ncol,1);
H=zeros(ncol,2);

% Create A and B
for i=1:ncol % System Matrices are made for each tau
    % Factor 1
    AffineG11=sqrt((kappa1+lambda1)^2+2*sigma1^2); % eqn a.8 & a.9 & a.10
    AffineG12=kappa1+lambda1+AffineG11;
    AffineG13=2*kappa1*theta1/sigma1^2;
    AffineG14=2*AffineG11+AffineG12*(exp(AffineG11*tau(i))-1);
    AffineA1=((2*AffineG11*exp(AffineG12*tau(i)/2))/AffineG14)^AffineG13;
    A1=-log(AffineA1)/tau(i);
    AffineB1=2*(exp(AffineG11*tau(i))-1)/AffineG14;
    B1=AffineB1/tau(i);
    % Factor 2
    AffineG21=sqrt((kappa2+lambda2)^2+2*sigma2^2); % % eqn a.8 & a.9 & a.10
    AffineG22=kappa2+lambda2+AffineG21;
    AffineG23=2*kappa2*theta2/sigma2^2;
    AffineG24=2*AffineG21+AffineG22*(exp(AffineG21*tau(i))-1);
    AffineA2=((2*AffineG21*exp(AffineG22*tau(i)/2))/AffineG24)^AffineG23;
    A2=-log(AffineA2)/tau(i);
    AffineB2=2*(exp(AffineG21*tau(i))-1)/AffineG24;
    B2=AffineB2/tau(i);

    A(i,1)=A1+A2;       % eqn b.2
    H(i,1)=B1;          % eqn b.2
    H(i,2)=B2;          % eqn b.2
end

%% Kalman filter
% Step 1
initx=[theta1;theta2];                      % eqn c.1
initV=[(sigma1^2*theta1)/(2*kappa1),0;...
        0,(sigma2^2*theta2)/(2*kappa2)];    % eqn c.2
% Starting values
AdjS=initx;
VarS=initV;
LL=zeros(nrow,1);          % log-likelihood vector initialization

for i=1:nrow
    PredS=C+F*AdjS;                 % eqn c.10
    Q=[(theta1*sigma1*sigma1*(1-exp(-kappa1*dt))^2/(2*kappa1)+sigma1*sigma1/kappa1*(exp(-kappa1*dt)-exp(-2*kappa1*dt)))*AdjS(1),0;...
        0,(theta2*sigma2*sigma2*(1-exp(-kappa2*dt))^2/(2*kappa2)+sigma2*sigma2/kappa2*(exp(-kappa2*dt)-exp(-2*kappa2*dt)))*AdjS(2)];    % eqn b.4
    VarS=F*VarS*F'+Q;               % eqn c.11
    % Step 2
    PredY=A+H*PredS;                    % eqn c.4
    VarY=H*VarS*H'+R;                   % eqn c.5
    % Step 3
    PredError=Y(i,:)'-PredY;            % eqn c.6
    InvVarY=inv(VarY);
    DetY=det(VarY);
    KalmanGain=VarS*H'*InvVarY;         % eqn c.8
    AdjS=PredS+KalmanGain*PredError;    % eqn c.7
    VarS=VarS*(1-KalmanGain*H);         % eqn c.9

    % Step 5
    LL(i)=-(ncol/2)*log(2*pi)-0.5*log(DetY)-0.5*PredError'*InvVarY*PredError; % eqn c.10
end
sumll=-sum(LL);
end


text: %MAX Dataset overload

author: marianux

repo: ecg-kit

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


text:
function [estimateclasstotal,model,D]=adaboost(mode,X,labels,model,itt,weak_type)
% This function AdaBoost, consist of two parts a simpel weak classifier and
% a boosting part:
% The weak classifier tries to find the best treshold in one of the data
% dimensions to sepparate the data into two classes -1 and 1
% The boosting part calls the clasifier iteratively, after every classification
% step it changes the weights of miss-classified examples. This creates a
% cascade of "weak classifiers" which behaves like a "strong classifier"
%
%
%   input -----------------------------------------------------------------
%
%       o mode   : string,  'train' or 'apply'.
%
%       o X      : (N x D), dataset of N samples of dimension D.
%
%       o labels : (N x 1), class labels {-1,+1}
%
%       o model  : struct of model learned, set ot [] when training.
%
%       o itt    :  number of iterations to perform which is equivalent to
%                   the number of resulting weak classifiers.
%
%       o weak_type : string, weak classifier type:
%
%               - decision_stump : splits space along the dimension
%
%   output ----------------------------------------------------------------
%
%  inputs/outputs:
%    datafeatures : An Array with size number_samples x number_features
%    dataclass : An array with the class off all examples, the class
%                 can be -1 or 1
%    itt : The number of training itterations
%    model : A struct with the cascade of weak-classifiers
%    estimateclass : The by the adaboost model classified data
%

if ~exist('weak_type','var'), weak_type = 'decision_stump'; end

switch(mode)
    case 'train'
        % Train the adaboost model

        % Set the data class labels
        labels    = labels(:);
        model     = struct;

        % Weight of training samples, first every sample is even important
        % (same weight)
        D=ones(length(labels),1)/length(labels);

        % This variable will contain the results of the single weak
        % classifiers weight by their alpha
        predict_label   =   zeros(size(labels));

        % Calculate max min of the data
        boundary=[min(X,[],1) max(X,[],1)];
        % Do all model training itterations
        for t=1:itt
            % Find the best treshold to separate the data in two classes

            if strcmpi(weak_type,'decision_stump')

                [estimateclass,err,h] = WeightedThresholdClassifier(X,labels,D);

            elseif strcmpi(weak_type,'linear')

                [B_hat,mapping] = train_linear_classifier(X,labels,D);
                [estimateclass] = predictor_linear_classifier(X,B_hat,mapping);
                err             = sum(D(:) .* (estimateclass ~= labels));
            end

            % Weak classifier influence on total result is based on the current
            % classification error
            alpha=1/2 * log((1-err)/max(err,eps));

            % Store the model parameters
            model(t).alpha      = alpha;

            if strcmpi(weak_type,'decision_stump')
                model(t).dimension  = h.dimension;
                model(t).threshold  = h.threshold;
                model(t).direction  = h.direction;
                model(t).boundary   = boundary;
            elseif strcmpi(weak_type,'linear')
                model(t).B_hat   = B_hat;
                model(t).mapping = mapping;
            end

            % We update D so that wrongly classified samples will have more weight
            D = D.* exp(-model(t).alpha.*labels.*estimateclass);
            D = D./sum(D);

            % Calculate the current error of the cascade of weak
            % classifiers
            predict_label       =   predict_label + estimateclass * model(t).alpha;
            estimateclasstotal  =   sign(predict_label);

            model(t).error      =   sum(estimateclasstotal~=labels)/length(labels);

            if(model(t).error==0), break; end

        end

    case 'apply'
        % Apply Model on the test data

        if strcmpi(weak_type,'decision_stump')

            % Limit datafeatures to orgininal boundaries
            if(length(model)>1);
                minb=model(1).boundary(1:end/2);
                maxb=model(1).boundary(end/2+1:end);
                X=bsxfun(@min,X,maxb);
                X=bsxfun(@max,X,minb);
            end

            % Add all results of the single weak classifiers weighted by their alpha
            predict_label=zeros(size(X,1),1);
            for t=1:length(model);
                predict_label=predict_label+model(t).alpha*ApplyClassTreshold(model(t), X);
            end
            % If the total sum of all weak classifiers
            % is less than zero it is probablly class -1 otherwise class 1;
        elseif strcmpi(weak_type,'linear')

            predict_label = zeros(size(X,1),1);

            for t=1:length(model);
                estimateclass = predictor_linear_classifier(X,model(t).B_hat,model(t).mapping);
                predict_label = predict_label + model(t).alpha * estimateclass;
            end

        end

        estimateclasstotal = sign(predict_label);

    otherwise
        error('adaboost:inputs','unknown mode');
end

author: epfl-lasa


text:
function [ ap, rcond, z, info ] = sppco ( ap, n )

%*****************************************************************************80
%
%% SPPCO factors a real symmetric positive definite matrix in packed form.
%
%  Discussion:
%
%    SPPCO also estimates the condition of the matrix.
%
%    If RCOND is not needed, SPPFA is slightly faster.
%
%    To solve A*X = B, follow SPPCO by SPPSL.
%
%    To compute inverse(A)*C, follow SPPCO by SPPSL.
%
%    To compute determinant(A), follow SPPCO by SPPDI.
%
%    To compute inverse(A), follow SPPCO by SPPDI.
%
%  Packed storage:
%
%    The following program segment will pack the upper triangle of
%    a symmetric matrix.
%
%      k = 0
%      do j = 1, n
%        do i = 1, j
%          k = k + 1
%          ap(k) = a(i,j)
%        end
%      end
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real AP(N*(N+1)/2), the packed form of a symmetric matrix A.
%    The columns of the upper triangle are stored sequentially in a
%    one-dimensional array.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real AP(N*(N+1)/2), an upper riangular matrix R, stored
%    in packed form, so that A = R'*R.  If INFO /= 0, the factorization
%    is not complete.
%
%    Output, real RCOND, an estimate of the reciprocal condition
%    of A.  For the system A*X = B, relative perturbations in A and B of size
%    EPSILON may cause relative perturbations in X of size EPSILON/RCOND.
%    If RCOND is so small that the logical expression
%      1.0 + RCOND == 1.0D+00
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate underflows.
%
%    Output, real Z(N), a work vector whose contents are usually
%    unimportant.  If A is singular to working precision, then Z is an
%    approximate null vector in the sense that
%      norm(A*Z) = RCOND * norm(A) * norm(Z).
%    If INFO /= 0, Z is unchanged.
%
%    Output, integer INFO, error flag.
%    0, for normal return.
%    K, signals an error condition.  The leading minor of order K is
%    not positive definite.
%

%
%  Find the norm of A.
%
  j1 = 1;
  for j = 1 : n
    z(j) = sasum ( j, ap(j1:j1+j-1), 1 );
    ij = j1;
    j1 = j1 + j;
    for i = 1 : j-1
      z(i) = z(i) + abs ( ap(ij) );
      ij = ij + 1;
    end
  end

  anorm = max ( z(1:n) );
%
%  Factor.
%
  [ ap, info ] = sppfa ( ap, n );

  if ( info ~= 0 )
    return
  end
%
%  RCOND = 1/(norm(A)*(estimate of norm(inverse(A)))).
%
%  Estimate = norm(Z)/norm(Y) where A * Z = Y and A * Y = E.
%
%  The components of E are chosen to cause maximum local
%  growth in the elements of W where R'*W = E.
%
%  The vectors are frequently rescaled to avoid overflow.
%
%  Solve R' * W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  kk = 0;

  for k = 1 : n

    kk = kk + k;

    if ( z(k) ~= 0.0 )
      ek = -abs ( ek ) * r4_sign ( z(k) );
    end

    if ( ap(kk) < abs ( ek - z(k) ) )
      s = ap(kk) / abs ( ek - z(k) );
      z(1:n) = s * z(1:n);
      ek = s * ek;
    end

    wk = ek - z(k);
    wkm = -ek - z(k);
    s = abs ( wk );
    sm = abs ( wkm );
    wk = wk / ap(kk);
    wkm = wkm / ap(kk);
    kj = kk + k;

    if ( k + 1 <= n )

      for j = k + 1 : n
        sm = sm + abs ( z(j) + wkm * ap(kj) );
        z(j) = z(j) + wk * ap(kj);
        s = s + abs ( z(j) );
        kj = kj + j;
      end

      if ( s < sm )

        t = wkm - wk;
        wk = wkm;
        kj = kk + k;

        for j = k+1 : n
          z(j) = z(j) + t * ap(kj);
          kj = kj + j;
        end

      end

    end

    z(k) = wk;

  end

  z(1:n) = z(1:n) / sasum ( n, z(1:n), 1 );
%
%  Solve R * Y = W.
%
  for k = n : -1 : 1

    if ( ap(kk) < abs ( z(k) ) );
      s = ap(kk) / abs ( z(k) );
      z(1:n) = s * z(1:n);
    end

    z(k) = z(k) / ap(kk);
    kk = kk - k;
    t = -z(k);
    z(1:k-1) = saxpy ( k-1, t, ap(kk+1:kk+k-1), 1, z(1:k-1), 1 );

  end

  z(1:n) = z(1:n) / sasum ( n, z(1:n), 1 );

  ynorm = 1.0;
%
%  Solve R' * V = Y.
%
  for k = 1 : n

    z(k) = z(k) - sdot ( k-1, ap(kk+1:kk+k-1), 1, z(1:k-1), 1 );
    kk = kk + k;

    if ( ap(kk) < abs ( z(k) ) )
      s = ap(kk) / abs ( z(k) );
      z(1:n) = s * z(1:n);
      ynorm = s * ynorm;
    end

    z(k) = z(k) / ap(kk);

  end

  s = 1.0 / sasum ( n, z(1:n), 1 );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;
%
%  Solve R * Z = V.
%
  for k = n : -1 : 1

    if ( ap(kk) < abs ( z(k) ) )
      s = ap(kk) / abs ( z(k) );
      z(1:n) = s * z(1:n);
      ynorm = s * ynorm;
    end

    z(k) = z(k) / ap(kk);
    kk = kk - k;
    t = -z(k);
    z(1:k-1) = saxpy ( k-1, t, ap(kk+1:kk+k-1), 1, z(1:k-1), 1 );

  end
%
%  Make ZNORM = 1.0.
%
  s = 1.0 / sasum ( n, z(1:n), 1 );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;

  if ( anorm ~= 0.0 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0;
  end

  return
end


text:
function output = callqpas(interfacedata)

% Retrieve needed data
options = interfacedata.options;
model = yalmip2quadprog(interfacedata);

if size(model.Aeq,1)==0
    model.Aeq = [];
    model.beq = [];
end

if options.savedebug
    save debugfile model
end

if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end
solvertime = tic;
[x,flag,lm] = qpas(model.Q, model.c, model.A, model.b, model.Aeq, model.beq, model.lb, model.ub, options.verbose);
solvertime = toc(solvertime);

% Internal format for duals
if ~isempty(lm)
    D_struc = [lm.equality;lm.inequality];
else
    D_struc = [];
end

if isempty(flag)
    problem = 9;
else
    % Check, currently not exhaustive...
    switch flag
        case 0
            problem = 0;
        case 2
            problem = 1;
        otherwise
            problem = 1;
    end
end

infostr = yalmiperror(problem,interfacedata.solver.tag);

% Save all data sent to solver?
if options.savesolverinput
    solverinput = model;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.x = x;
    solveroutput.flag = flag;
    solveroutput.lm = lm;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x(:),D_struc,[],problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [bestAlgo, inliers] = prtUtilRansac(data,algo,evalFn,varargin)
%[bestAlgo, inliers] = prtUtilRansac(data,algo,evalFn)
%   Perform RANSAC on prtDataSet data.
%
%[bestAlgo, inliers] = prtUtilRansac(data,algo,evalFn,param1,value1,...)
%   Enables inputs of parameter/value pairs as described below:
%
%   Parameters:
%       nIterations (100)
%       nBootstrapSamples (10) - the number of samples to use at each
%         iteration to fit a model.
%       fitErrorThreshold (0.2) - absolute error between truth and guess
%         should be less than fitErrorThreshold to be included inliers.
%
%   Example usage:
%
%       ds = prtDataGenRansac(100,0.2);
%       [~,inliers] = prtUtilRansac(ds,prtRegressLslr,@(a,b)abs(a-b));
%       plot(ds.X(inliers),ds.Y(inliers),'ro',...
%            ds.X(setdiff(1:100,inliers)),ds.Y(setdiff(1:100,inliers)),'bs')

% Copyright (c) 2014 Patrick Wang, Peter Torrione, Kenneth Morton
%
% Permission is hereby granted, free of charge, to any person obtaining a
% copy of this software and associated documentation files (the
% "Software"), to deal in the Software without restriction, including
% without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Software, and to permit
% persons to whom the Software is furnished to do so, subject to the
% following conditions:
%
% The above copyright notice and this permission notice shall be included
% in all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
% NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
% DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
% OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
% USE OR OTHER DEALINGS IN THE SOFTWARE.

%% parse input options
p = inputParser;

p.addParamValue('nIterations',100);
p.addParamValue('nBootstrapSamples',10);
p.addParamValue('fitErrorThreshold',0.2);

p.parse(varargin{:});
inputStructure = p.Results;

%% iterate
bestNFit = 0;
bestAlgo = [];
for iter = 1:inputStructure.nIterations
	% bootstrap
	indices = randperm(data.nObservations);
	indices = indices(1:inputStructure.nBootstrapSamples);
	dataIn = data.retainObservations(indices);
	dataOut = data.removeObservations(indices);

	% train on bootstrapped samples
	algoTrained = algo.train(dataIn);

	% test on left-out samples
	guessOut = algoTrained.run(dataOut);

	% find how many left-out samples fit model well
	error = evalFn(guessOut.X,dataOut.Y); % dataOut.Y?
	nFit = sum(error<inputStructure.fitErrorThreshold);

	% keep the best model
	if nFit > bestNFit
		bestNFit = nFit;
		bestAlgo = algoTrained;
	end
end

%% find all inliers and retrain model
guess = bestAlgo.run(data);
error = evalFn(guess.X,data.Y); % data.Y?
inliers = find(error<inputStructure.fitErrorThreshold);
bestAlgo = algo.train(data.retainObservations(inliers));

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
% Nonlinear model of two CSTRs in series from
%
% Henson, M.A. and Seborg, D.E., Feedback Linearizing Control, Chap. 4 of Nonlinear Process
%   Control, Edited by Hensen, M.A. and Seborg, D.E., Prentice Hall (1997)
%
%  t -- time (not used)
%  y -- state value vector
%  xdot -- set to the vector of state derivatives

function xdot = cstr5(t,y)

global u

% Input (1):
% Coolant Flowrate (L/min)
qc = u;

% States (4):
% Concentration of A in Reactor #1 (mol/L)
Ca1 = y(1);
% Temperature of Reactor #1 (K)
T1 = y(2);
% Concentration of A in Reactor #2 (mol/L)
Ca2 = y(3);
% Temperature of Reactor #2 (K)
T2 = y(4);

% Parameters
% Flowrate (L/min)
q = 100;
% Feed Concentration of A (mol/L)
Caf = 1.0;
% Feed Temperature (K)
Tf	= 350.0 ;
% Coolant Temperature (K)
Tcf = 350.0;
% Volume of Reactor #1 (L)
V1 = 100;
% Volume of Reactor #2 (L)
V2 = 100;
% UA1 or UA2 - Overall Heat Transfer Coefficient (J/min-K)
UA1 = 1.67e5;
UA2 = 1.67e5;
% Pre-exponential Factor for A->B Arrhenius Equation
k0 = 7.2e10;
% EoverR - E/R (K) - Activation Energy (J/mol) / Gas Constant (J/mol-K)
EoverR = 1e4;
% Heat of Reaction - Actually (-dH) for an exothermic reaction (J/mol)
dH = 4.78e4;
% Density of Fluid (g/L)
rho = 1000;
% Density of Coolant Fluid (g/L)
rhoc = 1000;
% Heat Capacity of Fluid (J/g-K)
Cp = 0.239;
% Heat Capacity of Coolant Fluid (J/g-K)
Cpc = 0.239;

%	Dynamic Balances

%	dCa1/dt
xdot(1,1) = q*(Caf-Ca1)/V1 - k0*Ca1*exp(-EoverR/T1);
%	dT1/dt
xdot(2,1) = q*(Tf-T1)/V1 + (dH*k0/rho/Cp)*Ca1*exp(-EoverR/T1) + ...
   rhoc*Cpc/rho/Cp/V1 * qc * (1-exp(-UA1/qc/rhoc/Cpc)) * (Tcf-T1);
%	dCa2/dt
xdot(3,1) = q*(Ca1-Ca2)/V2 - k0*Ca2*exp(-EoverR/T2);
%	dT2/dt
xdot(4,1)  = q*(T1-T2)/V2 + (dH*k0/rho/Cp)*Ca2*exp(-EoverR/T2) + ...
   rhoc*Cpc/rho/Cp/V2 * qc * (1-exp(-UA2/qc/rhoc/Cpc)) * ...
   (T1 - T2 + exp(-UA1/qc/rhoc/Cpc)*(Tcf-T1));

author: Sable

repo: mcbench-benchmarks


text:
% Evaluate fitting performance on different datasets.
%
% Author: Luca Vezzaro (elvezzaro@gmail.com)

clear all;
close all;

addpath('..');

datasets = { 'child', 'Cootes', 'mother', 'movie1', 'movie2'};
convergence_thresholds = {2, 7, 3, 5.5, 5.5};
test = {[1 2 3 4 5 10 13], [8 10 12 22 26 ], [2 3 14], [7 12 5 21], [1 3 22 31 ]};

% Selected dataset
sel = 3;
blur_training_set = 1;

base_path = sprintf('performance/%s', datasets{sel});

% Available landmark files
landmark_files = dir(sprintf('%s/training/*.mat', base_path));
ns = numel(landmark_files);

for j=1:ns
	% Prepare appearance data
	app = double(imread(sprintf('%s/training/%s', base_path, landmark_files(j).name(1:end-4)))) / 255;

	if blur_training_set
		app = convn(convn(app, exp(-(-5:5).^2)./sum(exp(-(-5:5).^2)), 'same'), exp(-(-5:5).^2)./sum(exp(-(-5:5).^2))', 'same');
	end

	appearances(:,:,:,j) = app;
end

for j=1:ns
	% Load associated landmarks
	load(sprintf('%s/training/%s', base_path, landmark_files(j).name));
	shapes(:,:,j) = xy2ij(annotations, size(app,1));
end

load(sprintf('%s/triangulation.mat', base_path));

np = size(shapes, 1);

% Sum of fitting times
tempi = 0;
% Number of fitting times considered
ntempi = 0;

% Number of times the fitting is evaluated
ntries = zeros(10,1);
% Number of times the fitting successfully produced a result (ie. the warp didn't diverge)
nfitted = zeros(10,1);
% Number of times the pt.pt error was below the convergence threshold
nconverged = zeros(10,1);
% Sum of ptpt errors for a given value of 'k'. The number of values considered here is given by 'nfitted'.
ptpt = zeros(10,1);
% Minimum pt.pt. error for a given value of 'k'
ptpt_min = ones(10,1) * inf;
% Maximum pt.pt. error for a given value of 'k'
ptpt_max = zeros(10,1);

% Carry out leave-one-out testing
for jj=1:numel(test{sel})
	j = test{sel}(jj);

	fprintf('\n\n\n****************************************************************************\n');
	fprintf('Carrying out %d/%d leave-one-out try on dataset: %s\n', jj, numel(test{sel}), datasets{sel});
	fprintf('****************************************************************************\n');

	one_out = [1:j-1,j+1:ns];
	AAM = build_model_2d(shapes(:,:,one_out), appearances(:,:,:,one_out), 'triangulation', triangulation);

	for k=1:10
		fprintf('%d/10\n', k);

		d1 = [rand() rand()] - [0.5 0.5];
		d1 = 3*d1 / norm(d1);

		d2 = [rand() rand()] - [0.5 0.5];
		d2 = 3*d2 / norm(d2);

		% dx dy
		%displacements = k*0.6*[0 0; 0 5; 5 0; 0 -5; -5 0; -3.5355 3.5355; 3.5355 3.5355; 3.5355 -3.5355; -3.5355 -3.5355];
		displacements = k*0.6*[0 5; 5 0; 0 -5; -5 0];
		%displacements = k*[d1; d2];
		rotations = k*[0 pi/80 -pi/80 ];
		scales = [1 1 1 ] + k*[0 0.02 -0.02];

		cog = repmat(mean(shapes(:,:,j), 1), [np 1]);

		for ti=1:size(displacements, 1)
			for rot=1:numel(rotations)
				for sc=1:numel(scales)
					ntries(k) = ntries(k) + 1;
					theta = rotations(rot);

					A = [cos(theta) sin(theta); -sin(theta) cos(theta)] * scales(sc);
					d = repmat(displacements(ti,:), [np, 1]);

					init_shape = (shapes(:,:,j) - cog) * A + d + cog;
					%init_shape = shapes(:,:,j) + d;

					%ptpt_init = mean(sqrt(sum((init_shape - shapes(:,:,j)).^2, 2)));
					%rms_init = sqrt(sum(sum((init_shape - shapes(:,:,j)).^2)) / (2*np));

					diverged = 0;
					try
						tic;
						[ fitted_shape ] = fit_2d(AAM, init_shape, appearances(:,:,:,j), 20);
						dt = toc;
						nfitted(k) = nfitted(k) + 1;
					catch
						dt = toc;
						diverged = 1;
					end

					if diverged == 0
						ptpt_err = mean(sqrt(sum((fitted_shape - shapes(:,:,j)).^2, 2)))
						%rms_err = sqrt(sum(sum((fitted_shape - shapes(:,:,j)).^2)) / (2*np));

						ptpt(k) =  ptpt(k) + ptpt_err;

						if ptpt_err < ptpt_min(k)
							ptpt_min(k) = ptpt_err;
						end

						if ptpt_err > ptpt_max(k)
							ptpt_max(k) = ptpt_err;
						end

						tempi = tempi + dt;
						ntempi = ntempi + 1;

						% Check if error is acceptable
						if ptpt_err < convergence_thresholds{sel}
							nconverged(k) = nconverged(k) + 1;
						end
					end

					%hold off;
					%imshow(appearances(:,:,:,j));
					%hold on;
					%triplot(triangulation, init_shape(:,2), init_shape(:,1), 'b', 'LineWidth', 2);
					%if diverged == 0
					%	triplot(triangulation, fitted_shape(:,2), fitted_shape(:,1), 'w', 'LineWidth', 2);
					%end
					%pause;
				end
			end
		end
	end
end

% Convergence rate as a function of 'k'
results.conv_rate = nconverged ./ ntries;
% Mean fitting time
results.time_mean = tempi / ntempi;
% Mean pt.pt error as a function of 'k'
results.ptpt_mean = ptpt ./ nfitted;

results.ptpt_min = ptpt_min;
results.ptpt_max = ptpt_max;

results.nfitted = nfitted;
results.ntries = ntries;
results.nconverged = nconverged;

disp 'Plotting convergence ratio as a function of k...'
figure;
plot(1:numel(nconverged), results.conv_rate);

result_path = sprintf('%s/results.mat', base_path);
save(result_path, 'results');
fprintf('\For the full results of the performance evaluation see: ./%s. See that file for more details.\n', result_path);

clear all;


text:
function imdb = setup_voc07_ssw(varargin)
% setup_voc07_ssw  Initialize PASCAL VOC2007 data with selective
% search windows

% Warning! boxes are in the format of ([y1 x1 y2 x2])

opts.dataDir = fullfile('data') ;
opts.proposalDir = fullfile(opts.dataDir,'SSW');
opts.loadTest = 1;
opts = vl_argparse(opts, varargin) ;

% -------------------------------------------------------------------------
%                                                 Load selective search win
% -------------------------------------------------------------------------
%% get selective search windows
files = {'SelectiveSearchVOC2007trainval.mat', ...
  'SelectiveSearchVOC2007test.mat'} ;

if ~exist(opts.proposalDir, 'dir')
  mkdir(opts.proposalDir) ;
end

for i=1:numel(files)
  if ~exist(fullfile(opts.proposalDir, files{i}), 'file')
    url = sprintf('http://koen.me/research/downloads/%s',files{i}) ;
    fprintf('downloading %s\n', url) ;
    urlwrite(url,[opts.proposalDir filesep files{i}]);
  end
end

if ~isempty(opts.proposalDir)
  t1 = load([opts.proposalDir,filesep,files{1}]);
  if opts.loadTest
    t2 = load([opts.proposalDir,filesep,files{2}]);
    ssw.id = [str2double(t1.images);str2double(t2.images)]';
    ssw.boxes = cat(2,t1.boxes,t2.boxes);
  else
    ssw.id = str2double(t1.images)';
    ssw.boxes = t1.boxes;
  end

  [~,si] = sort(ssw.id);
  ssw.id = ssw.id(si);
  ssw.boxes = ssw.boxes(si);
end

% -------------------------------------------------------------------------
%                                                  Load categories metadata
% -------------------------------------------------------------------------
cats = {'aeroplane','bicycle','bird','boat','bottle','bus','car',...
  'cat','chair','cow','diningtable','dog','horse','motorbike','person',...
  'pottedplant','sheep','sofa','train','tvmonitor'};

if ~exist(opts.dataDir,'dir')
  error('wrong data folder!');
end

if ~exist(opts.dataDir,'dir')
  error('wrong data folder!');
end

% Download VOC Devkit and data
if ~exist(fullfile(opts.dataDir,'VOCdevkit'),'dir')
  files = {'VOCtest_06-Nov-2007.tar',...
           'VOCtrainval_06-Nov-2007.tar',...
           'VOCdevkit_08-Jun-2007.tar'} ;
  for i=1:numel(files)
    if ~exist(fullfile(opts.dataDir, files{i}), 'file')
      outPath = fullfile(opts.dataDir,files{i}) ;
      url = sprintf('http://host.robots.ox.ac.uk/pascal/VOC/voc2007/%s',files{i}) ;
      fprintf('Downloading %s to %s\n', url, outPath) ;
      urlwrite(url,outPath) ;
      untar(outPath,opts.dataDir);
    end
  end
end
addpath(fullfile(opts.dataDir, 'VOCdevkit', 'VOCcode'));

traindata = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main','train.txt'));
valdata = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main','val.txt'));
testdata = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main','test.txt'));

assert(numel(traindata)==2501);
assert(numel(valdata)==2510);
assert(numel(testdata)==4952);

imdb.classes.name = cats ;
imdb.classes.description = cats ;
imdb.imageDir = fullfile(opts.dataDir, fullfile('VOCdevkit','VOC2007','JPEGImages')) ;

% -------------------------------------------------------------------------
%                                                           Training images
% -------------------------------------------------------------------------%
names = cell(1,numel(traindata));
labels = zeros(numel(traindata),numel(cats));

% load image names
for t=1:numel(traindata)
  names{t} = sprintf('%06d.jpg',traindata(t));
%   data{t} = imread(sprintf('%s/%s',imdb.imageDir,names{t}));
end

% load binary labels
for c=1:numel(cats)
  t = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main',[cats{c},'_train.txt']));
  labels(:,c) = t(:,2);
end

imdb.images.id = traindata';
imdb.images.name = names ;
imdb.images.set = ones(1, numel(names)) ;
imdb.images.label = labels' ;
% imdb.images.data = data;

% -------------------------------------------------------------------------
%                                                         Validation images
% -------------------------------------------------------------------------

names = cell(1,numel(valdata));
labels = zeros(numel(valdata),numel(cats));
% data = cell(1,numel(valdata));

% load image names
for t=1:numel(valdata)
  names{t} = sprintf('%06d.jpg',valdata(t));
%   data{t} = imread(sprintf('%s/%s',imdb.imageDir,names{t}));
end

% load binary labels
for c=1:numel(cats)
  t = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main',[cats{c},'_val.txt']));
  labels(:,c) = t(:,2);
end

imdb.images.id = horzcat(imdb.images.id, valdata') ;
imdb.images.name = horzcat(imdb.images.name, names) ;
imdb.images.set = horzcat(imdb.images.set, 2*ones(1,numel(names))) ;
imdb.images.label = horzcat(imdb.images.label, labels') ;
% imdb.images.data = horzcat(imdb.images.data, data) ;

% % -------------------------------------------------------------------------
% %                                                               Test images
% % -------------------------------------------------------------------------
%
%
if opts.loadTest
  names = cell(1,numel(testdata));
  labels = zeros(numel(testdata),numel(cats));
  % data = cell(1,numel(testdata));

  % load image names
  for t=1:numel(testdata)
    names{t} = sprintf('%06d.jpg',testdata(t));
  %   data{t} = imread(sprintf('%s/%s',imdb.imageDir,names{t}));
  end

  % load binary labels
  for c=1:numel(cats)
    t = importdata(fullfile(opts.dataDir,'VOCdevkit','VOC2007','ImageSets','Main',[cats{c},'_test.txt']));
    labels(:,c) = t(:,2);
  end

  imdb.images.id = horzcat(imdb.images.id, testdata') ;
  imdb.images.name = horzcat(imdb.images.name, names) ;
  imdb.images.set = horzcat(imdb.images.set, 3 * ones(1,numel(names))) ;
  imdb.images.label = horzcat(imdb.images.label, labels') ;
  % imdb.images.data = horzcat(imdb.images.data, data) ;
end
% -------------------------------------------------------------------------
%                                                            Postprocessing
% -------------------------------------------------------------------------
[~,sorti] = sort(imdb.images.id);

imdb.images.id = imdb.images.id(sorti);
imdb.images.name = imdb.images.name(sorti) ;
imdb.images.set = imdb.images.set(sorti) ;
imdb.images.label = single(imdb.images.label(:,sorti)) ;
imdb.images.size = zeros(numel(imdb.images.name),2);

if ~isempty(opts.proposalDir)
  imdb.images.boxes = ssw.boxes;
  assert(all(ssw.id==imdb.images.id));
end

% this is zero as scores of selective search windows are not much
% informative
if ~isempty(opts.proposalDir)
imdb.images.boxScores = cell(size(imdb.images.boxes));
for i=1:numel(imdb.images.boxes)
  imdb.images.boxes{i} = int16(imdb.images.boxes{i});
  imdb.images.boxScores{i} = zeros(size(imdb.images.boxes{i},1),1,'single');
  imf = imfinfo(fullfile(imdb.imageDir,imdb.images.name{i}));
  imdb.images.size(i,:) = [imf.Height,imf.Width];
end
end
end

author: hbilen

repo: WSDDN

save_path: github-repos/MATLAB/hbilen-WSDDN


text:
function [bvals, bvec, f_bvec, f_bvals] = scd_scheme2bvecsbvals(scheme, acq_basename)
% scd_scheme2bvecs(scheme, acq_basename)
% EXAMPLE:  scheme = scd_schemefile_read(ls('*.scheme'));
%           [bvals, bvec] = scd_scheme2bvecsbvals(scheme);

bvec = scheme(:,1:3);
gyro = 42.57; % kHz/mT
bvals = (2*pi*gyro*scheme(:,4).*scheme(:,6).*10^(3)).^2.*(scheme(:,5)-scheme(:,6)/3)*10^(-3);

if nargin>1
f_bvec=[acq_basename '.bvec.txt'];
f_bvals=[acq_basename '.bvals.txt'];
fid_bvec_tot = fopen(f_bvec,'w');
fid_bvals_tot = fopen(f_bvals,'w');

for i=1:size(scheme,1)
        % write bvecs
        fprintf(fid_bvec_tot, '%f %f %f\n',bvec(i,:));
        % write bvals
        fprintf(fid_bvals_tot, '%f\n',bvals(i));

end

fclose all;

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%GOBJ_REVOLVE Revolve face/geometry object.
%
%   [ GOBJ, STAT ] = GOBJ_REVOLVE( GOBJ, TH, P, V, PIT, FACE_WPL )
%   Revolve geometry object or face from GOBJ a an angle TH (degrees)
%   around axis defined by point P and direction vector V. PIT is
%   optionally the pitch (distance per revolution) to offset the
%   revolution in the axial direction (to form a coil/helix).
%   FACE_WPL can in 3D indicate a face to extrude or in 2D is a
%   workplane struct defined by (p, n, t).

% Copyright 2013-2022 Precise Simulation, Ltd.


text:
function h = cutfunc(x,a,b)
c1 = (x <= a);
c2 = (x>a & x<b);
c3 = (x >= b);
h = 0*c1+ (x-a)/(b-a).*c2 + 1*c3;

% for i=2:nx
%
%     h(1,i)= (x(1,i)-1)/(-1);
%
%     h(1,1)=1;
% end

repo: Aero-matlab