text:
function [newDietModel,pointsModel,roiFlux,pointsModelSln,menuChanges] = nutritionAlgorithmWBM(model,obj,objMinMax,rois,roisMinMax,options)
% Identifies the minimal changes to a diet necessary to get a desired
% change in one or more reactions of interest. One may enter a metabolite
% of the pointsModel instead of a reaction and the algorithm will optimize the
% diet with a sink or demand reaction for the corresponding metabolite of
% interest.
%
% USAGE:
%
%    [newDietModel,pointsModel,slnMin,slnMax,pointsModelSln,itemsRemoved,itemsAdded] = nutritionAlgorithmWBM(pointsModel,obj,objMinMax,rois,roisMinMax,options)
%
%     Example: [newDietModel,pointsModel,roiFlux,pointsModelSln,itemsRemoved,itemsAdded] = nutritionAlgorithmWBM(WBmodel,'Whole_body_objective_rxn','max',{},{})
%
% INPUTS:
%    pointsModel:          COBRA pointsModel structure with the fields:
%                      * .S
%                      * .b
%                      * .ub
%                      * .ub
%                      * .mets  (required if pointsModel.SIntRxnBool absent)
%                      * .rxns  (required if pointsModel.SIntRxnBool absent)
%
%   obj:           organism's objective function
%
%   objMinMax:     minimize ('min') or maximize ('max') objective function
%
%   rois:          cell array of all reactions of interest
%
%   roisMinMax:    cell array of 'min'/'max' entries for rois
%
% OPTIONAL INPUTS:
%   options:  Structure containing the optional specifications:
%
%       * .foodOrMets:  dictates if the algorithm adds individual
%       metabolites to the diet or food items. Default is food items.
%       "Food Cat" adjust algorithm to identify categories of food rather
%       than specific items."AllMets" allows any dietary metabolite into
%       the solution and "FoodMets" only allows metabolites that are in
%       the fdTable spreadsheet into the solution.
%       Possible inputs are: "Food Items", "Food Cat", "AllMets", "FoodMets".
%
%       * .roiWeights:   a vector of weights for each reaction of interest
%       default is equal to 1
%
%       * .weightedFoodItems: A cell vector that specifies any food items
%       or metabolites that should be weighted and the corresponding weight.
%
%       * .initObjSln: provide an initial solution for the objective
%       function. Output from optimizeWBmodel.
%
%       * .caloricRange: 1x2 vector defining boundries for diet calories
%
%       * .slnType: Specify if solution should be 'Detailed' or 'Quick'.
%                   Default setting is 'Detailed'
%
%       * .roiBound: 'Unbounded' or 'Bounded'. Default is 'Bounded'.
%
%       * .foodAddedLimit: Specify a limit for the units of food that can
%                          be added to the diet
%
%       * .foodRemovedLimit: Specify a limit for the units of food that can
%                          be removed from the diet
%
%       * .freeMets: Specifies any metabolites that should be freely
%       available to the model.
%
%       * .calorieWeight: set to 'True' to weight by caloric content rather
%       than servings. Default is 'False'
%
%       * .graphicalAnalysis: set to 'True' include graphical analysis and
%       'False' to not include. Default is 'False' if .slnType is set to
%       'Quick but is 'True; if .slnType is 'Detailed'.
%
% OUTPUT:
%    solution:       Structure containing the following fields:
%
% relaxedModel       pointsModel structure that admits a flux balance solution
%
% .. Authors: - Bronson R. Weston   2021-2022

disp('_____________________________________________________')
%is foodOrMets variable established in the options struct?
if exist('options','var')
    if isfield(options,'foodOrMets')
        foodOrMets=options.foodOrMets;
    else
        foodOrMets='Food Items';
    end
end

%Identify which tables should be loaded
if strcmp(foodOrMets,'Food Items')
    load('fdTable.mat')
elseif strcmp(foodOrMets,'Food Cat')
    try
        load('FoodCategories/fdCategoriesTable.mat')
    catch
        load('fdCategoriesTable.mat')
    end
    fdTable=fdCategoriesTable;
elseif ~strcmp(foodOrMets,'FoodMets') && ~strcmp(foodOrMets,'AllMets')
    error('foodOrMets invalid. Possible inputs are: "Food Items", "Food Cat", "AllMets", "FoodMets".')
else
    options.graphicalAnalysis
    try
        strcmp(options.graphicalAnalysis,'True')
        warning('graphicalAnalysis not available for metabolite based solutions at this time. Only for food items or categories')
        options.graphicalAnalysis='False';
    catch
    end
    load('fdTable.mat')
end

model = changeObjective(model,obj);
model.osenseStr = objMinMax;

% Determine if any rois are metabolites
metRois=[];
for i=1:length(rois)
    if any(strcmp(model.mets,rois{i}))
        metRois=[metRois,i];
        %         if strcmp(slnType,'Detailed')
        %             slnType='Quick';
        %             disp('slnType changed to Quick because one or more rois defined as a metabolite')
        %         end
        if strcmp(roisMinMax{i},'max')
            model=addDemandReaction(model,rois{i}); %adds demand reaction as 'DM_metabolite'
            rois{i}=['DM_',rois{i}];
            model=changeRxnBounds(model,rois{i},100000,'u');
        else
            model=addSinkReactions(model,rois(i),-100000,0);
            rois{i}=['sink_',rois{i}];
        end
    end
end

%initialize optional variables
roiWeights=10*ones(1,length(rois));
initObjSln=[];
weightedFoodItems={};
caloricRange=[0 1e6];
slnType='Detailed';
roiBound='Bounded';
foodAddedLimit=1000000;
foodRemovedLimit=1000000;
calorieWeight='False';
freeMets={};
try
    if strcmp(options.slnType,'Detailed')
        graphicalAnalysis='True';
    else
        graphicalAnalysis='False';
    end
catch
    graphicalAnalysis='True';
end

if exist('options','var')
    fn = fieldnames(options);
    for k=1:numel(fn)
        %         if( isnumeric(options.(fn{k})) )
        %             % do stuff
        %         end
        if strcmp(fn{k},'roiWeights')
            roiWeights=options.roiWeights;
            if length(roiWeights)~=length(rois)
                error('length of roiWeights vector must be the same as items rois')
            end
        elseif strcmp(fn{k},'foodOrMets')
            foodOrMets=options.foodOrMets;
        elseif strcmp(fn{k},'calorieWeight')
            try
                foodOrMets=options.foodOrMets;
            catch
            end
            if strcmp(foodOrMets,'Food Items') || strcmp(foodOrMets,'Food Cat')
                calorieWeight=options.calorieWeight;
            else
                if strcmp(options.calorieWeight,'True') || strcmp(options.calorieWeight,'true')
                    error('"calorieWeight" cannot be True for metabolite solutions')
                end
            end
        elseif strcmp(fn{k},'roiBound')
            roiBound=options.roiBound;
            if ~strcmp(roiBound,'Unbounded') && ~strcmp(roiBound,'Bounded')
                error('Invalid roiBound input. Must be "Unbounded" or "Bounded"')
            end
        elseif strcmp(fn{k},'graphicalAnalysis')
            graphicalAnalysis=options.graphicalAnalysis;
            if ~strcmp(graphicalAnalysis,'True') && ~strcmp(graphicalAnalysis,'False')
                error('Invalid graphicalAnalysis input. Must be "True" or "False"')
            end
        elseif strcmp(fn{k},'foodAddedLimit')
            foodAddedLimit=options.foodAddedLimit;
        elseif strcmp(fn{k},'foodRemovedLimit')
            foodRemovedLimit=options.foodRemovedLimit;
        elseif strcmp(fn{k},'initObjSln')
            initObjSln=options.initObjSln;
        elseif strcmp(fn{k},'freeMets')
            freeMets=options.freeMets;
        elseif strcmp(fn{k},'caloricRange')
            caloricRange=options.caloricRange;
        elseif strcmp(fn{k},'slnType')
            slnType=options.slnType;
            if ~strcmp(slnType,'Detailed') && ~strcmp(slnType,'Quick')
                error('Invalid slnType input. Must be "Detailed" or "Quick"')
            end
        elseif strcmp(fn{k},'weightedFoodItems')
            weightedFoodItems=options.weightedFoodItems;
            if any(contains(options.weightedFoodItems(:,1),'Any_'))
                for i=length(options.weightedFoodItems(:,1)):-1:1
                    word=strsplit(options.weightedFoodItems{i,1},'Any_');
                    if isempty(word{1}) %If the word starts with 'Any_'
                        tmp=fdTable.Properties.VariableNames(contains(fdTable.Properties.VariableNames,word{2})).';
                        tmp=[tmp,num2cell(weightedFoodItems{i,2}*ones(length(tmp),1))];
                        weightedFoodItems(i,:)=[];
                        weightedFoodItems=[weightedFoodItems;tmp];
                    end
                end
            end
        else
            error(['Invalid "options" field entered: ', fn{k}])
        end
    end
end

if any(roiWeights<=0)
    error('"roiWeights" variable must be greater than zero')
end

%adjust ub and lb if roiBound specifies 'Unbound'
if strcmp(roiBound, 'Unbounded')
    for i=1:length(rois)
        f=find(strcmp(model.rxns,rois{i}));
        if strcmp(roisMinMax{i},'max')
            if model.ub(f)~=0
                model.ub(f)=100000;
            end
        else
            if model.lb(f)~=0
                model.lb(f)=-100000;
            end
        end
    end
end

for i=1:length(freeMets)
    str=['Diet_EX_',freeMets{i},'[d]'];
    try
        if strcmp(freeMets{i},'h2o')
            model = changeRxnBounds(model, str, -1e7, 'l');
        else
            model = changeRxnBounds(model, str, -1e5, 'l');
        end
    catch
        error(['Invalid metabolite specified in freeMets: ', freeMets{i}])
    end
end

newDietModel=model; %Copy original instance of model for new diet pointsModel
pointsModel=model; %Copy original instance of model for points pointsModel

%Calculate newDietModel objective function and restrict obj in main pointsModel
objIndex=find(contains(model.rxns,obj));
%
% roiIndex=find(strcmp(newDietModel.rxns,roi));
% disp(['Reaction of Interest = ', newDietModel.rxns{roiIndex}])
roiIndexO=zeros(1,length(rois));
for i=1:length(roiIndexO)
    roiIndexO(i)=find(strcmp(newDietModel.rxns,rois{i}));
    disp(['Reaction of Interest ', num2str(i),' = ', newDietModel.rxns{roiIndexO(i)}])
end

% get flux of objective function
if ~isempty(initObjSln)
    model_Obj=initObjSln;
    f1=model_Obj.f;
    initRoiFlux=model_Obj.v(roiIndexO);
elseif model.ub(objIndex)~=model.lb(objIndex)
    model_Obj = optimizeWBModel(newDietModel);
    f1=model_Obj.f;
    initRoiFlux=model_Obj.v(roiIndexO);
else
    f1=model.ub(objIndex);
    initRoiFlux=NaN(1,length(rois));
end

% If sln type is detailed, check if roi is already min or maxed out and
% if not, define min max range for roi

if strcmp(slnType,'Detailed')
    for i=1:length(rois)
        if initRoiFlux(i)==newDietModel.lb(roiIndexO)
            OroiFluxMin(i)=newDietModel.lb(roiIndexO(i));
        else
            pointsModel = changeObjective(pointsModel,rois{i});
            pointsModel.osenseStr = 'min';
            sln = optimizeWBModel(pointsModel);
            OroiFluxMin(i)=sln.v(roiIndexO(i));
        end
        if initRoiFlux(i)==newDietModel.ub(roiIndexO)
            OroiFluxMax(i)=newDietModel.ub(roiIndexO(i));
        else
            pointsModel = changeObjective(pointsModel,rois{i});
            pointsModel.osenseStr = 'max';
            sln = optimizeWBModel(pointsModel);
            OroiFluxMax(i)=sln.v(roiIndexO(i));
        end
    end
end
pointsModel=addMetabolite(pointsModel, 'unitOfFoodAdded[dP]');
pointsModel=addMetabolite(pointsModel, 'unitOfFoodRemoved[dP]');
pointsModel=addMetabolite(pointsModel, 'unitOfFoodChange[dP]');
pointsModel=addMetabolite(pointsModel, 'roiPoint[roiP]');
pointsModel=addMetabolite(pointsModel, 'point[P]');

fdTableMod=fdTable;
pro_Dindex=find(contains(fdTableMod.Var1,'pro_D')); %for now, remove pro_D from table.
fdTableMod(pro_Dindex,:)=[];
% eIndex=find(contains(fdTableMod.Var1,'Energy in Kcal'));  %remove Energy row from table.
% fdTableMod(eIndex,:)=[];

%Now add row for unitOfFoodAdded variable
tableChange=num2cell(ones(1,length(fdTableMod.Properties.VariableNames)-1));
tableChange=['unitOfFoodAdded',tableChange];
tableChange=cell2table(tableChange,'VariableNames',fdTableMod.Properties.VariableNames);
fdTableMod=[fdTableMod;tableChange];

%Modify unitOfFoodAdded points based on weight from weightedFoodItems
if ~isempty(weightedFoodItems)
    fdTableMod(end, fdTableMod.Properties.VariableNames(weightedFoodItems(:,1)))=(weightedFoodItems(:,2).');
end

if strcmp(foodOrMets,'Food Items') || strcmp(foodOrMets,'Food Cat')
    %Specify all food reactions and food metabolites
    if strcmp(calorieWeight,'True')
        fdTableMod{end,2:end}=fdTableMod{end-1,2:end}/206.28;
    end
    foodRxns=fdTableMod.Properties.VariableNames(2:end);
    foodRxns=strcat('Food_Added_EX_',foodRxns);
    foodRxns=strcat(foodRxns,'[d]');
    foodMetabolites= fdTableMod.Var1;
    foodMetabolites= strcat(foodMetabolites.','[d]');
    uof=find(contains(foodMetabolites,'unitOfFoodAdded'));
    foodMetabolites{uof}='unitOfFoodAdded[dP]';
    %Include food added reaction to pointsModel
    sMatrix=-1*table2array(fdTableMod(1:length(foodMetabolites),2:end));
    pointsModel = addMultipleReactions(pointsModel, foodRxns, foodMetabolites, sMatrix, 'lb', -100000*ones(1,length(foodRxns)), 'ub', zeros(1,length(foodRxns)));
    pointsModel = addMultipleReactions(pointsModel, {'Point_EX_unitOfFoodRemoved2Change[dp]','Point_EX_unitOfFoodAdded2Change[dp]','Point_EX_unitOfFoodChange[dP]_[P]','Point_EX_Point[P]'}, {'unitOfFoodRemoved[dP]','unitOfFoodAdded[dP]','unitOfFoodChange[dP]','point[P]'}, [-1 0 0 0;0 -1 0 0;1 1 -1 0;0 0 1 -1], 'lb', [-1000000,-1000000, -1000000,-1000000], 'ub', [foodRemovedLimit,foodAddedLimit,1000000,1000000]);
    % pointsModel = addMultipleReactions(pointsModel, {'Point_EX_unitOfFoodChange[dP]_[P]','Point_EX_Point[P]','Excretion_EX_Energy'}, {'unitOfFoodChange[dP]','point[P]','Energy in Kcal[d]'}, [-1 0 0;1 -1 0; 0 0 -1], 'lb', [-1000000,-1000000,caloricRange(1)], 'ub', [1000000,1000000,caloricRange(2)]);
elseif strcmp(foodOrMets,'AllMets') || strcmp(foodOrMets,'allMets')
    foodRxns= find(contains(pointsModel.rxns,'Diet_EX_'));
    foodMetabolites= foodRxns;
    foodMetabolites=regexprep(pointsModel.rxns(foodMetabolites),'Diet_EX_','');
    foodRxns=pointsModel.rxns(foodRxns);
    foodRxns=regexprep(foodRxns,'Diet_EX_','Food_Added_EX_');
    foodMetabolites{end+1}='unitOfFoodAdded[dP]';

    %Include food added reaction to pointsModel
    sMatrix=-1*eye(length(foodRxns));
    sMatrix=[sMatrix;-1*ones(1,length(foodRxns))];
    pointsModel = addMultipleReactions(pointsModel, foodRxns, foodMetabolites, sMatrix, 'lb', -100000*ones(1,length(foodRxns)), 'ub', zeros(1,length(foodRxns)));
    pointsModel = addMultipleReactions(pointsModel, {'Point_EX_unitOfFoodAdded2Change[dp]','Point_EX_unitOfFoodChange[dP]_[P]','Point_EX_Point[P]'}, {'unitOfFoodAdded[dP]','unitOfFoodChange[dP]','point[P]'}, [-1 0 0;1 -1 0;0 1 -1], 'lb', [-1000000, -1000000,-1000000], 'ub', [foodAddedLimit,1000000,1000000]);
elseif strcmp(foodOrMets,'FoodMets') || strcmp(foodOrMets,'foodMets')
    foodMetabolites= fdTableMod.Var1;
    foodMetabolites= strcat(foodMetabolites.','[d]');
    uof=find(contains(foodMetabolites,'unitOfFoodAdded'));
    foodMetabolites{uof}='unitOfFoodAdded[dP]';
    foodRxns=foodMetabolites;
    foodRxns(uof)=[];
    f=find(contains(foodMetabolites,'Energy_in_Kcal'));
    foodRxns(f)=[];
%     foodRxns=regexprep(foodRxns,'Diet_EX_','Food_Added_EX_');
    foodRxns=strcat('Food_Added_EX_',foodRxns);
    foodMetabolites(f)=[];
    %     foodRxns= strcat('Food_Added_EX_',foodRxns);
    %Include food added reaction to pointsModel
    sMatrix=-1*eye(length(foodRxns));
    sMatrix=[sMatrix;-1*ones(1,length(foodRxns))];
    pointsModel = addMultipleReactions(pointsModel, foodRxns, foodMetabolites, sMatrix, 'lb', -100000*ones(1,length(foodRxns)), 'ub', zeros(1,length(foodRxns)));
    pointsModel = addMultipleReactions(pointsModel, {'Point_EX_unitOfFoodAdded2Change[dp]','Point_EX_unitOfFoodChange[dP]_[P]','Point_EX_Point[P]'}, {'unitOfFoodAdded[dP]','unitOfFoodChange[dP]','point[P]'}, [-1 0 0;1 -1 0;0 1 -1], 'lb', [-1000000, -1000000,-1000000], 'ub', [foodAddedLimit,1000000,1000000]);
else
    error('Invalid foodOrMets specification')
end

uof=find(contains(foodMetabolites,'unitOfFoodAdded[dP]'));
foodMetabolites{uof}='unitOfFoodRemoved[dP]';
if strcmp(foodOrMets,'Food Items') || strcmp(foodOrMets,'Food Cat')
    %Identify food items already in the diet
    foodDietIndex=find(contains(pointsModel.rxns,'Food_EX_'));
    foodInDietIndex=foodDietIndex(pointsModel.lb(foodDietIndex)<0);
    %Build reaction matrix for food to be removed from diet
    sMatrix=zeros(length(foodMetabolites),length(foodInDietIndex));
    for i=1:length(foodInDietIndex)
        foodItem = regexprep(pointsModel.rxns{foodInDietIndex(i)},'Food_EX_','');
        foodRxn=strcat('Food_Removed_EX_',foodItem);
%         foodItem = regexprep(foodItem,'\[d\]','');
%         sMatrix(:,i)=fdTableMod.(foodItem);
%         sMatrix(end,i)=-1/sMatrix(end,i);
%         sMatrix(:,i)
        RxnFormula=printRxnFormula(pointsModel,pointsModel.rxns{foodInDietIndex(i)},0);
        [metaboliteList, stoichCoeffList, ~]=parseRxnFormula(RxnFormula{1});
        metaboliteList=[metaboliteList,{'unitOfFoodRemoved[dP]'}];
        stoichCoeffList=-1*stoichCoeffList;
        stoichCoeffList=[stoichCoeffList,-1];
        pointsModel = addMultipleReactions(pointsModel, {foodRxn}, metaboliteList, stoichCoeffList, 'lb', pointsModel.ub(foodInDietIndex(i)), 'ub', 0);
    end

%     pointsModel = addMultipleReactions(pointsModel, foodRxns, foodMetabolites, sMatrix, 'lb', pointsModel.ub(foodInDietIndex), 'ub', zeros(1,length(foodRxns)));
else

    %Build reaction matrix for food to be removed from diet
    sMatrix=-1*sMatrix;
    sMatrix(end,:)=-1./sMatrix(end,:);
    foodRxns = regexprep(foodRxns,'Food_Added_EX_','Food_Removed_EX_');
    pointsModel = addMultipleReactions(pointsModel, foodRxns, foodMetabolites, sMatrix, 'lb', -1000000*ones(1,length(foodRxns)), 'ub', zeros(1,length(foodRxns)));

end

%note that previous line makes food removal lb equal to food consumption ub
%as more food should not be able to be removed than the food consumption
%capabilities. If removing all of a food item yeilds the optimum solution,
%the diet consumption of said food item should hug the ub and food removal
%will then be equivalent to said ub yeilding 0.

% Introduce any sink or demand reactions if necessary
% if ~isempty(metRois)
%     for i=length(metRois):-1:1
%         if strcmp(roisMinMax{metRois(i)},'max')
%             pointsModel = addMultipleReactions(pointsModel, {strcmp(rois{metRois(i)},'_Demand')}, {rois{metRois(i)},'roiPoint[roiP]'}, [-1; -1*roiWeights(metRois(i))], 'lb', 0, 'ub', 1000000);
%             rois{metRois(i)}=[];
%             roiIndexO(metRois(i))=[];
%         else
%             pointsModel = addMultipleReactions(pointsModel, {strcmp(rois{metRois(i)},'_Sink')}, {rois{metRois(i)},'roiPoint[roiP]'}, [-1; roiWeights(metRois(i))], 'lb', -1000000, 'ub', 0);
%             rois{metRois(i)}=[];
%             roiIndexO(metRois(i))=[];
%         end
%     end
% end

%Get roi Indexes
for i=1:length(rois)
    roiIndexP(i)=find(strcmp(pointsModel.rxns,rois{i}));
end
roiUB=pointsModel.ub(roiIndexP);
roiLB=pointsModel.lb(roiIndexP);

%replace roi function
stoich=pointsModel.S(:,roiIndexP);
if length(roiIndexP)>1
    metInd=find(any(stoich.'~=0));
    metsRoi=pointsModel.mets(any(stoich.'~=0)).';
    metsStoich=full(stoich(any(stoich.'~=0),:));
else
    metsRoi=pointsModel.mets(find(stoich~=0)).';
    metsStoich=full(stoich(find(stoich~=0)));
end

weightVector=zeros(1,length(roiIndexP));
weightVector(contains(roisMinMax,'max'))=-1;
weightVector(contains(roisMinMax,'min'))=1;
metsStoich=[metsStoich;weightVector.*roiWeights;zeros(1,length(roiIndexP))];

for i=1:length(rois)
    evalc('[pointsModel,~,~]= removeRxns(pointsModel, rois{i})');
% [pointsModel,~,~]= removeRxns(pointsModel, rois{i});
end

pointsModel = addMultipleReactions(pointsModel, [rois,'Point_EX_roiPoints[roiP]_[P]'], [metsRoi,'roiPoint[roiP]','point[P]'], [metsStoich,[zeros(length(metsStoich(:,1))-2,1);-1;1]], 'lb', [roiLB.',-1000000], 'ub', [roiUB.',1000000]);
caloriesRxn=find(contains(pointsModel.rxns,'EX_DietEnergy'));
pointsModel.lb(caloriesRxn)=caloricRange(1);
pointsModel.ub(caloriesRxn)=caloricRange(2);

%Find solution
pointsModel = changeObjective(pointsModel,'Point_EX_Point[P]');
pointsModel.osenseStr = 'min';
pointsModelSln = optimizeWBModel(pointsModel);
% pointsModelSln.v(roiIndex)

disp(['Solution points =',num2str(pointsModelSln.f)])
disp([num2str(pointsModelSln.v(find(strcmp(pointsModel.rxns,'Point_EX_unitOfFoodChange[dP]_[P]')))),' come from diet']);
disp([num2str(pointsModelSln.v(find(strcmp(pointsModel.rxns,'Point_EX_roiPoints[roiP]_[P]')))),' come from roi']);
foodAddedIndexes=find(contains(pointsModel.rxns,'Food_Added_EX_'));
foodRemovedIndexes=find(contains(pointsModel.rxns,'Food_Removed_EX_'));
slnIndexes1=foodAddedIndexes(pointsModelSln.v(foodAddedIndexes)<0);
slnIndexes2=foodRemovedIndexes(pointsModelSln.v(foodRemovedIndexes)<0);

disp('Food items of interest are:')
T=table([pointsModel.rxns(slnIndexes1);pointsModel.rxns(slnIndexes2)],pointsModelSln.v([slnIndexes1;slnIndexes2]),'VariableNames',{'Food Rxn', 'Flux'})
if ~strcmp(foodOrMets,'AllMets') && ~strcmp(foodOrMets,'FoodMets')
    disp(['Diet Energy = ',num2str(pointsModelSln.v(contains(pointsModel.rxns,'EX_DietEnergy')))])
end

%Add and remove relevant food items from diet in newDietModel
if strcmp(foodOrMets,'Food Items') || strcmp(foodOrMets,'Food Cat')
    foodItemsAdd= regexprep(pointsModel.rxns(slnIndexes1),'Food_Added_EX_','Food_EX_');
    foodItemsRemove= regexprep(pointsModel.rxns(slnIndexes2),'Food_Removed_EX_','Food_EX_');
else
    foodItemsAdd= regexprep(pointsModel.rxns(slnIndexes1),'Food_Added_EX_','Diet_EX_');
    foodItemsRemove= regexprep(pointsModel.rxns(slnIndexes2),'Food_Removed_EX_','Diet_EX_');
end
modelOindexAdd=zeros(1,length(foodItemsAdd));
sl2IndexAdd=zeros(1,length(foodItemsAdd));
modelOindexRemove=zeros(1,length(foodItemsRemove));
sl2IndexRemove=zeros(1,length(foodItemsRemove));
for i=1:length(foodItemsAdd)
    modelOindexAdd(i)=find(contains(newDietModel.rxns,foodItemsAdd(i)));
    sl2IndexAdd(i)=find(contains(pointsModel.rxns,foodItemsAdd(i)));
end
for i=1:length(foodItemsRemove)
    modelOindexRemove(i)=find(contains(newDietModel.rxns,foodItemsRemove(i)));
    sl2IndexRemove(i)=find(contains(pointsModel.rxns,foodItemsRemove(i)));
end
newDietModel.lb(modelOindexAdd)=(pointsModelSln.v(sl2IndexAdd)+pointsModelSln.v(slnIndexes1))*1.01;
newDietModel.ub(modelOindexAdd)=(pointsModelSln.v(sl2IndexAdd)+pointsModelSln.v(slnIndexes1))*0.99;
newDietModel.lb(modelOindexRemove)=(pointsModelSln.v(sl2IndexRemove)-pointsModelSln.v(slnIndexes2))*1.01;
newDietModel.ub(modelOindexRemove)=(pointsModelSln.v(sl2IndexRemove)-pointsModelSln.v(slnIndexes2))*0.99;

if strcmp(graphicalAnalysis,'True')
    [ogCatTable] = getDietComposition(model,'Off');
    [newCatTable] = getDietComposition(newDietModel,'Off');
    OgMacros=ogCatTable.('Mass (g)');
    NewCategories=newCatTable.('Category');
    NewMacros=newCatTable.('Mass (g)');
    figure()
    t = tiledlayout(1,2,'TileSpacing','compact');
    ax1 = nexttile;
    pie(ax1,OgMacros(1:4))
    title('Original Diet')
    ax2 = nexttile;
    pie(ax2,NewMacros(1:4))
    title('New Diet')
    legend(NewCategories(1:4),'Location','South')
end

menuChanges=T;
if strcmp(slnType,'Quick')
    for i=1:length(rois)
        ind=find(strcmp(pointsModel.rxns,rois{i}));
        disp([rois{i},' flux = ', num2str(pointsModelSln.v(ind))])
    end
    roiFlux(i)=pointsModelSln.v(ind);
    slnRanges=pointsModelSln.v(roiIndexP);
    %     menuChanges.itemsRemoved=[];
    %     menuChanges.itemsAdded=[];
    %     menuChanges.fullMenu=[];
    newDietModel = changeObjective(newDietModel,obj);
    newDietModel.osenseStr = objMinMax;
    return
end

%Find new obj flux with new diet

if model.ub(objIndex)~=model.lb(objIndex)
    model_Obj = optimizeWBModel(newDietModel);
    f2=model_Obj.f;
    newDietModel=changeRxnBounds(newDietModel,obj,f2,'b'); %constrain pointsModel obj flux
else
    f2=model.ub(objIndex);
end

disp(['f1 =',num2str(f1), ' & f2=', num2str(f2)])

%Compute new min max ranges for roi with new diet
%%
for i=1:length(rois)
    disp(rois{i})
    newDietModel = changeObjective(newDietModel,rois{i});
    newDietModel.osenseStr = 'min';
    tmp=optimizeWBModel(newDietModel);
    slnMin.(['Rxn',num2str(i)]) = tmp;
    NroiFluxMin(i)=slnMin.(['Rxn',num2str(i)]).v(roiIndexO(i));
    newDietModel.osenseStr = 'max';
    tmp=optimizeWBModel(newDietModel);
    slnMax.(['Rxn',num2str(i)]) = tmp;
    NroiFluxMax(i)=slnMax.(['Rxn',num2str(i)]).v(roiIndexO(i));
    disp(['Original Diet RoI range = ', num2str(OroiFluxMin(i)), ':', num2str(OroiFluxMax(i))])
    disp(['New Diet RoI range = ', num2str(NroiFluxMin(i)), ':', num2str(NroiFluxMax(i))])
end

% slnMin
roiFlux=[slnMin.',slnMax.'];
newDietModel.ub(objIndex)=model.ub(objIndex);
newDietModel.lb(objIndex)=model.lb(objIndex);
newDietModel = changeObjective(newDietModel,obj);
newDietModel.osenseStr = objMinMax;

disp('___________________________________________________________________')

end

%TO DO:
% -Impliment itemsRemoved and itemsAdded
% -Introduce option for metabolite adjustments
% -Include intelligent default weighting for roi
% -slnMin & slnMax turn to structs or convert to flux values

%Finished:
%-introduced Quick/Detailed functionality
%-renamed variables to be more obvious
%-commented code

% figure()
% t = tiledlayout(2,2,'TileSpacing','compact');
% ax1 = nexttile;
% pie(ax1,OgMacros(1:4))
% title('Original Diet')
% ax2 = nexttile;
% x=pie(ax2,NewMacros(1:4))
% title('New Diet')
% legend(NewCategories(1:4),'Location','South')
% ax3 = nexttile;
% bar(ax3,(NewMacros(1:4)-OgMacros(1:4))./OgMacros(1:4))
% ylim([-1 14])
% ax4= nexttile;
% legend(x,NewCategories(1:4))

author: opencobra

repo: cobratoolbox


text:
%% Machine Learning Online Class
%  Exercise 6 | Support Vector Machines
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     gaussianKernel.m
%     dataset3Params.m
%     processEmail.m
%     emailFeatures.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% =============== Part 1: Loading and Visualizing Data ================
%  We start the exercise by first loading and visualizing the dataset.
%  The following code will load the dataset into your environment and plot
%  the data.
%

fprintf('Loading and Visualizing Data ...\n')

% Load from ex6data1:
% You will have X, y in your environment
load('ex6data1.mat');

% Plot training data
plotData(X, y);

fprintf('Program paused. Press enter to continue.\n');
%pause;

%% ==================== Part 2: Training Linear SVM ====================
%  The following code will train a linear SVM on the dataset and plot the
%  decision boundary learned.
%

% Load from ex6data1:
% You will have X, y in your environment
load('ex6data1.mat');

fprintf('\nTraining Linear SVM ...\n')

% You should try to change the C value below and see how the decision
% boundary varies (e.g., try C = 1000)
C = 100;
model = svmTrain(X, y, C, @linearKernel, 1e-3, 20);
visualizeBoundaryLinear(X, y, model);

fprintf('Program paused. Press enter to continue.\n');
%pause;

%% =============== Part 3: Implementing Gaussian Kernel ===============
%  You will now implement the Gaussian kernel to use
%  with the SVM. You should complete the code in gaussianKernel.m
%
fprintf('\nEvaluating the Gaussian Kernel ...\n')

x1 = [1 2 1]; x2 = [0 4 -1]; sigma = 2;
sim = gaussianKernel(x1, x2, sigma);

fprintf(['Gaussian Kernel between x1 = [1; 2; 1], x2 = [0; 4; -1], sigma = 0.5 :' ...
         '\n\t%f\n(this value should be about 0.324652)\n'], sim);

fprintf('Program paused. Press enter to continue.\n');
%pause;

%% =============== Part 4: Visualizing Dataset 2 ================
%  The following code will load the next dataset into your environment and
%  plot the data.
%

fprintf('Loading and Visualizing Data ...\n')

% Load from ex6data2:
% You will have X, y in your environment
load('ex6data2.mat');

% Plot training data
plotData(X, y);

fprintf('Program paused. Press enter to continue.\n');
%pause;

%% ========== Part 5: Training SVM with RBF Kernel (Dataset 2) ==========
%  After you have implemented the kernel, we can now use it to train the
%  SVM classifier.
%
fprintf('\nTraining SVM with RBF Kernel (this may take 1 to 2 minutes) ...\n');

% Load from ex6data2:
% You will have X, y in your environment
load('ex6data2.mat');

% SVM Parameters
C = 1; sigma = 0.1;

% We set the tolerance and max_passes lower here so that the code will run
% faster. However, in practice, you will want to run the training to
% convergence.
model= svmTrain(X, y, C, @(x1, x2) gaussianKernel(x1, x2, sigma));
visualizeBoundary(X, y, model);

fprintf('Program paused. Press enter to continue.\n');
%pause;

%% =============== Part 6: Visualizing Dataset 3 ================
%  The following code will load the next dataset into your environment and
%  plot the data.
%

fprintf('Loading and Visualizing Data ...\n')

% Load from ex6data3:
% You will have X, y in your environment
load('ex6data3.mat');

% Plot training data
plotData(X, y);

fprintf('Program paused. Press enter to continue.\n');
%pause;

%% ========== Part 7: Training SVM with RBF Kernel (Dataset 3) ==========

%  This is a different dataset that you can use to experiment with. Try
%  different values of C and sigma here.
%

% Load from ex6data3:
% You will have X, y in your environment
load('ex6data3.mat');

% Try different SVM Parameters here
[C, sigma] = dataset3Params(X, y, Xval, yval);

% Train the SVM
model= svmTrain(X, y, C, @(x1, x2) gaussianKernel(x1, x2, sigma));
visualizeBoundary(X, y, model);

fprintf('Program paused. Press enter to continue.\n');
pause;

author: yhyap

repo: machine-learning-coursera


text:
%CODEGENERATOR.GENSLBLOCKJACOBIAN Generate Simulink block for robot Jacobians
%
% cGen.genslblockjacobian() generates a robot-specific Simulink block to compute
% robot Jacobians (world and tool frame).
%
% Notes::
% - Is called by CodeGenerator.genjacobian if cGen has active flag genslblock
% - The Simulink blocks are generated and stored in a robot specific block
%   library cGen.slib in the directory cGen.basepath.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.genjacobian.

% Copyright (C) 2012-2014, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function genslblockjacobian(CGen)

%% Open or create block library
bdclose('all')                                                              % avoid problems with previously loaded libraries
load_system('simulink');
if ~(exist([CGen.slibpath,simulinkext]) == 2)                                  % Create new block library if none exists
 CGen.createnewblocklibrary;
end
open_system(CGen.slibpath);
set_param(CGen.slib,'lock','off');

q = CGen.rob.gencoords;
%% Jacobian0
CGen.logmsg([datestr(now),'\tGenerating jacobian Embedded Matlab Function Block with respect to the robot base frame']);
%     [datestr(now),'\tGenerating jacobian Embedded Matlab Function Block with respect to the end-effector frame']);
symname = 'jacob0';
fname = fullfile(CGen.sympath,[symname,'.mat']);

if exist(fname,'file')
    tmpStruct = load(fname);
else
    error ('genSLBlockFkine:SymbolicsNotFound','Save symbolic expressions to disk first!')
end

blockaddress = [CGen.slib,'/',symname];          % treat intermediate transformations separately
if ~isempty(find_system(CGen.slib,'SearchDepth',1,'Name',symname))                    % Delete previously generated block
    delete_block(blockaddress);
    save_system;
end

symexpr2slblock(blockaddress,tmpStruct.(symname),'vars',{q});

CGen.logmsg('\t%s\n',' done!');

%% Jacobe
CGen.logmsg([datestr(now),'\tGenerating jacobian Embedded Matlab Function Block with respect to the end-effector frame']);
symname = 'jacobe';
fname = fullfile(CGen.sympath,[symname,'.mat']);

if exist(fname,'file')
    tmpStruct = load(fname);
else
    error ('genSLBlockFkine:SymbolicsNotFound','Save symbolic expressions to disk first!')
end

blockaddress = [CGen.slib,'/',symname];          % treat intermediate transformations separately
if ~isempty(find_system(CGen.slib,'SearchDepth',1,'Name',symname))                    % Delete previously generated block
    delete_block(blockaddress);
    save_system;
end

symexpr2slblock(blockaddress,tmpStruct.(symname),'vars',{q});
CGen.logmsg('\t%s\n',' done!');

%% Cleanup
% Arrange blocks
distributeblocks(CGen.slib);

% Lock, save and close library
set_param(CGen.slib,'lock','on');
save_system(CGen.slib,CGen.slibpath);
close_system(CGen.slib);

end

author: petercorke

repo: robotics-toolbox-matlab


text:
function val=fallingFactorial(x,n)
%%FALLINGFACTORIAL Evaluate the falling factorial (x)_n, which for n>=0 is
%                  also known as the binomial polynomial, the lower
%                  factorial, the factorial power, and the falling
%                  factorial power. It is defined
%                  (x)_n=x*(x-1)*(x-2)*...*(x-(n-1)). It can be expressed
%                  as gamma(x+1)/gamma(x-n+1), which makes it continuous.
%                  It is implemented using logarthmic functions to reduce
%                  problems with overflows in intermediate results.
%
%INPUTS: x A scalar or matrix of first values in the product series of
%          x*(x-1)*...*(x-(n-1)). x cannot be negative.
%        n The scalar or matrix n value in the product series of
%          x*(x-1)*...*(x-(n-1)). If x is a matrix, n can be a scalar or a
%          matrix of the same size. If x is a scalar, n can be a scalar or
%          a matrix.
%
%OUTPUTS: val The falling factorial (x)_n for each x and the corresponding
%             n/ for each n and a single x depending on what is a matrix
%             and what is scalar. Note that if n>x, the falling factorial
%             is defined to be 0.
%
%The falling factorial is discussed in [1]. Note that
%fallingFactorial(n,n)=factorial(n);
%
%REFERENCES:
%[1] Weisstein, Eric W. "Falling Factorial." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/FallingFactorial.html
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

val=exp(gammaln(x+1)-gammaln(max(0,x-n+1)));

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 [geo] = GetGeometryfromGUI(root, tip, wing, cruise, engine, hplot, output, airfoildata)
% Get parameters from the GUI and create geometry structure "geo"
% Gets the user defined parameters from the GUI
% Performs limited error checking and sets boolean to tell ExecCalc to continue or not

%Clear clock
set(output.calctime,'String','----');

%Set wing geometry from values from GUI
geo.t0 = clock; %Start time
geo.b = str2num(get(wing.span,'String'))/3.281;  %Convert to meters
geo.c_r = str2num(get(root.chord,'String'))/3.281;  %Convert to meters
geo.taper = str2num(get(tip.chord,'String'))/str2num(get(root.chord,'String'));
geo.i_r = str2num(get(root.angle,'String'))*pi/180;   %Convert to radians
geo.twist = (str2num(get(tip.angle,'String'))-str2num(get(root.angle,'String')))*pi/180;  %Convert to radians
geo.dih = str2num(get(wing.dihedral,'String'))*pi/180; %Convert to radians
geo.root = airfoildata{1}{get(root.airfoil,'Value')};
geo.tip = airfoildata{1}{get(tip.airfoil,'Value')};
geo.alpha = str2num(get(wing.AOA,'String'))*pi/180; %Convert to radians
geo.V = str2num(get(cruise.velocity,'String'))*0.5144; %Convert knots to m/s
[geo.density, geo.a geo.visc] = StandardAtmosphere(get(cruise.altitude,'String'));
geo.M = geo.V/geo.a;
geo.S_Sref = str2num(get(cruise.wettedarea,'String'));  %Dimensionless
geo.cf = str2num(get(cruise.skinfriction,'String'));  %Dimensionless
geo.ns = str2num(get(wing.ns,'String'));
geo.nc = str2num(get(wing.nc,'String'));
geo.sweep = str2num(get(wing.sweep,'String'))*pi/180;  %Convert to radians
geo.S = 0.5*(geo.c_r + geo.c_r*geo.taper)*geo.b;  %Square meters
geo.c_av = 0.5*(geo.c_r + geo.c_r*geo.taper);  %Meters
geo.Re_r = geo.V*geo.density*geo.c_r/geo.visc; %Reynolds number at root
geo.Re_t = geo.V*geo.density*geo.c_r*geo.taper/geo.visc; %Reynolds number at root
geo.rootindex = get(root.airfoil,'Value');      %index to selected root airfoil
geo.tipindex = get(tip.airfoil,'Value');         %index to selected tip airfoil
geo.propefficiency = str2num(get(engine.propeffic,'String'));  %from Anderson, Aircraft Performance and Design
geo.propSFC = str2num(get(engine.SFC,'String'))*1.657e-6; % convert lb fuel/hp-hr to m^-1 Wikipedia (Specific Fuel Consumption)
geo.jetTSFC = str2num(get(engine.TSFC,'String'))/3600; %convert to lb fuel per lb-s thrust; Anderson, Aircraft Performance and Design pg. 299
geo.emptyweight = str2num(get(cruise.weight,'String'))*4.44;%Convert weight from lbf to N
geo.fueldens = 0.840*1e3; %(kg/m^3) Density of Jet A-1 at 15deg C (Wikipedia)
geo.wingfuel = str2num(get(engine.wingfuel,'String'))/100;  %Percent of wing dedicated to carrying fuel
geo.withinconstraints = 1; % Set boolean flag stating that all geometry is within constraints

%Set the selected type of engine
if get(engine.panel,'SelectedObject')==engine.prop
    geo.engine = 'prop';
elseif get(engine.panel,'SelectedObject')==engine.jet
    geo.engine = 'jet';
else
    geo.engine = 0;
end

%Set initial output
set(wing.taper,'String',num2str(geo.taper));
set(wing.twist,'String',num2str(geo.twist));
set(root.Re,'String',num2str(round(geo.Re_r/1e3)/1e3));
set(tip.Re,'String',num2str(round(geo.Re_t/1e3)/1e3));
set(cruise.density,'String',num2str(round(geo.density*1000)/1000));
set(cruise.viscosity,'String',num2str(round(geo.visc*10000000)/10000000));
set(cruise.mach,'String',num2str(round(geo.M*100)/100));

%Limited error checking
if geo.sweep < 80*pi/180  %If sweep < 80 deg
    set(wing.sweep,'ForegroundColor','black')
else
    set(wing.sweep,'ForegroundColor','red')
    ZeroOutput(output)
    geo.withinconstraints = 0; % User defined geometry is not within constraints
end
if geo.dih < 80*pi/180  %If dihedral < 80 deg
    set(wing.dihedral,'ForegroundColor','black')
else
    set(wing.dihedral,'ForegroundColor','red')
    ZeroOutput(output)
    geo.withinconstraints = 0; % User defined geometry is not within constraints
end
if geo.M < 0.65  %If M < 0.65
    set(cruise.mach,'ForegroundColor','black')
else
    set(cruise.mach,'ForegroundColor','red')
    ZeroOutput(output)
    geo.withinconstraints = 0; % User defined geometry is not within constraints
end

author: Sable

repo: mcbench-benchmarks


text:
%compute velmag_head
function [data,units]=compute_velmaghead(trx,n)
larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
velmaghead=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    velmaghead{1,i}=bsxfun(@hypot,trx(larva).dxhead_mm,trx(larva).dyhead_mm);
end

units=parseunits('mm/s');
data=velmaghead;

author: kristinbranson

repo: JAABA


text:
function X = ankurAngle2XYZ(X)

if(size(X,2)~=54)
  error('Wrong Dimensionality of Data');
end

matrix_to_bvh(X','tmp.bvh');
X = bvh_to_3dmatrix('tmp.bvh');

% clean up
system('rm tmp.bvh');

return;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% testRadiomicsFeaturesWithPyRadiomics.m
%
% Script to compare CERR and pyRadiomics features
%
% APA, 2/22/2018

global planC
indexS = planC{end};

% Wavelet decomposition flag
wavDecompFlg = 1; % 0: original image, 1: wavelet pre-processing

% Specify discretization
fixedMinMaxGrLevFlag = 0; % 1: to input fixed min/max/grLevels, 0: min/max
                          % calculated for the structure and bin width of binWidth.
binWidth = 25;
structNum = 1;
scanNum = 1;
pyradScanType = 'original';

paramS.firstOrderParamS.offsetForEnergy = planC{indexS.scan}(scanNum).scanInfo(1).CTOffset;
paramS.firstOrderParamS.binWidth = binWidth;

% Get structure
[rasterSegments, planC, isError] = getRasterSegments(structNum,planC);
[mask3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
scanArray3M = getScanArray(planC{indexS.scan}(scanNum));
scanArray3M = double(scanArray3M) - ...
    planC{indexS.scan}(scanNum).scanInfo(1).CTOffset;
% Wavelet decomposition
if wavDecompFlg == 1
    dirString = 'HHL';
    wavType = 'coif1';
    scanArray3M = flip(scanArray3M,3);
    if mod(size(scanArray3M,3),2) > 0
        scanArray3M(:,:,end+1) = 0*scanArray3M(:,:,1);
    end
    scanArray3M = wavDecom3D(double(scanArray3M),dirString,wavType);
    if mod(size(scanArray3M,3),2) > 0
        scanArray3M = scanArray3M(:,:,1:end-1);
    end
    scanArray3M = flip(scanArray3M,3);
end

SUVvals3M = mask3M.*double(scanArray3M(:,:,uniqueSlices));
[minr, maxr, minc, maxc, mins, maxs] = compute_boundingbox(mask3M);
maskBoundingBox3M = mask3M(minr:maxr,minc:maxc,mins:maxs);

% Assign NaN to image outside mask
volToEval = SUVvals3M(minr:maxr,minc:maxc,mins:maxs);
volToEval(~maskBoundingBox3M) = NaN;

% Get x,y,z grid for the shape features (flip y to make it monotically
% increasing)
[xValsV, yValsV, zValsV] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));
yValsV = fliplr(yValsV);
xValsV = xValsV(minc:maxc);
yValsV = yValsV(minr:maxr);
zValsV = zValsV(mins:maxs);
PixelSpacingX = abs(xValsV(2)-xValsV(1));
PixelSpacingY = abs(yValsV(2)-yValsV(1));
PixelSpacingZ = abs(zValsV(2)-zValsV(1));
VoxelVol = PixelSpacingX*PixelSpacingY*PixelSpacingZ*1000; % convert cm to mm

if fixedMinMaxGrLevFlag
    numGrLevels = paramS.higherOrderParamS.numGrLevels;
    minIntensity = paramS.higherOrderParamS.minIntensity;
    maxIntensity = paramS.higherOrderParamS.maxIntensity;
    % Quantize using the number of bins
    quantizedM = imquantize_cerr(volToEval,numGrLevels,...
        minIntensity,maxIntensity);

else % Quantize using the binwidth
    minIntensity = [];
    maxIntensity = [];
    numGrLevels = [];
    quantizedM = imquantize_cerr(volToEval,numGrLevels,...
        minIntensity,maxIntensity,binwidth);
    numGrLevels = max(quantizedM(:));
    paramS.higherOrderParamS.numGrLevels = numGrLevels;
end

paramS.higherOrderParamS.numGrLevels = numGrLevels;
paramS.higherOrderParamS.patchRadius3dV = [1 1 1];
paramS.higherOrderParamS.imgDiffThresh = 0;

% Number of voxels (used in run percentage calculation)
numVoxels = sum(~isnan(quantizedM(:)));

%% First order features
featureS.firstOrderS = radiomics_first_order_stats...
    (volToEval(logical(maskBoundingBox3M)), VoxelVol, ...
    paramS.firstOrderParamS.offsetForEnergy, paramS.firstOrderParamS.binWidth);
firstOrderS = featureS.firstOrderS;
cerrFirstOrderV = [firstOrderS.energy, firstOrderS.totalEnergy, firstOrderS.interQuartileRange, ...
    firstOrderS.kurtosis+3, firstOrderS.max, firstOrderS.mean, firstOrderS.meanAbsDev, ...
    firstOrderS.median, firstOrderS.medianAbsDev, firstOrderS.min, ...
    firstOrderS.P10, firstOrderS.P90, firstOrderS.interQuartileRange, ...
    firstOrderS.robustMeanAbsDev, firstOrderS.rms, firstOrderS.skewness, ...
    firstOrderS.std, firstOrderS.var, firstOrderS.entropy];
pyradFirstorderNamC = {'Energy', 'TotalEnergy','InterquartileRange','Kurtosis',...
    'Maximum', 'Mean','MeanAbsoluteDeviation','Median','medianAbsDev',...
    'Minimum','10Percentile','90Percentile','InterquartileRange',...
    'RobustMeanAbsoluteDeviation','RootMeanSquared','Skewness',...
    'StandardDeviation','Variance','Entropy'};
pyradFirstorderNamC = strcat([pyradScanType,'_firstorder_'],pyradFirstorderNamC);
pyRadFirstOrderV = [];
for i = 1:length(pyradFirstorderNamC)
    if isfield(teststruct,pyradFirstorderNamC{i})
        pyRadFirstOrderV(i) = teststruct.(pyradFirstorderNamC{i});
    else
        pyRadFirstOrderV(i) = NaN;
    end
end

diffFirstOrderV = (cerrFirstOrderV - pyRadFirstOrderV) ./ cerrFirstOrderV * 100

%% Shape features
featureS.shapeS = getShapeParams(maskBoundingBox3M, ...
        {xValsV, yValsV, zValsV});
shapeS = featureS.shapeS;
cerrShapeV = [shapeS.majorAxis, shapeS.minorAxis, shapeS.leastAxis, ...
    shapeS.flatness, shapeS.elongation, shapeS.max3dDiameter, shapeS.max2dDiameterAxialPlane,...
    shapeS.max2dDiameterSagittalPlane', shapeS.max2dDiameterCoronalPlane, ...
    shapeS.Compactness1, shapeS.Compactness2, shapeS.spherDisprop, ...
    shapeS.sphericity, shapeS.surfToVolRatio/10,...
    shapeS.surfArea*100, shapeS.volume*1000];
pyradShapeNamC = {'MajorAxis', 'MinorAxis', 'LeastAxis', 'Flatness',  'Elongation', ...
    'Maximum3DDiameter', 'Maximum2DDiameterSlice', 'Maximum2DDiameterRow', ...
    'Maximum2DDiameterColumn', 'Compactness1','Compactness2','spherDisprop','Sphericity', ...
    'SurfaceVolumeRatio','SurfaceArea','Volume'};
pyradShapeNamC = strcat([pyradScanType,'_shape_'],pyradShapeNamC);
pyRadShapeV = [];
for i = 1:length(pyradShapeNamC)
    if isfield(teststruct,pyradShapeNamC{i})
        pyRadShapeV(i) = teststruct.(pyradShapeNamC{i});
    else
        pyRadShapeV(i) = NaN;
    end
end
shapeDiffV = (cerrShapeV - pyRadShapeV) ./ cerrShapeV * 100

%% GLCM features

% Define flags
glcmFlagS.energy = 1;
glcmFlagS.jointEntropy = 1;
glcmFlagS.jointMax = 1;
glcmFlagS.jointAvg = 1;
glcmFlagS.jointVar = 1;
glcmFlagS.contrast = 1;
glcmFlagS.invDiffMoment = 1;
glcmFlagS.sumAvg = 1;
glcmFlagS.corr = 1;
glcmFlagS.clustShade = 1;
glcmFlagS.clustProm = 1;
glcmFlagS.haralickCorr = 1;
glcmFlagS.invDiffMomNorm = 1;
glcmFlagS.invDiff = 1;
glcmFlagS.invDiffNorm = 1;
glcmFlagS.invVar = 1;
glcmFlagS.dissimilarity = 1;
glcmFlagS.diffEntropy = 1;
glcmFlagS.diffVar = 1;
glcmFlagS.diffAvg = 1;
glcmFlagS.sumVar = 1;
glcmFlagS.sumEntropy = 1;
glcmFlagS.clustTendency = 1;
glcmFlagS.autoCorr = 1;
glcmFlagS.invDiffMomNorm = 1;
glcmFlagS.firstInfCorr = 1;
glcmFlagS.secondInfCorr = 1;

dirctn      = 1;
cooccurType = 2;
featureS.harFeat3DdirS = get_haralick(dirctn, cooccurType, quantizedM, ...
numGrLevels, glcmFlagS);

% harlCombS = featureS.harFeat3DcombS.CombS;
harlCombS = featureS.harFeat3DdirS.AvgS;
cerrGlcmV = [harlCombS.autoCorr, harlCombS.jointAvg, harlCombS.clustPromin, harlCombS.clustShade, harlCombS.clustTendency, ...
harlCombS.contrast, harlCombS.corr, harlCombS.diffAvg, harlCombS.diffEntropy, harlCombS.diffVar, harlCombS.dissimilarity, ...
harlCombS.energy, harlCombS.jointEntropy, harlCombS.invDiff, harlCombS.invDiffMom, harlCombS.firstInfCorr, ...
harlCombS.secondInfCorr, harlCombS.invDiffMomNorm, harlCombS.invDiffNorm, harlCombS.invVar, ...
harlCombS.sumAvg, harlCombS.sumEntropy, harlCombS.sumVar];

pyradGlcmNamC = {'Autocorrelation', 'JointAverage', 'ClusterProminence', 'ClusterShade',  'ClusterTendency', ...
    'Contrast', 'Correlation', 'DifferenceAverage', 'DifferenceEntropy', 'DifferenceVariance', 'Dissimilarity', ...
    'JointEnergy', 'JointEntropy','Id','Idm', 'Imc1' , ...
    'Imc2', 'Idmn','Idn','InverseVariance', 'sumAverage', 'SumEntropy', 'sumVariance'};

pyradGlcmNamC = strcat([pyradScanType,'_glcm_'],pyradGlcmNamC);
pyRadGlcmV = [];
for i = 1:length(pyradGlcmNamC)
    if isfield(teststruct,pyradGlcmNamC{i})
        pyRadGlcmV(i) = teststruct.(pyradGlcmNamC{i});
    else
        pyRadGlcmV(i) = NaN;
    end
end

glcmDiffV = (cerrGlcmV - pyRadGlcmV) ./ cerrGlcmV * 100

%% RLM features

rlmFlagS.sre = 1;
rlmFlagS.lre = 1;
rlmFlagS.gln = 1;
rlmFlagS.glnNorm = 1;
rlmFlagS.rln = 1;
rlmFlagS.rlnNorm = 1;
rlmFlagS.rp = 1;
rlmFlagS.lglre = 1;
rlmFlagS.hglre = 1;
rlmFlagS.srlgle = 1;
rlmFlagS.srhgle = 1;
rlmFlagS.lrlgle = 1;
rlmFlagS.lrhgle = 1;
rlmFlagS.glv = 1;
rlmFlagS.rlv = 1;
rlmFlagS.re = 1;

numGrLevels = paramS.higherOrderParamS.numGrLevels;

% 3D Run-Length features from combined run length matrix
dirctn      = 1;
rlmType     = 2;
featureS.rlmFeat3DdirS = get_rlm(dirctn, rlmType, quantizedM, ...
    numGrLevels, numVoxels, rlmFlagS);

rlmCombS = featureS.rlmFeat3DdirS.AvgS;
cerrRlmV = [rlmCombS.gln, rlmCombS.glnNorm, rlmCombS.glv, rlmCombS.hglre, rlmCombS.lglre, rlmCombS.lre, rlmCombS.lrhgle, ...
    rlmCombS.lrlgle, rlmCombS.rln, rlmCombS.rlnNorm, rlmCombS.rlv, rlmCombS.rp, ...
    rlmCombS.sre, rlmCombS.srhgle, rlmCombS.srlgle, rlmCombS.re];

pyradRlmNamC = {'GrayLevelNonUniformity', 'GrayLevelNonUniformityNormalized',...
    'GrayLevelVariance', 'HighGrayLevelRunEmphasis',  'LowGrayLevelRunEmphasis', ...
    'LongRunEmphasis', 'LongRunHighGrayLevelEmphasis', 'LongRunLowGrayLevelEmphasis',...
    'RunLengthNonUniformity', 'RunLengthNonUniformityNormalized', 'RunVariance', ...
    'RunPercentage', 'ShortRunEmphasis','ShortRunHighGrayLevelEmphasis', ...
    'ShortRunLowGrayLevelEmphasis','RunEntropy'};

pyradRlmNamC = strcat([pyradScanType,'_glrlm_'],pyradRlmNamC);
pyRadRlmV = [];
for i = 1:length(pyradRlmNamC)
    if isfield(teststruct,pyradRlmNamC{i})
        pyRadRlmV(i) = teststruct.(pyradRlmNamC{i});
    else
        pyRadRlmV(i) = NaN;
    end
end

rlmDiffV = (cerrRlmV - pyRadRlmV) ./ cerrRlmV * 100

%% Size Zone features in 3d

flagS.sae = 1;
flagS.lae = 1;
flagS.gln = 1;
flagS.glv = 1;
flagS.szv = 1;
flagS.glnNorm = 1;
flagS.szn = 1;
flagS.sznNorm = 1;
flagS.zp = 1;
flagS.lglze = 1;
flagS.hglze = 1;
flagS.salgle = 1;
flagS.sahgle = 1;
flagS.lalgle = 1;
flagS.larhgle = 1;
flagS.ze = 1;

szmType = 1; % 1: 3d, 2: 2d
szmM = calcSZM(quantizedM, numGrLevels, szmType);
numVoxels = sum(~isnan(quantizedM(:)));
featureS.szmFeature3dS = szmToScalarFeatures(szmM,numVoxels, szmFlagS);
szmS = featureS.szmFeature3dS;
% cerrSzmV = [szmS.gln, szmS.glnNorm, szmS.glv, szmS.hglre, szmS.lglre, szmS.lre, szmS.lrhgle, ...
%     szmS.lrlgle, szmS.rln, szmS.rlnNorm, szmS.rlv, szmS.rp, ...
%     szmS.sre, szmS.srhgle, szmS.srlgle];
cerrSzmV = [szmS.gln, szmS.glnNorm, szmS.glv, szmS.hglze, szmS.lglze, szmS.lae, szmS.lahgle, ...
    szmS.lalgle, szmS.szn, szmS.sznNorm, szmS.szv, szmS.zp, ...
    szmS.sae, szmS.sahgle, szmS.salgle, szmS.ze];

pyradSzmNamC = {'GrayLevelNonUniformity', 'GrayLevelNonUniformityNormalized',...
    'GrayLevelVariance', 'HighGrayLevelZoneEmphasis',  'LowGrayLevelZoneEmphasis', ...
    'LargeAreaEmphasis', 'LargeAreaHighGrayLevelEmphasis', 'LargeAreaLowGrayLevelEmphasis',...
    'SizeZoneNonUniformity', 'SizeZoneNonUniformityNormalized', 'ZoneVariance', ...
    'ZonePercentage', 'SmallAreaEmphasis','SmallAreaHighGrayLevelEmphasis', ...
    'SmallAreaLowGrayLevelEmphasis','ZoneEntropy'};

pyradSzmNamC = strcat([pyradScanType,'_glszm_'],pyradSzmNamC);
pyRadSzmV = [];
for i = 1:length(pyradSzmNamC)
    if isfield(teststruct,pyradSzmNamC{i})
        pyRadSzmV(i) = teststruct.(pyradSzmNamC{i});
    else
        pyRadSzmV(i) = NaN;
    end
end
szmDiffV = (cerrSzmV - pyRadSzmV) ./ cerrSzmV * 100

%% NGLDM features

patchRadius3dV = paramS.higherOrderParamS.patchRadius3dV;
imgDiffThresh = paramS.higherOrderParamS.imgDiffThresh;
% 3d
ngldmM = calcNGLDM(quantizedM, patchRadius3dV, ...
    numGrLevels, imgDiffThresh);
featureS.ngldmFeatures3dS = ngldmToScalarFeatures(ngldmM,numVoxels);

ngldmS = featureS.ngldmFeatures3dS;
cerrNgldmV = [ngldmS.lde, ngldmS.hde, ngldmS.lgce, ngldmS.hgce, ...
    ngldmS.ldlge, ngldmS.ldhge, ngldmS.hdlge, ngldmS.hdhge, ...
    ngldmS.gln, ngldmS.glnNorm, ngldmS.dcn, ngldmS.dcnNorm,...
    ngldmS.dcp, ngldmS.glv, ngldmS.dcv, ngldmS.entropy, ngldmS.energy];

pyradNgldmNamC = {'SmallDependenceEmphasis', 'LargeDependenceEmphasis',...
    'LowGrayLevelCountEmphasis', 'HighGrayLevelCountEmphasis',  'SmallDependenceLowGrayLevelEmphasis', ...
    'SmallDependenceHighGrayLevelEmphasis', 'LargeDependenceLowGrayLevelEmphasis', ...
    'LargeDependenceHighGrayLevelEmphasis', 'GrayLevelNonUniformity', 'GrayLevelNonUniformityNorm', ...
    'DependenceNonUniformity', 'DependenceNonUniformityNormalized', ...
    'DependencePercentage', 'GrayLevelVariance', 'DependenceVariance', ...
    'DependenceEntropy', 'DependenceEnergy'};

pyradNgldmNamC = strcat([pyradScanType,'_gldm_'],pyradNgldmNamC);
pyRadNgldmV = [];
for i = 1:length(pyradNgldmNamC)
    if isfield(teststruct,pyradNgldmNamC{i})
        pyRadNgldmV(i) = teststruct.(pyradNgldmNamC{i});
    else
        pyRadNgldmV(i) = NaN;
    end
end

ngldmDiffV = (cerrNgldmV - pyRadNgldmV) ./ cerrNgldmV * 100

%% NGTDM features
numGrLevels = paramS.higherOrderParamS.numGrLevels;
patchRadius3dV = paramS.higherOrderParamS.patchRadius3dV;
% 3d
[s,p] = calcNGTDM(quantizedM, patchRadius3dV, numGrLevels);
featureS.ngtdmFeatures3dS = ngtdmToScalarFeatures(s,p,numVoxels);

ngtdmS = featureS.ngtdmFeatures3dS;
cerrNgtdmV = [ngtdmS.busyness ngtdmS.coarseness ngtdmS.complexity ngtdmS.contrast ngtdmS.strength];

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% Local Labeling of DOFs
node = [1,0,0; 0,1,0; 0,0,0; 0,0,1];
elem = [1 2 3 4];
localEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];
edge = zeros(20,2);
edge([1 12 5 20 11 4],:) = localEdge;
elem2dof = [1 12 5 20 11 4];
figure(1); clf;
set(gcf,'Units','normal');
set(gcf,'Position',[0,0,0.6,0.4]);
subplot(1,2,1)
showmesh3(node,elem);
view(-14,12);
findnode3(node);
findedge3(node,localEdge);

%% Face and face to local DOF mapping
face1 = [2 3 4];
f1dof = [6 5 4];
face2 = [1 4 3];
f2dof = [6 2 3];
face3 = [1 2 4];
f3dof = [5 3 1];
face4 = [1 3 2];
f4dof = [4 1 2];

%%

repo: Aero-matlab


text:
%% UNO Rosenbrock
clc
%Objective
nlprob.objective = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
x0 = [0 1]';

opts = [];
opts.algorithm = nloptSolver('LN_PRAXIS');
opts.maxfeval = 1000;
opts.maxtime = 1;
opts.display = 2;
nlprob.options = opts;

[x,f,e,i] = nlopt(nlprob,x0)

%% UNO Rosenbrock [GRAD]
clc
%Objective
nlprob = [];
nlprob.objective = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
nlprob.gradient = @(x) mklJac(nlprob.objective,x);
x0 = [0 1]';

opts = [];
opts.algorithm = nloptSolver('LD_LBFGS');
opts.maxfeval = 100;
opts.maxtime = 1;
opts.display = 2;
nlprob.options = opts;

[x,f,e,i] = nlopt(nlprob,x0)

%% BOUNDED Rosenbrock
clc
%Objective
nlprob = [];
nlprob.objective = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
nlprob.lb = [-5;-5];
nlprob.ub = [1;0.8];
x0 = [0 0]';

opts = [];
opts.algorithm = nloptSolver('LN_PRAXIS');
opts.maxfeval = 1000;
opts.maxtime = 10;
opts.display = 2;
nlprob.options = opts;

[x,f,e,i] = nlopt(nlprob,x0)

%% BOUNDED GRAD Rosenbrock
clc
%Objective
nlprob = [];
nlprob.objective = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
nlprob.gradient = @(x) mklJac(nlprob.objective,x);
nlprob.lb = [-5;-5];
nlprob.ub = [1;0.8];
x0 = [0 0]';

opts = [];
opts.algorithm = nloptSolver('LD_LBFGS');
opts.maxfeval = 1000;
opts.maxtime = 10;
opts.display = 2;
nlprob.options = opts;

[x,f,e,i] = nlopt(nlprob,x0)

%% NLCON HS71
clc
%Objective
obj = @(x) x(1)*x(4)*sum(x(1:3)) + x(3);
%Linear Constraints
lb = ones(4,1);
ub = 5*ones(4,1);
%Nonlinear Constraints
nlcon = @(x) [prod(x); sum(x.^2)];

nlprob =[];
nlprob.objective = obj;
nlprob.nlcon = nlcon;
nlprob.nlrhs = [25;40];
nlprob.nle = [1;0];
nlprob.lb = lb;
nlprob.ub = ub;

opts = [];
opts.algorithm = nloptSolver('LN_COBYLA');
opts.maxfeval = 1000;
opts.maxtime = 2;
opts.display = 2;
nlprob.options = opts;

x0 = zeros(4,1);

[x,f,e,i] = nlopt(nlprob,x0,opts)

%% NLCON GRAD HS71
%Objective & Gradient
obj = @(x) x(1)*x(4)*sum(x(1:3)) + x(3);
grad = @(x) [ x(1)*x(4) + x(4)*sum(x(1:3));
              x(1)*x(4);
              x(1)*x(4) + 1;
              x(1)*sum(x(1:3)) ];
%Linear Constraints
lb = ones(4,1);
ub = 5*ones(4,1);
%Nonlinear Constraints
nlcon = @(x) [prod(x); sum(x.^2)];
nljac = @(x) [prod(x)./x'; 2*x'];

nlprob =[];
nlprob.objective = obj;
nlprob.gradient = grad;
nlprob.nlcon = nlcon;
nlprob.nljac = nljac;
nlprob.nlrhs = [25;40];
nlprob.nle = [1;0];
nlprob.lb = lb;
nlprob.ub = ub;

opts = [];
opts.algorithm = nloptSolver('LD_SLSQP');
opts.maxfeval = 1000;
opts.maxtime = 2;
opts.display = 2;
nlprob.options = opts;

x0 = [1 5 5 1]';

[x,f,e,i] = nlopt(nlprob,x0,opts)

%% NLP1 Hock & Schittkowski #71
clc
%Objective & Gradient
obj = @(x) x(1)*x(4)*sum(x(1:3)) + x(3);
grad = @(x) [ x(1)*x(4) + x(4)*sum(x(1:3));
              x(1)*x(4);
              x(1)*x(4) + 1;
              x(1)*sum(x(1:3)) ];
%Linear Constraints
lb = ones(4,1);
ub = 5*ones(4,1);
%Nonlinear Constraints
nlcon = @(x) [ prod(x);
               sum(x.^2)];
nljac = @(x) [ prod(x)./x';
                2*x' ];
nlrhs = [25 40]';
nle = [1 0]'; % (>=, ==)
%Setup Options
opts = optiset('solver','nlopt','warnings','on','display','iter','solverOpts',nloptset('algorithm','LD_SLSQP'));
%Build & Solve
Opt = opti('obj',obj,'grad',grad,'nlmix',nlcon,nlrhs,nle,'nljac',nljac,'bounds',lb,ub,'options',opts)
x0 = [1 5 5 1]';
[x,fval,exitflag,info]= solve(Opt,x0)

%% LP [-31.4]
clc
% clear all
%Objective & Constraints
f = -[6 5]';
A = ([1,4; 6,4; 2, -5]);
b = [16;28;6];

nlprob = [];
nlprob.objective = @(x) f'*x;
nlprob.gradient = @(x) f;

nlprob.nlcon = @(x) A*x;
nlprob.nljac = @(x) A;
nlprob.nlrhs = b;
nlprob.nle = [-1;-1;-1];

nlprob.lb = [0;0];
nlprob.ub = [10;10];

opts = [];
opts.algorithm = nloptSolver('LD_SLSQP');
opts.maxfeval = 100;
opts.maxtime = 2;
opts.display = 2;
nlprob.options = opts;

x0 = [0;0];

[x,f,e,i] = nlopt(nlprob,x0)

%% LP1
clc
%Objective & Constraints
f = -[6 5]';
A = ([1,4; 6,4; 2, -5]);
b = [16;28;6];
%Build Object
Opt = opti('obj',@(x) f'*x,'grad',@(x) f,'ineq',A,b,'bounds',[0;0],[10;10],'options',optiset('solver','nlopt','solverOpts',nloptset('algorithm','LD_SLSQP')))
%Build & Solve
[x,fval,exitflag,info] = solve(Opt,[0;0])
%Plot
plot(Opt)
%Check Solution
[ok,msg] = checkSol(Opt)

%% Auglag NLP
clc
fun = @(x) log(1+x(1)^2) - x(2);
grad = @(x)[[(2*x(1))/(x(1)^2+1)];[-1]];
nlcon = @(x) (1 + x(1)^2)^2 + x(2)^2 - 4;
nljac = @(x)[[4*x(1)*(x(1)^2+1),2*x(2)]];
nlrhs = 0;
nle = 0;
x0 = [2;2];

nlprob = [];
nlprob.objective = fun;
nlprob.gradient = grad;

nlprob.nlcon = nlcon;
nlprob.nljac = nljac;
nlprob.nlrhs = nlrhs;
nlprob.nle = nle;

opts = [];
opts.algorithm = nloptSolver('LD_AUGLAG');
opts.maxfeval = 120;
opts.maxtime = 2;
opts.display = 2;

opts.local_optimizer.algorithm = nloptSolver('LD_LBFGS');
nlprob.options = opts;

[x,f,e,i] = nlopt(nlprob,x0)

%% Error checking
clc
clear
nlprob.objective = @(x) -1;
nlprob.lb = [-1;-1];
nlprob.ub = [1;1];

nlprob.nlcon = @(x) [-1;-1];
nlprob.nlrhs = [-1;-1];
nlprob.nle = [1;1];

nlprob.options.algorithm = 25;
x0 = [1;1];

[x,f,e,i] = nlopt(nlprob,x0)

%% Stop Val Checking
clc
%Objective
nlprob = [];
nlprob.objective = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
nlprob.lb = [-5;-5];
nlprob.ub = [1;1];
x0 = [0 0]';

opts = [];
opts.algorithm = nloptSolver('LN_PRAXIS');
opts.maxfeval = 1000;
opts.maxtime = 10;
opts.display = 2;
opts.stopval = 1e-6;
nlprob.options = opts;

[x,f,e,i] = nlopt(nlprob,x0)

%% Vector Storage Checking
clc
%Objective
nlprob = [];
nlprob.objective = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
nlprob.gradient = @(x) mklJac(nlprob.objective,x);
nlprob.lb = [-5;-5];
nlprob.ub = [1;0.8];
x0 = [0 0]';

opts = [];
opts.algorithm = nloptSolver('LD_LBFGS');
opts.maxfeval = 1000;
opts.maxtime = 10;
opts.display = 2;
opts.vector_storage = 1;
nlprob.options = opts;

[x,f,e,i] = nlopt(nlprob,x0)

%% Initial Pop Testing
clc
%Objective
nlprob = [];
nlprob.objective = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
nlprob.lb = [0;0];
nlprob.ub = [1;1];
x0 = [0 0]';

opts = [];
opts.algorithm = nloptSolver('GN_CRS2_LM');
opts.maxfeval = 1000;
opts.maxtime = 10;
opts.display = 1;
opts.initial_pop = 30;
nlprob.options = opts;

[x,f,e,i] = nlopt(nlprob,x0)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%% DEMO_additional_colormaps
% Below is a demonstration for:
%
% * Additional colormaps available in GIBBON

%%
clear; close all; clc;

%%
%Plot settings
fontSize=15;

% Create example data for visualizations
n=250;
s=1;
[X,Y]=ndgrid(linspace(-3*s,3*s,n));
Z=exp( -0.5.*((X./s).^2+(Y./s).^2));
Z=Z./max(Z(:));
Z(X<0)=-Z(X<0);
colorLim=[-1 1];
colormapset={'gjet','wjet','graygjet','grayjet',...
    'fire','ice','fireice','warmcold','iwarmcold','blood','bloodbone',...
    'che','wcbp','viridis','magma','inferno','turbo',...
    'plasma','gviridis','igviridis','kvw','cividis','spectral','ukr','flare'};

%% The gibbon color maps

%%
% Visualizing all additional color maps

n1=floor(sqrt(numel(colormapset)));
n2=ceil(numel(colormapset)/n1);

cFigure;

for q=1:1:numel(colormapset)

    h=subplot(n1,n2,q); hold on;
    title(colormapset{q},'FontSize',fontSize);
    imagesc(Z);
    colormap(h,colormapset{q}); colorbar;
    axis tight; axis equal; axis xy; box on;
    axis off;
    set(h,'FontSize',fontSize);
    clim(colorLim);
end
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [rotmat,cornerpoints,volume,surface,edgelength] = minboundbox(x,y,z,metric,level)
% minboundbox: Compute the minimal bounding box of points in 3d
% usage: [rotmat,cornerpoints,volume,surface,edgelength] = minboundbox(x,y,z,metric,level)
%
% arguments: (input)
%  x,y,z - vectors of points, describing points in 3d as
%        (x,y,z) pairs. x, y and z must be the same lengths.
%
%  metric - (OPTIONAL) - single letter character flag which
%        denotes the use of minimal volume, surface or sum of edges as the
%        metric to be minimized. metric may be either 'v', 's' or 'e',
%        capitalization is ignored.
%
%        DEFAULT: 'v'    ('volume')
%
%  level - (OPTIONAL) - either 1, 2, 3 or 4. This denotes the level of
%       reliability of the resulting minimal bounding box.
%       '1' denotes the search for a bounding box, where one side of the
%           box coincides with a face of the convex hull. (fast)
%       '2' like '1', but also incorporates a search through each pair of edges
%           which form a plane which coincides with one side of the box (slow)
%       '3' like '1', but also incorporates a search through each edge which is
%           parallel to an edge of the box (reasonably fast)
%       '4' like '1', '2' and '3' together. (slowest) (Never needed that.)
%
%       It depends on the application, what should be chosen here.
%       See the example at the end for the effects of this parameter.
%
%        DEFAULT: '3'
%
% arguments: (output)
%  rotmat - (3x3) rotation matrix for mapping of the pointcloud into a
%                 box which is axis-parallel (use inv(rotmat) for inverse
%                 mapping).
%
%  cornerpoints - (8x3) the cornerpoints of the bounding box.
%
%  volume - (scalar) volume of the minimal box itself.
%
%  surface - (scalar) surface of the minimal box as found.
%
%  edgelength - (scalar) sum of the edgelengths of the minimal box as found.
%
% Thanks to John d'Errico for providing the solution to the 2d version of
% this problem via minboundrect from the FEX. I took a huge amount of his
% syntax, organisation and comments for this submission. (Hope he wont sue
% me, though.)
%
% Thanks to Roger Stafford and John d'Errico for helpful discussions and
% suggestions.
%
% The algorithm is still not proven to produce the smallest possible
% enclosing box. But i never found any counterexamples for the 'best'
% algorithm with level 4, whilst the second example below shows, there are
% differences between the levels. I would be happy to have a proof that
% level 3 suffices for each case.
%
% Example usage:
%
% (1)
%      x = rand(100,1);
%      y = rand(100,1);
%      z = rand(100,1);
%      [rotmat,cornerpoints,volume,surface] = minboundbox(x,y,z);
%      plot3(x,y,z,'b.');hold on;plotminbox(cornerpoints,'r');
%
% (2)
%      x=[1,0,0.1,1];y=[0,1,0.1,1];z=[0,0,0.9,1];
%      [nerd,cornerpoints1,nerd,nerd] = minboundbox(x,y,z,'v',1);
%      [nerd,cornerpoints2,nerd,nerd] = minboundbox(x,y,z,'v',2);
%      [nerd,cornerpoints3,nerd,nerd] = minboundbox(x,y,z,'v',3);
%      plot3(x,y,z,'bo','LineWidth',5);hold on;
%      plotminbox(cornerpoints1,'b');hold on;
%      plotminbox(cornerpoints2,'m');hold on;
%      plotminbox(cornerpoints3,'r');axis equal;grid on;hold off;
%
%
% See also: minboundcircle, minboundtri, minboundsphere, minboundrect
%
%
% Author: Johannes Korsawe
% E-mail: johannes.korsawe@volkswagen.de
% Release: 1.0
% Release date: 09/01/2008

% default for metric
if (nargin<4) || isempty(metric)
  metric = 'v';
elseif ~ischar(metric)
  error 'metric must be a character flag if it is supplied.'
else
  % check for 'v', 's' or 'e'
  metric = lower(metric(:)');
  ind = strmatch(metric,{'volume','surface','edges'});
  if isempty(ind)
    error 'metric does not match either ''volume'', ''surface'' or ''edges''.'
  end
  % just keep the first letter.
  metric = metric(1);
end

% default for level
if (nargin<5 || isempty(level)),
    level=4;
elseif ~isnumeric(level) || (level~=1 && level~=2 && level~=3 && level~=4),
        error 'metric does not match either ''1'', ''2'', ''3'' or ''3''.'
end

% preprocess data
x = x(:);
y = y(:);
z = z(:);

% not many error checks to worry about
n1 = length(x);n2 = length(y);n3 = length(z);
if n1 ~= n2 || n1 ~= n3 || n2 ~= n3,
  error 'x, y and z must be the same sizes'
end

% start out with the convex hull of the points to
% reduce the problem dramatically. Note that any
% points in the interior of the convex hull are
% never needed, so we drop them.

try
 K = convhulln([x,y,z],{'Qt'});  % 'Pp' will silence the warnings

  % exclude those points inside the hull as not relevant
  % also sorts the points into their convex hull as a
  % closed polygon

 Ki = unique(K(:));
 [tf,loc] = ismember(K(:),Ki);K = reshape(loc,size(K));
 x = x(Ki);
 y = y(Ki);
 z = z(Ki);
 n1 =length(x);

catch
    error 'The number and/or distribution of given points does not allow the construction of a convex hull.'
end

% now we must find the bounding box of those points that remain

% get the angle of each face of the hull polygon.
% calculate local coordinate system
% x,y and z values of faces' cornerpoints
fx = [x(K(:,1)),x(K(:,2)),x(K(:,3))];
fy = [y(K(:,1)),y(K(:,2)),y(K(:,3))];
fz = [z(K(:,1)),z(K(:,2)),z(K(:,3))];
% faces' edges
v1 = [fx(:,2)-fx(:,1),fy(:,2)-fy(:,1),fz(:,2)-fz(:,1)];
v1 = v1./(repmat(sqrt(sum(v1.^2,2)),1,3));
v2 = [fx(:,3)-fx(:,1),fy(:,3)-fy(:,1),fz(:,3)-fz(:,1)];
v2 = v2-repmat(dot(v1,v2,2),1,3).*v1;
v2 = v2./(repmat(sqrt(sum(v2.^2,2)),1,3));
% faces' normals
nv = cross(v1,v2,2);
% faces' orientations measured in Euler angles
[alpha,beta,gamma]=euler123(v1,v2,nv);

nang = size(alpha,1);
d = inf;        % this will be the minimal value
rotmat=[];minmax=[];
xyz = [x,y,z];
if level==1 || level==3,    % check each face of the hull
    for i = 1:nang,
      % check current orientation
      % receive minimal value, rotation matrix and dimensions of the enclosing box
      [d, rotmat, minmax] = checkbox(alpha(i),beta(i),gamma(i),xyz,metric,d,rotmat,minmax);
    end
end

if level>1,
    e = edgelist(K);ne = size(e,1); % get a sorted list of edges
end

if level==2 || level==4,    % check each set of two edges of the convex hull (which obviously includes case '1' above)
    for i = 1:ne,   % go through edge list
        va = xyz(e(i,1),:)-xyz(e(i,2),:);
        va = va/norm(va);
        for j = i+1:ne, % go through remaining edge list
            vb = xyz(e(j,1),:)-xyz(e(j,2),:);
            vb = vb - dot(va,vb)*va; % orthogonaylize second edge wrt first
            nv = cross(va,vb);  % normal to plane formed by the two edges
            if sum(abs(nv))>0,
                vb = vb/norm(vb);nv = nv/norm(nv);
                [alp,bet,gam]=euler123(va,vb,nv);
                [d, rotmat, minmax] = checkbox(alp,bet,gam,xyz,metric,d,rotmat,minmax);
            end
        end
    end
end

if level==3 || level==4,    % check edge as parallel to one edge of the bounding box
    for i = 1:ne, % go through edge list
        % calculate rhs with edge as one of the axes
        va = xyz(e(i,1),:)-xyz(e(i,2),:);
        vb = [va(2),-va(1),0];if sum(abs(vb))==0,vb = [va(3),0,-va(1)];end
        va = va/norm(va);vb = vb/norm(vb);
        nv = cross(va,vb);
        % check all combinations of possible rhs
        [alp,bet,gam]=euler123(va,vb,nv);
        [d, rotmat, minmax] = checkbox(alp,bet,gam,xyz,metric,d,rotmat,minmax);
        [alp,bet,gam]=euler123(vb,nv,va);
        [d, rotmat, minmax] = checkbox(alp,bet,gam,xyz,metric,d,rotmat,minmax);
        [alp,bet,gam]=euler123(nv,va,vb);
        [d, rotmat, minmax] = checkbox(alp,bet,gam,xyz,metric,d,rotmat,minmax);
    end
end

% get the output values
h = [minmax(1,1),minmax(1,2),minmax(1,3);    % xmin,ymin,zmin
     minmax(2,1),minmax(1,2),minmax(1,3);    % xmax,ymin,zmin
     minmax(2,1),minmax(2,2),minmax(1,3);    % xmax,ymax,zmin
     minmax(1,1),minmax(2,2),minmax(1,3);    % xmin,ymax,zmin
     minmax(1,1),minmax(1,2),minmax(2,3);    % xmin,ymin,zmax
     minmax(2,1),minmax(1,2),minmax(2,3);    % xmax,ymin,zmax
     minmax(2,1),minmax(2,2),minmax(2,3);    % xmax,ymax,zmax
     minmax(1,1),minmax(2,2),minmax(2,3);    % xmin,ymax,zmax
     ];

cornerpoints = h*inv(rotmat); % minimal boxes' cornerpoints
h = minmax(2,:)-minmax(1,:);
volume = h(1)*h(2)*h(3);
surface = 2*(h(1)*h(2)+h(2)*h(3)+h(3)*h(1));
edgelength = 4*sum(h);

% all done

end % mainline end

function [alpha,beta,gamma]=euler123(v1,v2,nv)
% calculate Euler angles for the x, y', z'' Euler sequence
% see also
% http://en.wikipedia.org/wiki/Euler_angles
% http://de.wikipedia.org/wiki/Eulersche_Winkel

% --> beta
beta = asin(sign(nv(:,1)).*min(1,abs(nv(:,1))));
% --> alpha
alpha=0*beta;
i1 = find(nv(:,1)==1);i2 = setdiff(1:size(nv,1),i1);
if ~isempty(i1),
    alpha(i1) = asin(sign(v2(i1,3)).*min(1,abs(v2(i1,3))));
end
if ~isempty(i2),
    alpha(i2) = acos(sign(nv(i2,3)./cos(beta(i2))).*min(1,abs(nv(i2,3)./cos(beta(i2)))));
    i3 = find(sign(nv(i2,2)) ~= sign(-sin(alpha(i2)).*cos(beta(i2))));
    alpha(i2(i3)) = -alpha(i2(i3));
end
% --> gamma
gamma = 0*alpha;
if ~isempty(i2),
    singamma = -v2(i2,1)./cos(beta(i2));
    i21 = find(v1(i2,1)>=0);i22=setdiff(1:length(i2),i21);
    gamma(i2(i21)) = asin(sign(singamma(i21)).*min(1,abs(singamma(i21))));
    gamma(i2(i22)) = -pi-asin(sign(singamma(i22)).*min(1,abs(singamma(i22))));
end

end % function euler123

function [d, rotmat, minmax] = checkbox(alpha,beta,gamma,xyz,metric,d,rotmat,minmax)

% will need a 3x3 rotation matrix through (Euler X, Y', Z'')-angles alpha, beta and gamma
Rmat = @(alpha,beta,gamma) [cos(beta)*cos(gamma) -cos(beta)*sin(gamma) sin(beta)
                            sin(alpha)*sin(beta)*cos(gamma)+cos(alpha)*sin(gamma) -sin(alpha)*sin(beta)*sin(gamma)+cos(alpha)*cos(gamma) -sin(alpha)*cos(beta)
                            -cos(alpha)*sin(beta)*cos(gamma)+sin(alpha)*sin(gamma) cos(alpha)*sin(beta)*sin(gamma)+sin(alpha)*cos(gamma) cos(alpha)*cos(beta)];

rot = Rmat(alpha,beta,gamma);
xyz_i = xyz*rot;                      % now the actual face is in the x-y plane

x_i = xyz_i(:,1);y_i = xyz_i(:,2);          % .. so take only the x and y values
rot2 = minrect(x_i,y_i,metric);             % find the optimal rotation around z-axis
rot = rot*[[rot2,[0;0]];[0,0,1]];           % combine that with the formar rotation
xyz_i = xyz*rot;                            % now again the xyz_i values, but in optimal axisparallel shape

xyzmin = min(xyz_i,[],1);
xyzmax = max(xyz_i,[],1);
h = xyzmax-xyzmin;

if metric == 'v',       % smallest volume
    d_i = h(1)*h(2)*h(3);
elseif metric == 's',   % smallest surface
    d_i = h(1)*h(2)+h(2)*h(3)+h(3)*h(1);
else,                   % smallest sum of edges
    d_i = sum(h);
end

if d_i < d,
    d = d_i;
    rotmat = rot;
    minmax = [xyzmin;xyzmax];
end

end % function checkbox

function rot = minrect(x,y,metric)
% see comments from minboundrect of John d'Errico
% i am only interested in the additional rotation matrix around [0,0,1]
edges = convhull(x,y,{'Qt'});
x = x(edges);y = y(edges);
ind = 1:length(x)-1;
Rmat = @(theta) [cos(theta) sin(theta);-sin(theta) cos(theta)];
edgeangles = atan2(y(ind+1) - y(ind),x(ind+1) - x(ind));edgeangles = unique(mod(edgeangles,pi/2));
nang = length(edgeangles);
area = inf;perimeter = inf;met = inf;xy = [x,y];
for i = 1:nang
  rot_i = Rmat(-edgeangles(i));xyr = xy*rot_i;
  xymin = min(xyr,[],1);xymax = max(xyr,[],1);
  A_i = prod(xymax - xymin);P_i = 2*sum(xymax-xymin);
  if metric=='v', M_i = A_i;else, M_i = P_i;end
  if M_i<met, rot = rot_i;met = M_i;end
end

end % function minrect

function e = edgelist(K)
% generate an edgelist from the convex body K
e = [K(:,1),K(:,2)];e = [e;K(:,2),K(:,3)];e = [e;K(:,1),K(:,3)];
e = [min(e,[],2),max(e,[],2)];ne = size(e,1);
e = [e,e(:,2)+ne*e(:,1)];
[nerd,I] = sort(e(:,3));
e = e(I(1:2:end),:);

end % edgelist


text:
%% 3D Model Demo
% This is short demo that loads and renders a 3D model of a human femur. It
% showcases some of MATLAB's advanced graphics features, including lighting and
% specular reflectance.

% Copyright 2011 The MathWorks, Inc.

%% Load STL mesh
% Stereolithography (STL) files are a common format for storing mesh data. STL
% meshes are simply a collection of triangular faces. This type of model is very
% suitable for use with MATLAB's PATCH graphics object.

% Import an STL mesh, returning a PATCH-compatible face-vertex structure
fv = stlread('femur.stl');

%% Render
% The model is rendered with a PATCH graphics object. We also add some dynamic
% lighting, and adjust the material properties to change the specular
% highlighting.

patch(fv,'FaceColor',       [0.8 0.8 1.0], ...
         'EdgeColor',       'none',        ...
         'FaceLighting',    'gouraud',     ...
         'AmbientStrength', 0.15);

% Add a camera light, and tone down the specular highlighting
camlight('headlight');
material('dull');

% Fix the axes scaling, and set a nice view angle
axis('image');
view([-135 35]);


text:
function camBfloat2PDB(bfloat_filename,pdb_filename,ctr_filename,fa_filename)

opt_ctr = [' -i ' ctr_filename];
opt_bfloat = [' ' bfloat_filename];
opt_pdb = [' -p ' pdb_filename];

if ~exist(ctr_filename,'file') && exist(fa_filename,'file')
    % Create contrack parameters script if it doesn't exist
    % Bfloat conversion currently uses contrack so we need ctr file
    ctr = ctrCreate();
    % Must set pdb image, but can be anything
    % XXX
    % Must be absolute pathway or something is going wrong with relative
    % naming
    % XXX
    [pathstr, name, ext] = fileparts(fa_filename);
    ctr.pdf_filename = [name ext];
    ctr.image_directory = [pathstr '/'];
    ctrSave(ctr,ctr_filename);

end

cmd = ['contrack_score.glxa64' opt_ctr opt_pdb '  --thresh 100000 --seq --bfloat_no_stats' opt_bfloat];
display(cmd);
system(cmd,'-echo');

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef SVC1V1 < Algorithm
    %SVC1V1 Support Vector Classifier using one-vs-one approach
    %classification by predicting class labels as a regression problem.
    %It uses libSVM-weight SVM implementation.
    %
    %   SVC1V1 methods:
    %      fitpredict               - runs the corresponding algorithm,
    %                                   fitting the model and testing it in a dataset.
    %      fit                        - Fits a model from training data
    %      predict                    - Performs label prediction
    %
    %   References:
    %     [1] C.-W. Hsu and C.-J. Lin
    %         A comparison of methods for multi-class support vector machines
    %         IEEE Transaction on Neural Networks,vol. 13, no. 2, pp. 415–425, 2002.
    %         https://doi.org/10.1109/72.991427
    %     [2] P.A. Gutiérrez, M. Pérez-Ortiz, J. Sánchez-Monedero,
    %         F. Fernández-Navarro and C. Hervás-Martínez
    %         Ordinal regression methods: survey and experimental study
    %         IEEE Transactions on Knowledge and Data Engineering, Vol. 28. Issue 1
    %         2016
    %         http://dx.doi.org/10.1109/TKDE.2015.2457911
    %     [3] LibSVM website: https://www.csie.ntu.edu.tw/~cjlin/libsvm
    %
    %   This file is part of ORCA: https://github.com/ayrna/orca
    %   Original authors: Pedro Antonio Gutiérrez, María Pérez Ortiz, Javier Sánchez Monedero
    %   Citation: If you use this code, please cite the associated paper http://www.uco.es/grupos/ayrna/orreview
    %   Copyright:
    %       This software is released under the The GNU General Public License v3.0 licence
    %       available at http://www.gnu.org/licenses/gpl-3.0.html
    properties
        description = 'Support Vector Machine Classifier with 1vs1 paradigm';
        parameters = struct('C', 0.1, 'k', 0.1);
    end
    properties (Access = private)
        algorithmMexPath = fullfile(fileparts(which('Algorithm.m')),'libsvm-weights-3.12','matlab');
    end

    methods

        function obj = SVC1V1(varargin)
            %SVC1V1 constructs an object of the class SVC1V1 and sets its default
            %   characteristics
            %   OBJ = SVC1V1(KERNEL) builds SVC1V1 with RBF as kernel function
            obj.parseArgs(varargin);
        end

        function [projectedTrain, predictedTrain]= privfit( obj, train , param)
            %PRIVFIT trains the model for the SVC1V1 method with TRAIN data and
            %vector of parameters PARAMETERS.
            if isempty(strfind(path,obj.algorithmMexPath))
                addpath(obj.algorithmMexPath);
            end

            weights = ones(size(train.targets));
            options = ['-t 2 -c ' num2str(param.C) ' -g ' num2str(param.k) ' -q'];
            model.libsvmModel = svmtrain(weights, train.targets, train.patterns, options);
            model.parameters = param;
            obj.model = model;
            [predictedTrain, acc, projectedTrain] = svmpredict(train.targets,train.patterns,model.libsvmModel, '');
            if ~isempty(strfind(path,obj.algorithmMexPath))
                rmpath(obj.algorithmMexPath);
            end
        end

        function [projected, predicted]= privpredict(obj,test)
            %PREDICT predicts labels of TEST patterns labels. The object needs to be fitted to the data first.
            if isempty(strfind(path,obj.algorithmMexPath))
                addpath(obj.algorithmMexPath);
            end
            [predicted, acc, projected] = svmpredict(ones(size(test,1),1),test,obj.model.libsvmModel, '');
            if ~isempty(strfind(path,obj.algorithmMexPath))
                rmpath(obj.algorithmMexPath);
            end

        end
    end
end

author: ayrna

repo: orca

save_path: github-repos/MATLAB/ayrna-orca


text:
classdef EKF_filter_jekeli76 < handle
    properties (Hidden)
        type = 'ekf';
        tag  = 'EKF_crude';  % ID tag
        covDim=3; % the dimension of covariance matrix, Rx RY RZ in earth centered
        % n-frame, vn, RPY, acc bias, gyro bias, acc scale and gyro scale
    end
    % The following properties can be set only by class methods
    properties (SetAccess = private)
        dt; % sampling interval of IMU, unit sec
        % the translation from antenna to mems frame, i.e., the antenna's
        % position in the mems frame.
        xi;
        vi;
        bias=0;
        p_k_k; % the covariance of the entire state vector
        psdq; % process noise PSD
    end
    methods
        function filter = EKF_filter_jekeli76(initX, initV, deltat, procnoiseq)
            filter.xi=initX;
            filter.vi=initV;
            filter.dt=deltat;
            filter.psdq=procnoiseq;
            filter.p_k_k=eye(filter.covDim);
        end
        %===============================================================================================
        %-- State transition function
        % propogate state with accelerometer and gyro input at time k-1 to state at k,
        % i.e., X(k|k-1), from state at k-1. U1 contains IMU measurement
        % U1 is the measured acceleration
        function ffun_state(filter, U1)

            ahat=U1-filter.bias;
            vi1=filter.vi+ahat*filter.dt;
           filter.xi=filter.xi+(vi1+filter.vi)/2*filter.dt;
           filter.vi=vi1;
        end
        function ffun_covariance(filter)
            %propagate the covariance corresponds to states, rs in e, vs in e, q s2e,
            deltat=filter.dt;
            STM=[1, 0, deltat; deltat, 1, .5*deltat^2; 0,0, 1];
            Qd=filter.psdq*[deltat, .5*deltat^2, 0; .5*deltat^2, deltat^3/3, 0; 0, 0,0];
            filter.p_k_k=STM*filter.p_k_k*STM'+Qd;  % the covariance of the navigation states and imu error terms
        end
        %==============================================================================================
        %update both the covariance and state
        function correctstates(filter, predict,measure, H,R)
            p_km1_k=filter.p_k_k;
            inno= predict-measure;
            %Kalman
            K=p_km1_k*H'/(H*p_km1_k*H'+R);
            deltaX=K*inno;
            % update covariance
            filter.p_k_k=(eye(filter.covDim)-K*H)*p_km1_k*(eye(filter.covDim)-K*H)'+K*R*K';
            % compute updated state
            filter.xi=filter.xi-deltaX(1);
            filter.vi=filter.vi-deltaX(2);
            filter.bias = filter.bias + deltaX(3);
        end
        function SaveToFile(filter, preimutime, ffilres)
            %Write result to the files
            fprintf(ffilres,'%6.2f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\t%6.3f\n', ...
                [preimutime;filter.vi;filter.xi;filter.bias;sqrt(diag(filter.p_k_k))]);
        end
    end
end

repo: ekfmonoslam


text:
function [ x, seed ] = normal_01_sample ( seed )

%*****************************************************************************80
%
%% NORMAL_01_SAMPLE samples the standard normal probability distribution.
%
%  Discussion:
%
%    The standard normal probability distribution function (PDF) has
%    mean 0 and standard deviation 1.
%
%    The Box-Muller method is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, X, a sample of the standard normal PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ r1, seed ] = r8_uniform_01 ( seed );
  [ r2, seed ] = r8_uniform_01 ( seed );
  x = sqrt ( -2.0 * log ( r1 ) ) * cos ( 2.0 * pi * r2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = drawBezierCurve(varargin)
%DRAWBEZIERCURVE Draw a cubic bezier curve defined by 4 control points.
%
%   drawBezierCurve(POINTS)
%   Draw the Bezier curve defined by the 4 control points stored in POINTS.
%   POINTS is either a 4-by-2 array (vertical concatenation of control
%   points coordinates), or a 1-by-8 array (horizontal concatenation of
%   control point coordinates).
%
%   drawBezierCurve(..., PARAM, VALUE)
%   Specifies additional drawing parameters, see the line function for
%   details.
%
%   drawBezierCurve(AX, ...);
%   Spcifies the handle of the axis to draw on.
%
%   H = drawBezierCurve(...);
%   Return a handle to the created graphic object.
%
%
%   Example
%     drawBezierCurve([0 0;5 10;10 5;10 0]);
%     drawBezierCurve([0 0;5 10;10 5;10 0], 'linewidth', 2, 'color', 'g');
%
%   See also
%     drawPolyline, cubicBezierToPolyline
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2011-03-16, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

% extract handle of axis to draw on
if isAxisHandle(varargin{1})
    ax = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

points = varargin{1};
varargin(1) = [];

% default number of discretization steps
N = 64;

% check if discretization step is specified
if ~isempty(varargin)
    var = varargin{1};
    if length(var) == 1 && isnumeric(var)
        N = round(var);
        varargin(1) = [];
    end
end

% convert control coordinates to polyline
poly = cubicBezierToPolyline(points, N);

% draw the curve
h = drawPolyline(ax, poly, varargin{:});

% eventually return a handle to the created object
if nargout > 0
    varargout = {h};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [] = plotKspaceVolumes(kspaceVolumes, scales)
%PLOTKSPACEVOLUMES Plot the kspace 3D volumes.
%   volumes: Cell array of volumes to plot
%   scales: Cell array of scales ([min, max]) for the image plots.
%
%   Code refractored from Berkin Bilgic's scripts: "script_Laplacian_unwrap_Sharp_Fast_TV_gre3D.m"
%   and "script_Laplacian_unwrap_Sharp_Fast_TV_gre3D.m"
%   Original source: https://martinos.org/~berkin/software.html
%
%   Original reference:
%   Bilgic et al. (2014), Fast quantitative susceptibility mapping with
%   L1-regularization and automatic parameter selection. Magn. Reson. Med.,
%   72: 1444-1459. doi:10.1002/mrm.25029

    for volumeIndex = 1:length(kspaceVolumes)
        figure(), subplot(1,3,1), imagesc( kspaceVolumes{volumeIndex}(:,:,1+end/2), scales{volumeIndex} ), axis square off, colormap gray
        figure(get(gcf,'Number')), subplot(1,3,2), imagesc( squeeze(kspaceVolumes{volumeIndex}(:,1+end/2,:)), scales{volumeIndex} ), axis square off, colormap gray
        figure(get(gcf,'Number')), subplot(1,3,3), imagesc( squeeze(kspaceVolumes{volumeIndex}(1+end/2,:,:)), scales{volumeIndex} ), axis square off, colormap gray
    end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
classdef LaunchVehicleViewProfileVectorData < matlab.mixin.SetGet
    %LaunchVehicleViewProfileVectorData Summary of this class goes here
    %   Detailed explanation goes here

    properties
        %these are for vehicle trajectory data and not the point itself!
        timesArr(1,:) cell = {};
        xInterps(1,:) cell = {};
        yInterps(1,:) cell = {};
        zInterps(1,:) cell = {};
        vxInterps(1,:) cell = {};
        vyInterps(1,:) cell = {};
        vzInterps(1,:) cell = {};

        markerPlot = matlab.graphics.GraphicsPlaceholder.empty(1,0)
        vector AbstractGeometricVector
        viewFrame AbstractReferenceFrame
    end

    methods
        function obj = LaunchVehicleViewProfileVectorData(vector, viewFrame)
            obj.vector = vector;
            obj.viewFrame = viewFrame;
        end

        function addData(obj, times, rVects, vVects)
            obj.timesArr{end+1} = times;

            if(length(times) >= 3)
                method = 'pchip';
            else
                method = 'linear';
            end

            obj.xInterps{end+1} = griddedInterpolant(times, rVects(:,1), method, 'linear');
            obj.yInterps{end+1} = griddedInterpolant(times, rVects(:,2), method, 'linear');
            obj.zInterps{end+1} = griddedInterpolant(times, rVects(:,3), method, 'linear');
            obj.vxInterps{end+1} = griddedInterpolant(times, vVects(:,1), method, 'linear');
            obj.vyInterps{end+1} = griddedInterpolant(times, vVects(:,2), method, 'linear');
            obj.vzInterps{end+1} = griddedInterpolant(times, vVects(:,3), method, 'linear');
        end

        function plotVectorAtTime(obj, time, hAx)
            [origin, vect] = obj.getVectorAtTime(time);

            if(isempty(obj.markerPlot))
                hold(hAx,'on');
                for(i=1:size(vect,2))
                    obj.markerPlot(i) = quiver3(hAx, origin(1,i),origin(2,i),origin(3,i), vect(1,i),vect(2,i),vect(3,i), 'AutoScale','off', 'Color',obj.vector.lineColor.color, 'LineStyle',obj.vector.lineSpec.linespec);
                end
                hold(hAx,'off');
            else
                for(i=1:size(vect,2))
                    obj.markerPlot(i).XData = origin(1,i);
                    obj.markerPlot(i).YData = origin(2,i);
                    obj.markerPlot(i).ZData = origin(3,i);

                    obj.markerPlot(i).UData = vect(1,i);
                    obj.markerPlot(i).VData = vect(2,i);
                    obj.markerPlot(i).WData = vect(3,i);
                end
            end
        end
    end

    methods(Access=private)
        function [origin, vect] = getVectorAtTime(obj, time)
            origin = [];
            vect = [];
            for(i=1:length(obj.timesArr))
                times = obj.timesArr{i};

                if(time >= min(floor(times)) && time <= max(ceil(times)))
                    xInterp = obj.xInterps{i};
                    x = xInterp(time);

                    yInterp = obj.yInterps{i};
                    y = yInterp(time);

                    zInterp = obj.zInterps{i};
                    z = zInterp(time);

                    vxInterp = obj.vxInterps{i};
                    vx = vxInterp(time);

                    vyInterp = obj.vyInterps{i};
                    vy = vyInterp(time);

                    vzInterp = obj.vzInterps{i};
                    vz = vzInterp(time);

                    vehElemSet = CartesianElementSet(time, [x;y;z], [vx;vy;vz], obj.viewFrame);
                    vect(:,end+1) = obj.vector.getVectorAtTime(time, vehElemSet, obj.viewFrame); %#ok<AGROW>

                    origin(:,end+1) = obj.vector.getOriginPointInViewFrame(time, vehElemSet, obj.viewFrame); %#ok<AGROW>
                end
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [croppedPeakImages] = compareSignalToMovie(inputMovie, inputImages, inputSignal, varargin)
	% Shows a cropped version of inputMovie for each inputImages and aligns it to inputSignal peaks to make sure detection is working.
	% Biafra Ahanonu
	% started: 2013.11.04 [18:40:45]
	% inputs
		% inputMovie - matrix dims are [X Y t] - where t = number of time points
		% inputImages - matrix dims are [X Y n] - where n = number of filters, NOTE THE DIFFERENCE
		% inputSignal - matrix dims are [n t] - where n = number of signals, t = number of time points
	% outputs
		% none, this is a display function

	[croppedPeakImages] = ciapkg.view.compareSignalToMovie(inputMovie, inputImages, inputSignal,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function cnn_imagenet_minimal()
% CNN_IMAGENET_MINIMAL   Minimalistic demonstration of how to run an ImageNet CNN model

% setup toolbox
run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', '..', 'matlab', 'vl_setupnn.m')) ;

% download a pre-trained CNN from the web
if ~exist('imagenet-vgg-f.mat', 'file')
  fprintf('Downloading a model ... this may take a while\n') ;
  urlwrite('http://www.vlfeat.org/matconvnet/models/imagenet-vgg-f.mat', ...
    'imagenet-vgg-f.mat') ;
end
net = load('imagenet-vgg-f.mat') ;

% obtain and preprocess an image
im = imread('peppers.png') ;
im_ = single(im) ; % note: 255 range
im_ = imresize(im_, net.meta.normalization.imageSize(1:2)) ;
im_ = im_ - net.meta.normalization.averageImage ;

% run the CNN
res = vl_simplenn(net, im_) ;

% show the classification result
scores = squeeze(gather(res(end).x)) ;
[bestScore, best] = max(scores) ;
figure(1) ; clf ; imagesc(im) ;
title(sprintf('%s (%d), score %.3f',...
   net.meta.classes.description{best}, best, bestScore)) ;

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function rgb3d=label2rgb3d(varargin)
[label,map,zerocolor,order,fcnflag] = parse_inputs(varargin{:});
%label= 3d image with labels
% map= specified color map
%==================================
numregion = double(max(label(:)));

% If MAP is a function, evaluate it.  Make sure that the evaluated function
% returns a valid colormap.
if  fcnflag == 1
    cmap = feval(map, numregion);
    if ~isreal(cmap) || any(cmap(:) > 1) || any(cmap(:) < 0) || ...
            ~isequal(size(cmap,2),3) || size(cmap,1) < 1
        eid = sprintf('Images:%s:functionReturnsInvalidColormap',mfilename);
        msg = 'function handle MAP must return a n x 3 colormap array';
        error(eid,'%s',msg);
    end
else
    cmap = map;
end

% If ORDER is set to 'shuffle', save original state.  The SHUFFLE keyword
% uses the same state every time it is called. After shuffling is completed,
% go back to original state.
if isequal(order,'shuffle')
    S = rand('state');
    rand('state', 0);
    index = randperm(numregion);
    cmap = cmap(index,:,:);
    rand('state', S);
end

% Issue a warning if the zerocolor (boundary color) matches the color of one
% of the regions.
for i=1:numregion
  if isequal(zerocolor,cmap(i,:))
    wid= sprintf('Images:%s:zerocolorSameAsRegionColor',mfilename);
    msg= sprintf('Region number %d has the same color as the ZEROCOLOR.',i);
    warning(wid,'%s',msg);
  end
end
cmap = [zerocolor;cmap];

% if label is of type double, need to pass 'label + 1' into IND2RGB.
% IND2RGB does not like double arrays containing zero values.
if isa(label, 'double')
    rgb3d = ind2rgb3d(label + 1, cmap);
else
    rgb3d = ind2rgb3d(label, cmap);
end

%======================================================================
%======================================================================
function [L, Map, Zerocolor, Order, Fcnflag] = parse_inputs(varargin)
% L         label 3d matrix: matrix containing non-negative values.
% Map       colormap: name of standard colormap, user-defined map, function
%           handle.
% Zerocolor RGB triple or Colorspec
% Order     keyword if specified: 'shuffle' or 'noshuffle'
% Fcnflag   flag to indicating that Map is a function

valid_order = {'shuffle', 'noshuffle'};
iptchecknargin(1,4,nargin,mfilename);

% set defaults
L = varargin{1};
Map = 'jet';
Zerocolor = [1 1 1];
Order = 'noshuffle';
Fcnflag = 0;

% parse inputs
if nargin > 1
    Map = varargin{2};
end
if nargin > 2
    Zerocolor = varargin{3};
end
if nargin > 3
    Order = varargin{4};
end

% error checking for L
iptcheckinput(L,{'numeric' 'logical'}, ...
              {'real' 'nonsparse' 'finite' 'nonnegative' 'integer'}, ...
              mfilename,'L',1);

% error checking for Map
[fcn, fcnchk_msg] = fcnchk(Map);
if isempty(fcnchk_msg)
    Map = fcn;
    Fcnflag = 1;
else
    if isnumeric(Map)
        if ~isreal(Map) || any(Map(:) > 1) || any(Map(:) < 0) || ...
                    ~isequal(size(Map,2), 3) || size(Map,1) < 1
          eid = sprintf('Images:%s:invalidColormap',mfilename);
          msg = 'Invalid entry for MAP.';
          error(eid,'%s',msg);
        end
    else
        eid = sprintf('Images:%s:invalidFunctionforMAP',mfilename);
        error(eid,'%s',fcnchk_msg);
    end
end

% error checking for Zerocolor
if ~ischar(Zerocolor)
    % check if Zerocolor is a RGB triple
    if ~isreal(Zerocolor) || ~isequal(size(Zerocolor),[1 3]) || ...
                any(Zerocolor> 1) || any(Zerocolor < 0)
      eid = sprintf('Images:%s:invalidZerocolor',mfilename);
      msg = 'Invalid RGB triple entry for ZEROCOLOR.';
      error(eid,'%s',msg);
    end
else
    [cspec, msg] = cspecchk(Zerocolor);
    if ~isempty(msg)
        eid = sprintf('Images:%s:notInColorspec',mfilename);
        error(eid,'%s',msg);
    else
        Zerocolor = cspec;
    end
end

% error checking for Order
idx = strmatch(lower(Order), valid_order);
eid = sprintf('Images:%s:invalidEntryForOrder',mfilename);
if isempty(idx)
    msg = 'Valid entries for ORDER are ''shuffle'' or ''noshuffle''.';
    error(eid,'%s',msg);
elseif length(idx) > 1
    msg = sprintf('Ambiguous string for ORDER: %s.', Order);
    error(eid,'%s',msg);
else
    Order = valid_order{idx};
end

%================================================================
%=================================================================
function [rout,g,b] = ind2rgb3d(a,cm)
%IND2RGB Convert indexed image to RGB image.
%   RGB = IND2RGB(X,MAP) converts the 3d matrix X and corresponding
%   colormap MAP to RGB (truecolor) format.
%
%   Class Support
%   -------------
%   X can be of class uint8, uint16, or double. RGB is an
%   M-by-N-by-3 array of class double.
%
%   See also IND2GRAY, RGB2IND (in the Image Processing Toolbox).

if ~isa(a, 'double')
    a = double(a)+1;    % Switch to one based indexing
end

error(nargchk(2,2,nargin));

% Make sure A is in the range from 1 to size(cm,1)
a = max(1,min(a,size(cm,1)));

% Extract r,g,b components
r = zeros(size(a)); r(:) = cm(a,1);
g = zeros(size(a)); g(:) = cm(a,2);
b = zeros(size(a)); b(:) = cm(a,3);

if nargout==3,
  rout = r;
else
  rout = zeros([size(r),3]);
  rout(:,:,:,1) = r;
  rout(:,:,:,2) = g;
  rout(:,:,:,3) = b;
end