text:
function OutputFile = dipoles_merge(DipoleFiles)
% DIPOLES_MERGE: Merge multiple dipoles files

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Elizabeth Bock, Francois Tadel, Jeremy Moreau, 2015

% Load in all dipole files
for i = 1:length(DipoleFiles)
    % Check file path file
    if ~file_exist(DipoleFiles{i})
        DipoleFiles{i} = file_fullpath(DipoleFiles{i});
    end
    if ~file_exist(DipoleFiles{i})
        error(['File not found: ' DipoleFiles{i}]);
    end
    % Load file
    DipoleMat = load(DipoleFiles{i});
    % Do not accept subsets of dipoles
    if (DipoleMat.Subset > 1)
        error('TODO: Update this function for merging files with subsets.')
    end
    % First file: template
    if (i == 1)
        MergeMat = DipoleMat;
        MergeMat.DipoleNames = {};
        MergeMat.Dipole = [];

        % Calculate sampling rate of first dipole file
        if (length(DipoleMat.Time) == 1)
            DipOneSamplingRate = 0;  % Set sampling rate to 0 if only one dipole
        elseif (length(DipoleMat.Time) > 1)
            DipOneSamplingRate = DipoleMat.Time(2) - DipoleMat.Time(1);
        end

    % Following files: check that they are compatible
    else
        if ((DipOneSamplingRate ~= 0) && (length(DipoleMat.Time) > 1))
            if (abs(DipOneSamplingRate - (DipoleMat.Time(2) - DipoleMat.Time(1))) > 0.00001)
                error('Only files with equal sampling rate can be merged.');
            end
        elseif (DipoleMat.Subset ~= MergeMat.Subset)
            error('Only files with equal number of channel subsets can be merged.');
        end
    end

    % Names just have the group number (i.e. Group #1, Group #2)
    lastGroupNumber = length(MergeMat.DipoleNames);
    nGroups = length(DipoleMat.DipoleNames);
    % Loop through the groups
    for g = 1:nGroups
        % Add the name
        MergeMat.DipoleNames{end+1} = ['Group #' num2str(lastGroupNumber + g)];
        % Update the index numbers
        ind = find([DipoleMat.Dipole.Index] == g);
        if ~isempty(ind)
            [DipoleMat.Dipole(ind).Index] = deal(lastGroupNumber + g);
        end
    end
    % Merge the Dipole structures
    MergeMat.Dipole = [MergeMat.Dipole, DipoleMat.Dipole];
    MergeMat.DataFile = '';
    % Merge history
    MergeMat = bst_history('add', MergeMat, 'merge', ['Merged file: ' file_short(DipoleFiles{i})]);
end

% Set new Time sampling for merged dipoles
MergeMat.Time = unique([MergeMat.Dipole.Time]);

% ===== SAVE NEW FILE =====
% Update the comment to reflect the number of merged files
MergeMat.Comment = ['Merge: ' num2str(length(DipoleFiles)) ' files'];
% Create output filename
OutputFile = file_unique(bst_fullfile(fileparts(DipoleFiles{1}), 'dipoles_merged.mat'));
% Save new file in Brainstorm format
bst_save(OutputFile, MergeMat, 'v7');

% ===== UPDATE DATABASE =====
% Get study of the first file
[sStudy,iStudy] = bst_get('DipolesFile', file_short(DipoleFiles{1}));
% Create structure
BstDipolesMat = db_template('Dipoles');
BstDipolesMat.FileName = file_short(OutputFile);
BstDipolesMat.Comment  = MergeMat.Comment;
BstDipolesMat.DataFile = MergeMat.DataFile;
% Add to study
sStudy = bst_get('Study', iStudy);
iDipole = length(sStudy.Dipoles) + 1;
sStudy.Dipoles(iDipole) = BstDipolesMat;
% Save study
bst_set('Study', iStudy, sStudy);
% Update tree
panel_protocols('UpdateNode', 'Study', iStudy);
% Select node
panel_protocols('SelectNode', [], BstDipolesMat.FileName);
% Save database
db_save();

author: brainstorm-tools


text:
function springDampers = springDampersVisualise(springDamperVals, startPoints, ...
                                         staticPoints, widths, springConst, ...
                                         damperConst)

% SPRINGDAMPERSVISUALISE Helper code for showing an spring dampers during 2-D visualisation.
% FORMAT
% DESC is a helper function for plotting spring damper locations using latent
% variable models.
% ARG springDamperValues : the values to set the springDamper data to.
% ARG startPoints : the fixed points of the spring damers.
% ARG widths : the widths of the spring dampers.
% ARG springConst : the spring constants of the spring dampers.
% ARG damperConst : the damper constants of the spring dampers.
% RETURN springDampers : contains array of spring damper objects.
%
% COPYRIGHT : Neil D. Lawrence, 2007
%
% SEEALSO : springDampersModify, fgplvmResultsDynamic

% MLTOOLS

for i = 1:length(springDamperVals)
  springDampers(i) = springDamperCreate(startPoints(i, :), ...
                                        staticPoints(i, :), ...
                                        widths(i), springConst(i), ...
                                        damperConst(i));
  springDampers(i).selected = 0;
end
springDampers = springDampersModify(springDampers, springDamperVals, ...
                                       startPoints, staticPoints, widths, ...
                                       springConst, damperConst);
for i = 1:length(springDamperVals)
  set(springDampers(i).handle, 'linewidth', 2)
  set(springDampers(i).spring.handle, 'linewidth', 2)
  set(springDampers(i).damper.handle, 'linewidth', 2)
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function hf_out = lhs_operation_joint(hf, samplesf, reg_filter, feature_reg, init_samplef, XH, init_hf, proj_reg)

% This is the left-hand-side operation in Conjugate Gradient

hf_out = cell(size(hf));

% Extract projection matrix and filter separately
P = cellfun(@real, hf(2,1,:), 'uniformoutput',false);
hf = hf(1,1,:);

% size of the padding
num_features = length(hf);
output_sz = [size(hf{1},1), 2*size(hf{1},2)-1];

% Compute the operation corresponding to the data term in the optimization
% (blockwise matrix multiplications)
%implements: A' diag(sample_weights) A f

% sum over all features in each block
sh_cell = cell(1,1,num_features);
for k = 1:num_features
    sh_cell{k} = mtimesx(samplesf{k}, permute(hf{k}, [3 4 1 2]), 'speed');
end

% sum over all feature blocks
sh = sh_cell{1};    % assumes the feature with the highest resolution is first
pad_sz = cell(1,1,num_features);
for k = 2:num_features
    pad_sz{k} = (output_sz - [size(hf{k},1), 2*size(hf{k},2)-1]) / 2;

    sh(:,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) = ...
        sh(:,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) + sh_cell{k};
end

% weight all the samples
% sh = bsxfun(@times,sample_weights,sh);

% multiply with the transpose
hf_out1 = cell(1,1,num_features);
hf_out1{1} = permute(conj(mtimesx(sh, 'C', samplesf{1}, 'speed')), [3 4 2 1]);
for k = 2:num_features
    hf_out1{k} = permute(conj(mtimesx(sh(:,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end), 'C', samplesf{k}, 'speed')), [3 4 2 1]);
end

% compute the operation corresponding to the regularization term (convolve
% each feature dimension with the DFT of w, and the tramsposed operation)
% add the regularization part
% hf_conv = cell(1,1,num_features);
for k = 1:num_features
    reg_pad = min(size(reg_filter{k},2)-1, size(hf{k},2)-1);

    % add part needed for convolution
    hf_conv = cat(2, hf{k}, conj(rot90(hf{k}(:, end-reg_pad:end-1, :), 2)));

    % do first convolution
    hf_conv = convn(hf_conv, reg_filter{k});

    % do final convolution and put toghether result
    hf_out1{k} = hf_out1{k} + convn(hf_conv(:,1:end-reg_pad,:), reg_filter{k}, 'valid');
end

% Stuff related to the projection matrix

% Set the filter to the current one instead (test)
% init_hf = cellfun(@(hf) permute(hf, [3 4 1 2]), hf, 'uniformoutput', false);

% B * P
BP_cell = cell(1,1,num_features);
for k = 1:num_features
    BP_cell{k} = mtimesx(mtimesx(init_samplef{k}, P{k}, 'speed'), init_hf{k}, 'speed');
end

BP = BP_cell{1};
for k = 2:num_features
    BP(1,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) = ...
        BP(1,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end) + BP_cell{k};
end

% multiply with the transpose: A^H * BP
hf_out{1,1,1} = hf_out1{1} +  permute(bsxfun(@times, BP, conj(samplesf{1})), [3 4 2 1]);

% B^H * BP
fBP = cell(1,1,num_features);
fBP{1} = reshape(bsxfun(@times, conj(init_hf{1}), BP), size(init_hf{1},1), []).';

% Compute proj matrix part: B^H * A_m * f
shBP = cell(1,1,num_features);
shBP{1} = reshape(bsxfun(@times, conj(init_hf{1}), sh), size(init_hf{1},1), []).';

for k = 2:num_features
    % multiply with the transpose: A^H * BP
    hf_out{1,1,k} = hf_out1{k} +  permute(bsxfun(@times, BP(1,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end), conj(samplesf{k})), [3 4 2 1]);

    % B^H * BP
    fBP{k} = reshape(bsxfun(@times, conj(init_hf{k}), BP(1,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end)), size(init_hf{k},1), []).';

    % Compute proj matrix part: B^H * A_m * f
    shBP{k} = reshape(bsxfun(@times, conj(init_hf{k}), sh(1,1,1+pad_sz{k}(1):end-pad_sz{k}(1), 1+pad_sz{k}(2):end)), size(init_hf{k},1), []).';
end

% hf_out2 = cell(1,1,num_features);
for k = 1:num_features
    fi = size(hf{k},1) * (size(hf{k},2)-1) + 1;  % index where the last frequency column starts

    % B^H * BP
    hf_out2 = 2*real(XH{k} * fBP{k} - XH{k}(:,fi:end) * fBP{k}(fi:end,:)) + proj_reg * P{k};

    % Compute proj matrix part: B^H * A_m * f
    hf_out{2,1,k} = hf_out2 + (2*real(XH{k} * shBP{k} - XH{k}(:,fi:end) * shBP{k}(fi:end,:)));
end
end

repo: CFWCR


text:
%kcxinspector 'Interactively inspects results of classification '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros cxinspector.pane file
%
% Parameters:
% InputFile: iorig 'Input Original Image', required: 'Input Original Image'
% InputFile: iclass 'Input Classified Image', required: 'Input Classified Image'
% InputFile: idist 'Input Probabilities', optional: 'input distances (probabilities) file'
% InputFile: icsf 'Class Specification File', required: 'input class specification file'
%
% Example: kcxinspector({iorig, iclass, idist, icsf}, {'iorig','';'iclass','';'idist','';'icsf',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% cxinspector - Interactively inspects results of classification
%
%  DESCRIPTION
% Inspector (cxinspector xvroutine) allows the user to interactively inspect the results for a classification using the classifiers in the Classify Toolbox. It can be used to verify pixels in the image for classification errors and results. The user can click on pixels on the classified result to see its original values and classification results.
% This xvroutine expects as inputs the original image [-iorig], the classified result [-iclass], and the Class Specification File for that classification task [-icsf], which will be used to get the classes' labels.
% Optionally the Inspector xvroutine can use the output for the Output Probabilities that most classifiers can create which are related with the discriminant functions for that classifier. If these values are presented as input for the parameter [-idist], their values will be presented in the output window.
% The xvroutine will have a workspace for display of images that can display either the original image or the classification result. When the mouse is clicked on the workspace, a small text report will be output'd on the Output Window (which is shown when the Output button is pressed) with basic information about the classifier, the position where the user clicked the mouse, the original input data at that position and the classification result for that pixel.
% There is a Image Display button which will bring a panel where the user can set whether the original image or the classification result will be displayed, whether the routine will use a private colormap for display and if the images are multiband (i.e. images with depth and/or time > 1) which band will be displayed.
%
%
%
%  EXAMPLES
% All examples for the Classify toolbox are listed on the Classify Toolbox Manual. For an example of this routine, please see the Classify:workspaces:MINDIST example workspace.
%
%  "SEE ALSO"
% cthematicmap, Class Specification File format in the Classify Toolbox Manual.
%
%  RESTRICTIONS
% Expect the parameter [-icsf] to point to a valid Class Specification file, and the parameters [-iclass] to be a result of classification with one of the classifiers in this toolbox.
%
%  REFERENCES
% All references for the Classify toolbox are listed on the Classify Toolbox Manual.
%
%  COPYRIGHT
% Copyright (C) 1997 Rafael Santos. Khoros (C) Khoral Research, Inc.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function res = parallelPlane(plane, point)
%PARALLELPLANE Parallel to a plane through a point or at a given distance.
%
%   PL2 = parallelPlane(PL, PT)
%   Constructs the plane parallel to plane PL and containing the point PT.
%
%   PL2 = parallelPlane(PL, D)
%   Constructs the plane parallel to plane PL, and located at the given
%   signed distance D.
%
%   Example
%     % Create a plane normal to the 3D vector DIR
%     dir = [3 4 5];
%     plane = createPlane([3 4 5], dir);
%     % Create plane at a specific distance
%     plane2 = parallelPlane(plane, 5);
%     % Create a line perpendicular to both planes
%     line = [2 4 1 3 4 5];
%     pi1 = intersectLinePlane(line, plane);
%     pi2 = intersectLinePlane(line, plane2);
%     % check the distance between intersection points
%     distancePoints3d(pi1, pi2)
%     ans =
%         5
%
%   See also
%   geom3d, parallelLine3d

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2012-08-22, using Matlab 7.9.0.529 (R2009b)
% Copyright 2012-2022 INRA - Cepia Software Platform

if size(point, 2) == 1
    % use a distance. Compute position of point located at distance DIST on
    % the line normal to the plane.
    normal = normalizeVector3d(planeNormal(plane));
    point = plane(:, 1:3) + bsxfun(@times, point, normal);
end

% change origin, and keep direction vectors
res = [point plane(:, 4:9)];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function state = ...
    psoboundssoft(state,Aineq,bineq,Aeq,beq,LB,UB,nonlcon,options)

x = state.Population ;
% v = state.Velocities ;

for i = 1:size(state.Population,1)
    lowindex = [] ; highindex = [] ;
    if ~isempty(LB), lowindex = x(i,:) < LB ; end
    if ~isempty(UB), highindex = x(i,:) > UB ; end

    outofbounds = any([lowindex,highindex]) ;
    if ~outofbounds && ~isempty(Aineq) % Check linear inequalities
        outofbounds = any(Aineq*x(i,:)' - bineq > options.TolCon) ;
    end % if ~isempty
    if ~outofbounds && ~isempty(Aeq) % Check linear equalities
        outofbounds = any(abs(Aeq*x(i,:)' - beq) > options.TolCon) ;
    end % if ~isempty
    if ~outofbounds && ~isempty(nonlcon) % Nonlinear constraint check
        [c,ceq] = nonlcon(x(i,:)) ;
        outofbounds = any(c > options.TolCon) ;
        outofbounds = outofbounds || any(abs(ceq) > options.TolCon) ;
    end

    if outofbounds
        state.Score(i) = inf ;
    end % if outofbounds

    state.OutOfBounds(i) = outofbounds ;
end % for i

state.Population = x ;
% state.Velocities = v ;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout=lpccc2ar(varargin)
%
% For calling details please see v_lpccc2ar.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpccc2ar.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpccc2ar.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpccc2ar(varargin{:});
else
    v_lpccc2ar(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [tmpTrace] = plotSignalsGraph(IcaTraces,varargin)
	% Plots signals, offsetting by a fixed amount.
	% Biafra Ahanonu
	% started: 2013.11.02
	% inputs
		%
	% outputs
		%

	[tmpTrace] = ciapkg.view.plotSignalsGraph(IcaTraces,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function airy_bi_values_test ( )

%*****************************************************************************80
%
%% AIRY_BI_VALUES_TEST demonstrates the use of AIRY_BI_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'AIRY_BI_VALUES_TEST:\n' );
  fprintf ( 1, '  AIRY_BI_VALUES stores values of \n' );
  fprintf ( 1, '  the Airy function Bi(X).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           Bi(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, bi ] = airy_bi_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, bi );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mean = gamma_mean ( a, b, c )

%*****************************************************************************80
%
%% GAMMA_MEAN returns the mean of the Gamma PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    0.0 < B,
%    0.0 < C.
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = a + b * c;

  return
end


text:
function result = calcMapTopkMapTopkPreTopkRecLabel(queryLabel, retrievalLabel, qB, rB, topk)
%% Function: calcMapTopkMapTopkPreTopkRecLabel:
%   calculate Map, Topk map, Topk precision, Topk recall for hamming ranking task.
% Input:
%   queryLabel: 0-1 label matrix (numQuery * numLabel) for query set.
%   retrievalLabel: 0-1 label matrix (numQuery * numLabel) for retrieval set.
%   qB: compressed binary code for query set.
%   rB: compressed binary code for retrieval set.
%   topk (optional): vector for different (non-zero and ascending) topk.
% Output:
%   result.map: map for whole retrieval set
%   result.topkMap: vector. topk-Map for different topk
%   result.topkPre: vector. topk-Precision for different topk
%   result.topkRec: vector. topk-Recall for different topk

flag = false;
if exist('topk', 'var')
    flag = true;
end

numQuery = size(qB, 1);
map = 0;
if flag
    nt = numel(topk);
    topkPre = zeros(1, nt);
    topkRec = zeros(1, nt);
    topkMap = zeros(1, nt);
end

for ii = 1: numQuery
    gnd = queryLabel(ii, :) * retrievalLabel' > 0;
    tsum = sum(gnd);
    if tsum == 0
        continue;
    end
    hamm = hammingDist(qB(ii, :), rB);
    [~, index] = sort(hamm);
    gnd = gnd(index);
    count = 1: tsum;
    tindex = find(gnd == 1);
    map = map + mean(count ./ tindex);
    if flag
        for jj = 1: nt
            tgnd = gnd(1: topk(jj));
            if sum(tgnd) == 0
                continue;
            end
            tcount = 1: sum(tgnd);
            tindex = find(tgnd == 1);
            topkMap(jj) = topkMap(jj) + mean(tcount ./ tindex);
            topkPre(jj) = topkPre(jj) + sum(tgnd) / topk(jj);
            topkRec(jj) = topkRec(jj) + sum(tgnd) / tsum;
        end
    end
end

result.map = map / numQuery;
if flag
    result.topkMap = topkMap / numQuery;
    result.topkPre = topkPre / numQuery;
    result.topkRec = topkRec / numQuery;
end
end

author: jiangqy


text:
function [solver,problem,forced_choice,failureMode] = selectsolver(options,ProblemClass,solvers,socp_are_really_qc,allsolvers);
%SELECTSOLVER Internal function to select solver based on problem category

problem = 0;
failureMode = '';

% UNDOCUMENTED
force_solver = yalmip('solver');
if length(force_solver)>0
    options.solver = force_solver;
end

% YALMIP has discovered in an previous call that the model isn't a GP, and
% now searches for a non-GP solver
if options.thisisnotagp
    ProblemClass.gppossible = 0;
end

% ***************************************************
% Maybe the user is stubborn and wants to pick solver
% ***************************************************
forced_choice = 0;
if length(options.solver)>0 & isempty(strfind(options.solver,'*'))

    if strfind(options.solver,'+')
        forced_choice = 1;
        options.solver = strrep(options.solver,'+','');
    end
    % Create tags with version also
    temp = expandSolverName(solvers);

    opsolver = lower(options.solver);
    splits = strfind(opsolver,',');
    if isempty(splits)
        names{1} = opsolver;
    else
        start = 1;
        for i  = 1:length(splits)
            names{i} = opsolver(start:splits(i)-1);
            start = splits(i)+1;
        end
        names{end+1} = opsolver(start:end);
    end

    index1 = [];
    index2 = [];
    for i = 1:length(names)
        index1 = [index1 find(strcmpi({solvers.tag},names{i}))];
        index2 = [index1 find(strcmpi({temp.tag},names{i}))];
    end
    if isempty(index1) & isempty(index2)
        % Specified solver not found among available solvers
        % Is it even a supported solver
        temp = expandSolverName(allsolvers);
        for i = 1:length(names)
            index1 = [index1 find(strcmp(lower({allsolvers.tag}),names{i}))];
            index2 = [index1 find(strcmp(lower({temp.tag}),names{i}))];
        end
        if isempty(index1) & isempty(index2)
            problem = -9;
        else
            problem = -3;
        end
        solver = [];
        return;
    else
        solvers = solvers(union(index1,index2));
    end
end

% ************************************************
% Prune based on forced global solver
% ************************************************
if  ~forced_choice && options.forceglobal
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).global;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'global optimization');
end
% Prune based on avoiding global solver (in BMIBNB we don't want that as
% upper solver)
if isfield(options,'avoidGlobalSolver')
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = ~solvers(i).global;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'global optimization');
end

% ************************************************
% Prune based on objective
% ************************************************
if ProblemClass.objective.sigmonial & ~forced_choice & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).objective.sigmonial;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'signomial terms in objective');
end

% ************************************************
% Prune based on objective
% ************************************************
if ProblemClass.objective.sigmonial & ~forced_choice & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).objective.sigmonial;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'signomial terms in objective');
end
if ProblemClass.objective.polynomial & ~forced_choice & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).objective.polynomial | solvers(i).objective.sigmonial;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'polynomial terms in objective');
end
if ProblemClass.objective.quadratic.nonconvex & ~forced_choice & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).objective.polynomial | solvers(i).objective.sigmonial | solvers(i).objective.quadratic.nonconvex;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'nonconvex quadratic terms in objective');
end
if ProblemClass.objective.quadratic.convex & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        direct = solvers(i).objective.polynomial | solvers(i).objective.sigmonial | solvers(i).objective.quadratic.nonconvex | solvers(i).objective.quadratic.convex;
        indirect = solvers(i).constraint.inequalities.semidefinite.linear | solvers(i).constraint.inequalities.secondordercone.linear;
        if direct | indirect
            keep(i)=1;
        else
            keep(i)=0;
        end
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'quadratic terms in objective');
end
if ProblemClass.objective.linear & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).objective.polynomial | solvers(i).objective.sigmonial | solvers(i).objective.quadratic.nonconvex | solvers(i).objective.quadratic.convex | solvers(i).objective.linear;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'linear terms in objective');
end

if ProblemClass.objective.maxdet.convex & ~ProblemClass.objective.linear & ~ProblemClass.objective.quadratic.convex & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).objective.maxdet.convex | solvers(i).constraint.inequalities.semidefinite.linear;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'maxdet terms in objective');
end

if ProblemClass.objective.maxdet.convex & (ProblemClass.objective.linear | ProblemClass.objective.quadratic.convex) & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).objective.maxdet.convex;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'maxdet terms in objective');
end

if ProblemClass.objective.maxdet.nonconvex & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).objective.maxdet.nonconvex;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'nonconvex maxdet terms in objective');
end

% ******************************************************
% Prune based on rank constraints
% ******************************************************
if ProblemClass.constraint.inequalities.rank & ~forced_choice & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.rank;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'rank constraints');
end

% ******************************************************
% Prune based on semidefinite constraints
% ******************************************************
if ProblemClass.constraint.inequalities.semidefinite.sigmonial & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.semidefinite.sigmonial;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'signomial semidefinite constraints');
end
if ProblemClass.constraint.inequalities.semidefinite.polynomial & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.semidefinite.sigmonial |  solvers(i).constraint.inequalities.semidefinite.polynomial;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'polynomial semidefinite constraints');
end
if ProblemClass.constraint.inequalities.semidefinite.quadratic & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.semidefinite.sigmonial |  solvers(i).constraint.inequalities.semidefinite.polynomial | solvers(i).constraint.inequalities.semidefinite.quadratic;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'quadratic semidefinite constraints');
end
if ProblemClass.constraint.inequalities.semidefinite.linear & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.semidefinite.sigmonial |  solvers(i).constraint.inequalities.semidefinite.polynomial | solvers(i).constraint.inequalities.semidefinite.quadratic | solvers(i).constraint.inequalities.semidefinite.linear;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'semidefinite constraints');
end

% If user has specified a, e.g., LP solver for an SDP when using OPTIMIZER,
% we must bail out, as there is no chance this model instantiates as an LP.
if forced_choice &  (ProblemClass.constraint.inequalities.semidefinite.linear | ProblemClass.constraint.inequalities.semidefinite.quadratic | ProblemClass.constraint.inequalities.semidefinite.polynomial | ProblemClass.constraint.inequalities.semidefinite.sigmonial)  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.semidefinite.sigmonial |  solvers(i).constraint.inequalities.semidefinite.polynomial | solvers(i).constraint.inequalities.semidefinite.quadratic | solvers(i).constraint.inequalities.semidefinite.linear;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'semidefinite constraints');
end
% Similarily, we have a SOCP by definition. We must support that
if forced_choice & ~socp_are_really_qc & (ProblemClass.constraint.inequalities.secondordercone.linear | ProblemClass.constraint.inequalities.secondordercone.nonlinear)  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.secondordercone.linear;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'second-order cone constraints');
end

% ******************************************************
% Prune based on cone constraints
% ******************************************************
if ProblemClass.constraint.inequalities.secondordercone.linear & ~socp_are_really_qc & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
         keep(i) = solvers(i).constraint.inequalities.secondordercone.linear | solvers(i).constraint.inequalities.semidefinite.linear | solvers(i).constraint.inequalities.elementwise.quadratic.nonconvex;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'second-order cone constraints');
end
if ProblemClass.constraint.inequalities.secondordercone.nonlinear & ~socp_are_really_qc & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
         keep(i) = solvers(i).constraint.inequalities.secondordercone.nonlinear | solvers(i).constraint.inequalities.elementwise.quadratic.nonconvex;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'nonlinear second-order cone constraints');
end
if ProblemClass.constraint.inequalities.rotatedsecondordercone.linear & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.rotatedsecondordercone.linear | solvers(i).constraint.inequalities.secondordercone.linear | solvers(i).constraint.inequalities.semidefinite.linear;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'rotated second-order cone constraints');
end
if ProblemClass.constraint.inequalities.rotatedsecondordercone.nonlinear & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.rotatedsecondordercone.nonlinear | solvers(i).constraint.inequalities.secondordercone.linear | solvers(i).constraint.inequalities.semidefinite.linear;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'nonlinear rotated second-order cone constraints');
end
% if ProblemClass.constraint.inequalities.powercone && ~forced_choice  && length(solvers)>0
%     keep = ones(length(solvers),1);
%     for i = 1:length(solvers)
%         keep(i) = solvers(i).powercone;
%     end
%     [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'power cone constraints');
% end

% ******************************************************
% Prune based on element-wise inequality constraints
% ******************************************************
if ProblemClass.constraint.inequalities.elementwise.sigmonial & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.elementwise.sigmonial;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'signomial constraints');
end
if ProblemClass.constraint.inequalities.elementwise.polynomial & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.elementwise.sigmonial |  solvers(i).constraint.inequalities.elementwise.polynomial;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'polynomial constraints');
end
if ProblemClass.constraint.inequalities.elementwise.quadratic.nonconvex & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.elementwise.sigmonial |  solvers(i).constraint.inequalities.elementwise.polynomial | solvers(i).constraint.inequalities.elementwise.quadratic.nonconvex;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'nonconvex quadratic constraints');
end
if ProblemClass.constraint.inequalities.elementwise.quadratic.convex | (ProblemClass.constraint.inequalities.secondordercone.linear & socp_are_really_qc) & ~forced_choice & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.elementwise.sigmonial |  solvers(i).constraint.inequalities.elementwise.polynomial | solvers(i).constraint.inequalities.elementwise.quadratic.nonconvex | solvers(i).constraint.inequalities.elementwise.quadratic.convex | solvers(i).constraint.inequalities.secondordercone.linear | solvers(i).constraint.inequalities.semidefinite.linear;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'quadratic constraints');
end
if ProblemClass.constraint.inequalities.elementwise.linear & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.elementwise.sigmonial |  solvers(i).constraint.inequalities.elementwise.polynomial | solvers(i).constraint.inequalities.semidefinite.quadratic | solvers(i).constraint.inequalities.elementwise.linear;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'constraints');
end

% ******************************************************
% Prune based on element-wise constraints
% ******************************************************
if ProblemClass.constraint.equalities.sigmonial & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.inequalities.elementwise.sigmonial | solvers(i).constraint.equalities.sigmonial;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'signomial equality constraints');
end
if ProblemClass.constraint.equalities.polynomial & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        indirect = solvers(i).constraint.inequalities.elementwise.sigmonial | solvers(i).constraint.inequalities.elementwise.polynomial;
        indirect = indirect | solvers(i).constraint.inequalities.elementwise.sigmonial | solvers(i).constraint.inequalities.elementwise.polynomial;
        direct = solvers(i).constraint.equalities.sigmonial | solvers(i).constraint.equalities.polynomial;
        keep(i) = direct | indirect;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'polynomial equality constraints');
end
if ProblemClass.constraint.equalities.quadratic & ~forced_choice & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        indirect = solvers(i).constraint.inequalities.elementwise.sigmonial | solvers(i).constraint.inequalities.elementwise.polynomial | solvers(i).constraint.inequalities.elementwise.quadratic.nonconvex;
        direct = solvers(i).constraint.equalities.sigmonial |  solvers(i).constraint.equalities.polynomial |  solvers(i).constraint.equalities.quadratic;
        keep(i) = direct | indirect;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'quadratic equality constraints');
end
if ProblemClass.constraint.equalities.linear & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        indirect = solvers(i).constraint.inequalities.elementwise.linear | solvers(i).constraint.inequalities.elementwise.sigmonial |  solvers(i).constraint.inequalities.elementwise.polynomial;
        direct = solvers(i).constraint.equalities.linear | solvers(i).constraint.equalities.sigmonial |  solvers(i).constraint.equalities.polynomial;
        keep(i) = direct | indirect;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'equality constraints');
end

% ******************************************************
% Discrete data
% ******************************************************
if ProblemClass.constraint.integer & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).constraint.integer;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'integer variables');
end
if ProblemClass.constraint.binary & ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
         keep(i) = solvers(i).constraint.integer | solvers(i).constraint.binary;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'binary variables');
end
if ProblemClass.constraint.sos1 & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
         %keep(i) = solvers(i).constraint.integer | solvers(i).constraint.binary | solvers(i).constraint.sos2;
         keep(i) = solvers(i).constraint.sos2;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'sos1 constraints');
end
if ProblemClass.constraint.sos2 & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
         keep(i) = solvers(i).constraint.integer | solvers(i).constraint.binary | solvers(i).constraint.sos2;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'sos2 constraints');
end
if ProblemClass.constraint.semicont  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
         %keep(i) = solvers(i).constraint.integer | solvers(i).constraint.binary | solvers(i).constraint.sos2;
         keep(i) = solvers(i).constraint.semivar;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'semicontinuous variables');
end

% ******************************************************
% Equalities with multiple monomoials (rule out GP)
% ******************************************************
if ProblemClass.constraint.equalities.multiterm  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
         keep(i) = solvers(i).constraint.equalities.multiterm;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'multi-term monomial equalities');
end
% FIXME
% No support for multiterm is YALMIPs current way of saying "GP solver". We
% use this flag to prune GPs based on objective too
if ~ProblemClass.gppossible  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
         keep(i) = solvers(i).constraint.equalities.multiterm;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'multi-term monomial equalities');
end

% ******************************************************
% Complementarity constraints
% ******************************************************
if ProblemClass.constraint.complementarity.linear & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
         keep(i) = solvers(i).constraint.complementarity.linear | solvers(i).constraint.integer | solvers(i).constraint.binary | solvers(i).constraint.equalities.polynomial | solvers(i).constraint.equalities.quadratic;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'complementarity constraints');
end

% ******************************************************
% Interval data
% ******************************************************
if ProblemClass.interval  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = solvers(i).interval;
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'interval data');
end

% ******************************************************
% Parametric problem
% ******************************************************
if ~forced_choice  & length(solvers)>0
    keep = ones(length(solvers),1);
    for i = 1:length(solvers)
        keep(i) = (ProblemClass.parametric == solvers(i).parametric);
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'multiparametric models');
end

% ******************************************************
% Exponential cone representable (exp, log,...)
% ******************************************************
keep = ones(length(solvers),1);
if ~forced_choice  & length(solvers)>0
    for i = 1:length(solvers)
        keep(i) = (ProblemClass.exponentialcone <= solvers(i).exponentialcone) || (ProblemClass.exponentialcone <= solvers(i).evaluation);
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'exponential cones');
end
% ******************************************************
% Power cone representable
% ******************************************************
keep = ones(length(solvers),1);
if ~forced_choice  & length(solvers)>0
    for i = 1:length(solvers)
        keep(i) = (ProblemClass.constraint.inequalities.powercone <= solvers(i).powercone) || (ProblemClass.constraint.inequalities.powercone <= solvers(i).evaluation);
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'exponential cones');
end

% ******************************************************
% General functions (sin, cos,...)
% ******************************************************
keep = ones(length(solvers),1);
if ~forced_choice  & length(solvers)>0
    for i = 1:length(solvers)
        keep(i) = (ProblemClass.evaluation <= solvers(i).evaluation) || (ProblemClass.exponentialcone && solvers(i).exponentialcone);
    end
    [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,'general nonlinearities');
end

% FIX : UUUUUUGLY
if isempty(solvers)
    solver = [];
else
    if length(options.solver)>0
        solver = [];

        % FIX : Re-use from above
        opsolver = lower(options.solver);
        splits = strfind(opsolver,',');
        if isempty(splits)
            names{1} = opsolver;
        else
            names = {};
            start = 1;
            for i  = 1:length(splits)
                names{i} = opsolver(start:splits(i)-1);
                start = splits(i)+1;
            end
            names{end+1} = opsolver(start:end);
        end

        temp = solvers;
        for i = 1:length(temp)
        if length(temp(i).version)>0
            temp(i).tag = lower([temp(i).tag '-' temp(i).version]);
        end
        end

        for i = 1:length(names)
            if isequal(names{i},'*')
                solver = solvers(1);
                break
            else
                j = find(strcmpi(lower({solvers.tag}),names{i}));
                if ~isempty(j)
                    solver = solvers(j(1));
                    break
                end
                j = find(strcmpi(lower({temp.tag}),names{i}));
                if ~isempty(j)
                    solver = solvers(j(1));
                    break
                end
            end
        end
    else
        solver = solvers(1);
    end
end

if isempty(solver)
    if length(options.solver)>0 % User selected available solver, but it is not applicable
        problem = -4;
    else
        problem = -2;
    end
end

% FIX : Hack when chosing the wrong fmincon thingy
if ~isempty(solver)
    c1 = (length(options.solver)==0 | isequal(lower(options.solver),'fmincon')) & isequal(lower(solver.tag),'fmincon') & isequal(solver.version,'geometric');
    c2 = (length(options.solver)==0 | isequal(lower(options.solver),'snopt')) & isequal(lower(solver.tag),'snopt') & isequal(solver.version,'geometric');
    if c1 | c2
        if ~(ProblemClass.objective.sigmonial | ProblemClass.constraint.inequalities.elementwise.sigmonial)
            solver.version = 'standard';
            solver.call    = strrep(solver.call,'gp','');
            solver.objective.linear = 1;
            solver.objective.quadratic.convex = 1;
            solver.objective.quadratic.nonconvex = 1;
            solver.objective.polynomial = 1;
            solver.objective.sigmonial = 1;
            solver.constraint.equalities.elementwise.linear = 1;
            solver.constraint.equalities.elementwise.quadratic.convex = 1;
            solver.constraint.equalities.elementwise.quadratic.nonconvex = 1;
            solver.constraint.equalities.elementwise.polynomial = 1;
            solver.constraint.equalities.elementwise.sigmonial = 1;
            solver.constraint.inequalities.elementwise.linear = 1;
            solver.constraint.inequalities.elementwise.quadratic.convex = 1;
            solver.constraint.inequalities.elementwise.quadratic.nonconvex = 1;
            solver.constraint.inequalities.elementwise.polynomial = 1;
            solver.constraint.inequalities.elementwise.sigmonial = 1;
            solver.constraint.inequalities.secondordercone.linear = 1;
            solver.constraint.inequalities.semidefinite.linear = 1;
            solver.constraint.inequalities.semidefinite.quadratic = 1;
            solver.constraint.inequalities.semidefinite.polynomial = 1;
            solver.constraint.inequalities.semidefinite.nonlinear = 1;
            solver.exponentialcone = 1;
            solver.powercone = 1;
            solver.dual = 1;
            solver.evaluation = 1;
        end
    end
end

function temp = expandSolverName(temp)
for i = 1:length(temp)
    if length(temp(i).version)>0
        temp(i).tag = lower([temp(i).tag '-' temp(i).version]);
    end
end

function [solvers,failureMode] = pruneBasedOn(solvers,keep,failureMode,reason)

solvers = solvers(find(keep));
if ~any(keep)
    failureMode = reason;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [h,Hf,Hh] = fromFrameHmg(f,hf)

% FROMFRAMEHMG  Fom-frame transformation for homogeneous coordinates
%   P = FROMFRAMEHMG(F,PF) transforms homogeneous point PF from frame F to
%   the global frame.
%
%   [p,Pf,Ppf] = ... returns the Jacobians wrt F and PF.

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

F = homogeneous(f)    ;
[t,q] = splitFrame(f) ;

h = F*hf;

if nargout > 1

    [a,b,c,d] = split(q);
    [hx,hy,hz,ht] = split(hf);

    Ht = [...
        [ ht,  0,  0]
        [  0, ht,  0]
        [  0,  0, ht]
        [  0,  0,  0]];
    Hq = [...
        [  2*a*hx-2*d*hy+2*c*hz,  2*b*hx+2*c*hy+2*d*hz, -2*c*hx+2*b*hy+2*a*hz, -2*d*hx-2*a*hy+2*b*hz]
        [  2*d*hx+2*a*hy-2*b*hz,  2*c*hx-2*b*hy-2*a*hz,  2*b*hx+2*c*hy+2*d*hz,  2*a*hx-2*d*hy+2*c*hz]
        [ -2*c*hx+2*b*hy+2*a*hz,  2*d*hx+2*a*hy-2*b*hz, -2*a*hx+2*d*hy-2*c*hz,  2*b*hx+2*c*hy+2*d*hz]
        [                     0,                     0,                     0,                     0]];

    Hf = [Ht Hq];
    Hh = F;

end

return

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
classdef KTA2 < ALGORITHM
% <multi/many> <real/integer> <expensive>
% Kriging-assisted Two_Arch2
% tau  --- 0.75 --- Proportion of one type noninfluential points in training data
% phi  ---  0.1 --- Number of randomly selected individuals
% wmax ---   10 --- Number of generations before updating CA and DA
% mu   ---    5 --- Number of re-evaluated solutions at each generation

%------------------------------- Reference --------------------------------
% Z. Song, H. Wang, C. He and Y. Jin, A Kriging-assisted two-archive
% evolutionary algorithm for expensive many-objective optimization, IEEE
% Transactions on Evolutionary Computation, 2021, 25(6): 1013-1027.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Zhenshou Song
% Email:zssong@stu.xidian.edu.cn

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [tau,phi,wmax,mu] = Algorithm.ParameterSet(0.75,0.1,10,5);

            %% Initialization
            p = 1/Problem.M;
            CAsize = Problem.N;
            N = Problem.N;
            P = UniformPoint(N, Problem.D, 'Latin');
            Population = Problem.Evaluation(repmat(Problem.upper-Problem.lower,N,1).*P+repmat(Problem.lower,N,1));
            All_Population = Population;
            Ho_Population = All_Population;
            CA = UpdateCA([],Population,CAsize);
            DA = Population;
            THETA_S = 5.*ones(Problem.M,Problem.D);
            THETA_IS =  5.*ones(2,Problem.M,Problem.D);
            Model_sensitive = cell(1,Problem.M);
            Model_insensitive = cell(2,Problem.M);
            %% Optimization
            while Algorithm.NotTerminated(All_Population)
                %***** Building influential point-insensitive model********
                % build sensitive model
                Dec = All_Population.decs;
                Obj = All_Population.objs;
                for i = 1:Problem.M
                    dmodel     = dacefit(Dec,Obj(:,i),'regpoly0','corrgauss',THETA_S(i,:),1e-5.*ones(1,Problem.D),100.*ones(1,Problem.D));
                    Model_sensitive{i}   = dmodel;
                    THETA_S(i,:) = dmodel.theta;
                end
                % build insensitive models
                Centers = zeros(Problem.M,2);
                for i = 1 : Problem.M
                    [~,N1] = sort(Obj(:,i));
                    num = ceil(length(All_Population).*tau);
                    mean_index{1} = N1(1:num);
                    mean_index{2} = N1(end-num:end);
                    for j = 1:2
                        Centers(i,j) = mean(Obj(mean_index{j},i));  % lambda and miu
                    end
                    for j = 1 : 2
                        train_X = Dec(mean_index{j},:);
                        train_Y = Obj(mean_index{j},i);
                        dmodel  = dacefit(train_X,train_Y,'regpoly0','corrgauss',THETA_IS(j,i,:),1e-5.*ones(1,Problem.D),100.*ones(1,Problem.D));
                        Model_insensitive{j,i} = dmodel;
                        THETA_IS(j,i,:)        = dmodel.theta;
                    end
                end
                % Set the CCA and CDA as the current CA and DA
                CAobj = CA.objs; CAdec = CA.decs;
                DAobj = DA.objs; DAdec = DA.decs;
                w = 1;
                while w <= wmax   % this part is same as Two_Arch2
                    [~,ParentCdec,~,ParentMdec] = MatingSelection_KTA2(CAobj,CAdec,DAobj,DAdec,Problem.N);
                    OffspringDec = [OperatorGA(Problem,ParentCdec,{1,20,0,0});OperatorGA(Problem,ParentMdec,{0,0,1,20})];
                    PopDec = [DAdec;CAdec;OffspringDec];
                    N      = size(PopDec,1);
                    PopObj = zeros(N,Problem.M);
                    MSE    = zeros(N,Problem.M);
                    %****** Using influential point-insensitive model *****
                    for i = 1:N
                        for j = 1:Problem.M
                            [PopObj(i,j),~,~] = predictor(PopDec(i,:),Model_sensitive{j});
                            if abs(PopObj(i,j)- Centers(j,1)) <= abs(PopObj(i,j)- Centers(j,2))
                                model = Model_insensitive{1,j};
                            else
                                model = Model_insensitive{2,j};
                            end
                            [PopObj(i,j),~,MSE(i,j)] = predictor(PopDec(i,:),model);
                        end
                    end
                    [CAobj,CAdec,~] = K_UpdateCA(PopObj,PopDec,MSE,CAsize);
                    [DAobj,DAdec,DAvar] = K_UpdateDA(PopObj,PopDec,MSE,Problem.N,p);
                    w = w + 1;
                end

                % Adaptive sampling
                Offspring01 = Adaptive_sampling(CAobj,DAobj,CAdec,DAdec,DAvar,DA,mu,p,phi);

                [~,index] = unique(Offspring01 ,'rows');
                PopNew = Offspring01(index,:);
                Offspring02 = [];
                for i = 1:size(PopNew,1)
                    dist2 = pdist2(real( PopNew(i,:)),real(All_Population.decs));
                    if min(dist2) > 1e-5
                        Offspring02 = [Offspring02;PopNew(i,:)];
                    end
                end
                if ~isempty(Offspring02)
                    Offspring = Problem.Evaluation(Offspring02);

                    temp =  All_Population.decs;
                    for i = 1:size(Offspring,2)
                        dist2 = pdist2(real(Offspring(i).decs),real(temp));
                        if min(dist2) > 1e-5
                            All_Population = [All_Population,Offspring(i)];
                        end
                        temp = All_Population.decs;
                    end
                    CA = UpdateCA(CA,Offspring,CAsize);
                    DA = UpdateDA(DA,Offspring,Problem.N,p);
                    Ho_Population = [Ho_Population,Offspring];
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function rpatch = patchResize(opatch,r)

% PATCHRESIZE  Patch resize.
%   PATCHRESIZE(PATCH,R)  resizes the image of PATCH.I by a
%   factor R and crops the result to size(PATCH.I). Works only
%   for square patch images. See PIX2PATCH for information about
%   the structure PATCH.
%
%   See also PIX2PATCH, PATCHSCAN

% (c) 2006 Joan Sola @ LAAS-CNRS

I    = single(opatch.I);

% J = single(imresize2(I,r));
J = imresize2(I,r);

rpatch.I    = J;
rpatch.SI   = sum(sum(J));
rpatch.SII  = sum(sum(J.*J));
rpatch.bias = opatch.bias;

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE


text:
function [] = plot_nodes(ax, parsed_osm, only_node_indices, show_id)
% plot (selected) nodes and label each with index and id
%
% usage
%   PLOT_NODES(ax, parsed_osm, only_node_indices, show_id)
%
% input
%   ax = axes object handle where to plot the nodes.
%   parsed_osm = MATLAB structure containing the OpenStreetMap XML data
%                after parsing, as returned by function
%                parse_openstreetmap.
%   only_node_indices = selected node indices in the global node matrix.
%   show_id = select whether to show or not the ID numbers of nodes as text
%             labes within the plot
%           = 0 (do not show labels) | 1 (show labels)
%
% 2012.04.24 (c) Ioannis Filippidis, jfilippidis@gmail.com
%
% See also PARSE_OPENSTREETMAP, ROUTE_PLANNER.

% do not show node id (default)
if nargin < 4
    show_id = 0;
end

nodes = parsed_osm.node;
node_ids = nodes.id;
node_xys = nodes.xy;

% which nodes to plot ?
n = size(node_xys, 2);
if nargin < 3
    only_node_indices = 1:n;
end

%% plot
held = takehold(ax);

% nodes selected exist ?
if n < max(only_node_indices)
    warning('only_node_indices contains node indices which are too large.')
    return
end

% plot nodes
xy = node_xys(:, only_node_indices);
plotmd(ax, xy, 'yo')

% label plots
for i=only_node_indices
    node_id_txt = num2str(node_ids(1, i) );
    if show_id
        curtxt = {['index=', num2str(i) ], ['id=', node_id_txt] }.';
    else
        curtxt = ['index=', num2str(i) ];
    end
    textmd(node_xys(:, i), curtxt, 'Parent', ax)
end

restorehold(ax, held)

author: johnyf

repo: openstreetmap


text:
function [init_factors, init_factors_opts] = generate_init_factors(V, rank, options)
% Initialization algorithm
%
% Created by H.Kasai on May 21, 2019
%
% Change log:
%
%
%   May  21, 2019 (Hiroyuki Kasai): Created the initial version
%
%   Oct. 14, 2020 (H.Huangi): Added 'LPinit'
%
%   June. 20, 2022 (Hiroyuki Kasai): Modified to handle the case where one
%                                   factorization matrix is input, and
%                                   the other is not.

    % V = WH + R

    m = size(V, 1);
    n = size(V, 2);

    init_factors_opts = [];

    % initialize
    generate_wh_init = true;
    generate_r_init = true;

    %% generate W and H
    if isfield(options, 'x_init')
        if isfield(options.x_init, 'W') && isfield(options.x_init, 'H')
            init_factors.W = options.x_init.W;
            init_factors.H = options.x_init.H;
            generate_wh_init = false;
        elseif isfield(options.x_init, 'W') && ~isfield(options.x_init, 'H')
            init_factors.W = options.x_init.W;
            init_factors.H = rand(rank, n);
            generate_wh_init = false;
        elseif ~isfield(options.x_init, 'W') && isfield(options.x_init, 'H')
            init_factors.W = rand(m, rank);
            init_factors.H = options.x_init.H;
            generate_wh_init = false;
        end

        if isfield(options.x_init, 'R')
            init_factors.R = options.x_init.R;
            generate_r_init = false;
        end
    end

    if isfield(options, 'special_init_factors')

        if generate_wh_init || generate_r_init
            [init_factors, init_factors_opts] = options.special_init_factors(V, rank, generate_wh_init, init_factors_opts, options);
        end

    else

        if generate_wh_init

            if ~isfield(options, 'init_alg')
                alg = 'random';
            else
                alg = options.init_alg;
            end

            switch(alg)

                case 'LPinit'

                    H = LPinitSemiNMF(V, rank);
                    W = V * pinv(H);

                    init_factors.W = W;
                    init_factors.H = H;

                case 'random'

                    W = rand(m, rank);
                    H = rand(rank, n);

                    init_factors.W = W;
                    init_factors.H = H;

                case 'ones'

                    W = ones(m, rank);
                    H = ones(rank, n);

                    init_factors.W = W;
                    init_factors.H = H;

                case 'semi_random' % for SemiNMF

                    H = rand(rank, n);
                    W = V / H; % V * inv(H)

                    init_factors.W = W;
                    init_factors.H = H;

                case 'symm'

                    W = ones(m, rank);

                    init_factors.W = W;
                    init_factors.H = W';

                case 'symm_mean'

                    % make sure that entries of H fall into the interval [0, 2*sqrt(m/k)],
                    % where 'm' is the average of all entries of V.
                    % See https://github.com/dakuang/symnmf

                    W = 2 * full(sqrt(mean(mean(V)) / rank)) * rand(m, rank);

                    init_factors.W = W;
                    init_factors.H = W';

                case 'NNDSVD'

                    [W, H] = NNDSVD(abs(V), rank, 0);

                    init_factors.W = W;
                    init_factors.H = H;

                case 'kmeans'

                    [label, center] = litekmeans(V',rank, 'maxIter', 10);
                    center = max(0,center);
                    W = center';
                    WTW = W'*W;
                    WTW = max(WTW, WTW');
                    WTX = W'*V;
                    H = max(0, WTW\WTX);

                    init_factors.W = W;
                    init_factors.H = H;

                case {'prob_expectation', 'prob_random'}

                    T = tn_vector();

                    if strcmp(alg, 'prob_random')
                        E = exponential_dist();
                    end

                    exp_U   = options.exp_U;
                    exp_V   = options.exp_V;
                    var_U   = options.var_U;
                    var_V   = options.var_V;
                    mu_U    = options.mu_U;
                    mu_V    = options.mu_V;
                    tau_U   = options.tau_U;
                    tau_V   = options.tau_V;

                    for k = 1 : rank
                        for i = 1 : m
                            tau_U(i, k) = 1;

                            if options.ard
                                hyperparam = options.exp_lambdak(k);
                            else
                                hyperparam = options.hyperparams.lambdaU(i, k);
                            end

                            if strcmp(alg, 'prob_random')
                                mu_U(i, k) = E.exponential_draw(hyperparam);
                            else
                                mu_U(i, k) = 1.0/hyperparam;
                            end
                        end
                    end

                    for k = 1 : rank
                        for j = 1 : n
                            tau_V(j, k) = 1;

                            if options.ard
                                hyperparam = options.exp_lambdak(k);
                            else
                                hyperparam = options.hyperparams.lambdaV(j, k);
                            end

                            if strcmp(alg, 'prob_random')
                                mu_V(j, k) = E.exponential_draw(hyperparam);
                            else
                                mu_V(j, k) = 1.0/hyperparam;
                            end
                        end
                    end

                    for k = 1 : rank
                        [exp_U, var_U] = update_exp(T, exp_U, var_U, mu_U, tau_U, k);
                    end

                    for k = 1 : rank
                        [exp_V, var_V] = update_exp(T, exp_V, var_V, mu_V, tau_V, k);
                    end

                    init_factors.W = exp_U;
                    init_factors.H = exp_V';
                    init_factors_opts.var_W = var_U;
                    init_factors_opts.var_H = var_V';
                    init_factors_opts.mu_W  = mu_U;
                    init_factors_opts.mu_H  = mu_V';
                    init_factors_opts.tau_W = tau_U;
                    init_factors_opts.tau_H = tau_V';

                otherwise

                    % do random initialization

                    W = rand(m, rank);
                    H = rand(rank, n);

                    init_factors.W = W;
                    init_factors.H = H;

            end

        end

        if isfield(options, 'norm_w')
            if options.norm_w ~= 0
                % normalize W
                init_factors.W = normalize_W(init_factors.W, options.norm_w);
            end
        end

        if isfield(options, 'norm_h')
            if options.norm_h ~= 0
                % normalize H
                init_factors.H = normalize_H(init_factors.H, options.norm_h);
            end
        end

        %% generate R
        if generate_r_init
            if isfield(options, 'x_init_robust')
                if options.x_init_robust
                    if isfield(options, 'x_init')
                        if isfield(options.x_init, 'R')
                            init_factors.R = options.x_init.R;
                            generate_r_init = false;
                        end
                    end

                    if generate_r_init
                        init_factors.R = rand(m, n);
                    end
                else
                    init_factors.R = zeros(m, n);
                end
            else
                init_factors.R = zeros(m, n);
            end
        end

    end
end

function [exp_A, var_A] = update_exp(T, exp_A, var_A, mu_A, tau_A, k)

    exp_A(:,k) = T.expectation(mu_A(:,k), tau_A(:,k));
    var_A(:,k) = T.variance(mu_A(:,k), tau_A(:,k));

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function def = spm_get_def(Bx,By,Bz,beta)
% Calculating deformation field from coefficient vector.
%
% FORMAT [def] = get_def(Bx,By,Bz,beta);
% Bx, By, Bz  - Separable basis sets such that B=kron(Bz,kron(By,Bx));
% beta        - Coefficient vector for basis set.
% or
% FORMAT [def] = get_def(dim,order,beta);
% dim         - Dimensionality of image in x,y and z directions.
% order       - Order of DCT set in x,y and z directions.
% beta        - Coefficient vector for basis set.
%_______________________________________________________________________
%
% Calculates the equivalent to kron(Bz,kron(By,Bx))*beta
% in a faster and more efficient way. Note that this routine
% can be used to calculate AtY efficiently as well since
% A'*y = (diag(dfdy)*kron(Bz,kron(By,Bx)))'*y =
% kron(Bz',kron(By',Bx'))*(diag(dfdy)*y) = get_def(Bx',By',Bz',dfdy.*y)
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jesper Andersson
% $Id: spm_get_def.m 1143 2008-02-07 19:33:33Z spm $

if nargin == 4
   [nx,mx] = size(Bx);
   [ny,my] = size(By);
   [nz,mz] = size(Bz);
   if mx*my*mz ~= size(beta,1)
      warning('get_def: Size mismatch between beta and basis-set');
      return
   end
elseif nargin == 3
   if numel(Bx) ~= 3 || numel(By) ~= 3
      warning('get_def: Wrong dimensionality on input');
   elseif prod(By) ~= size(Bz,1)
      warning('get_def: Size mismatch between beta and basis-set');
      return
   end
   nx = Bx(1); ny = Bx(2); nz = Bx(3);
   mx = By(1); my = By(2); mz = By(3);
   beta = Bz;
   Bx = spm_dctmtx(nx,mx);
   By = spm_dctmtx(ny,my);
   Bz = spm_dctmtx(nz,mz);
end

def = zeros(nx*ny*nz,1);
for bf = 1:mz
   mbeta = reshape(beta((bf-1)*my*mx+1:bf*my*mx),mx,my);
   tmp = reshape(Bx*mbeta*By',nx*ny,1);
   for sl = 1:nz
      def((sl-1)*ny*nx+1:sl*ny*nx) = def((sl-1)*ny*nx+1:sl*ny*nx) + Bz(sl,bf)*tmp;
   end
end

return

author: spm


text:
function y = huber_pos( x, M, t )

%HUBER_POS   Monotonic huber penalty function.
%   For a vector X, HUBER_POS(X) computes the monotonic Huber-style function
%
%                      0     if 0>=X
%       HUBER_POS(X) = X^2   if 0<=X<=1
%                      2*X-1 if    X>=1
%
%   HUBER_POS(X,M) is the monotonic Huber-style penalty function of
%   halfwidth M, M.^2.*HUBER_POS(X./M). M must be real and positive.
%
%   HUBER_POS(X,M,T) computes the monotonic Huber-style penalty function
%   with halfwidth M and concomitant scale T:
%
%       HUBER_POS(X,M,T) = T.*HUBER_POS(X./T,M) if T > 0
%                          +Inf                 if T <= 0
%
%   See the help file for HUBER for information about this usage.
%
%   For matrices and N-D arrays, the penalty function is applied to each
%   element of X independently. M and T must be compatible with X in the same
%   sense as .*: one must be a scalar, or they must have identical size.
%
%   Disciplined convex programming information:
%       HUBER_POS is jointly convex in X and T. It is nondecreasing in X and
%       nonincreasing in T. Therefore, when used in CVX specifications, X
%       must be convex and T must be concave (or affine). Both must be real.

%
% Check arguments
%

narginchk(1,3);
if ~isreal( x ),
    error( 'First argument must be real.' );
end
if nargin < 2,
    M = 1;
elseif ~isreal( M ) || any( M( : ) <= 0 ),
    error( 'Second argument must be real and positive.' );
end
if nargin < 3,
    t = 1;
elseif ~isreal( t ),
    error( 'Third argument must be real.' );
end
sz = cvx_size_check( x, M, t );
if isempty( sz ),
    error( 'Sizes are incompatible.' );
end

%
% Compute result
%

y = max( x, 0 );
z = min( y, M );
y = t .* z .* ( 2 * y - z );
q = t <= 0;
if nnz( q ),
    if length( t ) == 1,
        y = Inf * ones( sy );
    else
        y( q ) = Inf;
    end
end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [t,u]=fpreadh(fname)
% FPREADH reads the unformatted data exported by FlexPDE HISTORY plot command.
%
% Description:
% This utility function reads the unformatted
% data EXPORTED by FlexPDE using the HISTORY
% plot command. This code has been tested under
% Matlb R13, R14, and FlexPDE 4.x.
%
% Syntax:
%       [sv,u]=fpreadh(fname)
%              fpreadh(fname)
%
% Input arguments:
% fname: a string contains the data file name (and its path).
%
% Output arguments:
% t: vector of time data points.
% u: matrix of values of u at exported grid points
%
% If fpreadh is called without output arguments it simply
% will draw data in a 2D plot.
%
% Limitaions:
%	1. Only one variable can be exported to each file.
%   2. Supported FlexPDE output files limited to unformatted data file.
%
%
% See also FPREADE, FPREADSC

%{

___________________________________________________
File Properties:
Name:		fpreadh.m
Subject:	Reading History plot data from FlexPDE
Category:	FPMat: FlexPDE Interface
Author:		Mohammad Rahmani
Created:	Dec. 3, 2004
Version: 	1.0
Comments:	See FPMat Interface.pdf for more information
______________________________________________________

Mohammad Rahmani
Chemical Eng. Dep.
Amir-Kabir Uni. of Technology
Tehran, IRAN
Mohammad.Rahmani@Gmail.com

%}

%% 1 - Check for input/output arguments
if nargin<1
    error('Not enough input arguments')
end
if nargout==0
    outputfcn='Plot';
elseif nargout~=2
    error('Wrong number of output arguments')
else
    outputfcn='SendOut';
end

%% 2- Open the Data file and catch the error
[fid,errmsg] = fopen(fname,'r');
if fid == -1 % An error has been occured
    error([fname, '!??  ',errmsg])
end;

%% 3- Bypass the header of data file
% This header normally is a FlexPDE header enclosed in {}
while (feof(fid) ~= 1)	% ignore lines till }
  tline = fgetl(fid);
  if (tline == '}')
    break
    end
end

%% 4- Read the rest of data  and close the file
svar=textscan(fid,'%f','headerLines',2);
nc=0;
while (feof(fid) ~=1)
    nc=nc+1;
    data(nc)=textscan(fid,'%f','headerLines',1);
end
fclose(fid);

%% 5- Convert cell data to double data matrix
svar=svar{1};
data=cell2mat(data);

%% 6 - Prepare output arguments
switch outputfcn
    case 'Plot'
        figure('Name',['History Plot for: ',fname]);
        plot(svar,data);
        return
    case 'SendOut'
        t=svar;
        u=data;
        return
end

author: Sable

repo: mcbench-benchmarks


text:
function reg_filter = get_reg_filter(sz, target_sz, params, reg_window_edge)

% Compute the spatial regularization function and derive the corresponding
% filter operation used for the optimization

if nargin < 3 || isempty(reg_window_edge)
    reg_window_edge = params.reg_window_edge;
end

if params.use_reg_window
    % create weight window
    reg_window_power = params.reg_window_power;

    % normalization factor
    reg_scale = 0.5 * target_sz;

    % construct grid
    wrg = -(sz(1)-1)/2:(sz(1)-1)/2;
    wcg = -(sz(2)-1)/2:(sz(2)-1)/2;
    [wrs, wcs] = ndgrid(wrg, wcg);

    % construct the regukarization window
    reg_window = (reg_window_edge - params.reg_window_min) * (abs(wrs/reg_scale(1)).^reg_window_power + abs(wcs/reg_scale(2)).^reg_window_power) + params.reg_window_min;

    % compute the DFT and enforce sparsity
    reg_window_dft = fft2(reg_window) / prod(sz);
    reg_window_dft(abs(reg_window_dft) < params.reg_sparsity_threshold * max(abs(reg_window_dft(:)))) = 0;

    % do the inverse transform, correct window minimum
    reg_window_sparse = real(ifft2(reg_window_dft));
    reg_window_dft(1,1) = reg_window_dft(1,1) - prod(sz) * min(reg_window_sparse(:)) + params.reg_window_min;
    reg_window_dft = fftshift(reg_window_dft);

    % find the regularization filter by removing the zeros
    reg_filter = cast(real(reg_window_dft(~all(reg_window_dft==0,2), ~all(reg_window_dft==0,1))), 'like', params.data_type);
else
    % else use a scaled identity matrix
    reg_filter = cast(params.reg_window_min, 'like', params.data_type);
end

author: martin-danelljan

repo: ECO


text:
function [metList,noMolMetList] = mol2sdf(mets,molfileDir,sdfFileName,includeRs)
% Concatenates molfiles in molfileDir into one SDF file.
%
% metList = mol2sdf(mets,molfileDir,sdfFileName)
%
% INPUTS
% mets              mx1 cell array of metabolite identifiers (e.g., BiGG
%                   abbreviations)
% molfileDir        Path to directory containing molfiles for metabolites
%                   in mets. Molfile names should match the metabolite
%                   identifiers in mets without compartment assignments.
%
% OPTIONAL INPUTS
% sdfFileName       Name of SDF file. Default is MetStructures.sdf.
% includeRs         {0,(1)}. If 0, variable structures such as R groups and
%                   repeat units will not be included in SDF.
%
% OUTPUTS
% metList           Cell array listing metabolites in SDF.
% noMolMetList      Cell array listing metabolites without mol file.
%
% WRITTEN OUTPUTS
% sdfFileName.sdf   SDF with metabolite structures in same order as in
%                   metList. Metabolite identifiers in the SDF are the same
%                   as in MetList.
%
% Ronan M.T. Fleming
% Hulda SH, Nov. 2012   Changed output format from CDF to SDF. Renamed
%                       function. Simplified code.

% Format inputs
if ~strcmp(molfileDir(end),filesep)
    molfileDir = [molfileDir filesep];
end

if ~exist('sdfFileName','var')
    sdfFileName = 'MetStructures.sdf';
end
if isempty(sdfFileName)
    sdfFileName = 'MetStructures.sdf';
end
if ~strcmp(sdfFileName(end-3:end),'.sdf')
    sdfFileName = [sdfFileName '.sdf'];
end

mets = reshape(mets,length(mets),1);
omets = mets; % store original metabolite identifiers
if ischar(mets)
    mets = strtrim(cellstr(mets));
end
if isnumeric(mets)
    mets = strtrim(cellstr(num2str(mets)));
end
mets = regexprep(mets,'(\[\w\])$',''); % Remove compartment assignment

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

% List all molfiles in molfile directory
d = dir(molfileDir);
dirbool = cat(1,d.isdir);
d = d(~dirbool);
molfileNames = {d.name};
molfileNames = regexprep(molfileNames,'(\.mol)$','');

molfileNames = molfileNames(ismember(molfileNames,mets)); % Only include molfiles for metabolites in mets

noMolMetList = unique(mets(ismember(mets,molfileNames)));
noMolMetList = reshape(noMolMetList,length(noMolMetList),1);

metList = molfileNames;

% List all elements
elements = {'H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', 'K', 'Ca',...
            'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', 'Rb', 'Sr', 'Y',...
            'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', 'Cs', 'Ba', 'La', 'Ce',...
            'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir',...
            'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', 'Ra', 'Ac', 'Th', 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm',...
            'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr', 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt', 'Ds', 'Rg', 'Cn', 'Uut', 'Fl', 'Uup', 'Lv', 'Uus', 'Uuo'};

% Concatenate molfiles into single SDF file
fid = fopen(sdfFileName,'w+'); % Create SDF
first = true;
for m = 1:length(molfileNames)
    molfile = fileread([molfileDir molfileNames{m} '.mol']); % Read in molfile
    molfile = reshape(molfile, 1, length(molfile)); % for compatibility with Octave (where fileread returns a column vector of chars)
    molfile = regexprep(molfile,'\r',''); % Remove carriage returns
    molfile = regexprep(molfile,'[^\n]*\n',sprintf('%s\n',molfileNames{m}),'once'); % Replace top line with metabolite identifier
    molfile = regexprep(molfile,'\n+$',''); % Ensure there is only one newline caracter at end of each molfile in the SDF

    if includeRs == 0
        % Get atom list
        lines = regexp(molfile,'\n','split');
        try
            atomCount = str2double(lines{4}(1:3));
        catch
            disp(lines)
        end
        if atomCount > 0
            atomBlock = lines(5:5 + atomCount - 1);
            pat = '[^a-z_A-z\*]+(?<atom>[a-z_A-Z\*]+)[^a-z_A-z\*]+';
            atoms = regexp(atomBlock,pat,'names');
            atoms = [atoms{:}]';
            atoms = {atoms.atom}';

            % If structure contains a non-elemental atom or repeat units it
            % will not be appended to SDF
            if ~all(ismember(atoms,elements)) || ~isempty(strmatch('M  STY',lines))
                metList = setdiff(metList,molfileNames{m});
                continue;
            end
        end
    end

    % Add molfile to SDF
    if first
        fprintf(fid,'%s\n',molfile);
        first = false;
    else
        fprintf(fid,'%s\n','$$$$');
        fprintf(fid,'%s\n',molfile);
    end
end
fclose(fid);

noMolFileCount = sum(~ismember(mets,molfileNames));
fprintf('Percentage of metabolites without mol files: %.1f%%\n', 100*noMolFileCount/length(mets));

metList = reshape(metList,length(metList),1);
if isnumeric(omets)
    metList = str2double(metList);
end

author: opencobra

repo: cobratoolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Released under the MIT License.
% If you use this code, please cite the following paper:
% Mahmoud Afifi, Abdelrahman Abdelhamed, Abdullah Abuolaim, Abhijith
% Punnappurath, and Michael S Brown. CIE XYZ Net: Unprocessing Images for
% Low-Level Computer Vision Tasks. arXiv preprint, 2020.
%
% Author: Mahmoud Afifi | Email: mafifi@eecs.yorku.ca, m.3afifi@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

classdef resizeLayer_ < nnet.layer.Layer

    properties
        target_size
    end

    methods
        function layer = resizeLayer_(name, target_size)
            layer.Name = name;
            layer.Description = "Image resizing";
            layer.target_size = target_size;
        end

        function Z = predict(layer, X)
            sz = size(X);
            if length(sz) < 3
                sz = [sz ones(1, 3 - length(sz))];

            end
            poolsize_d1 = floor(sz(1)/layer.target_size);
            poolsize_d2 = floor(sz(2)/layer.target_size);

            if isempty(X) || sz(1) == 0
                Z = X;
            else
                L = length(X(:))/(sz(1) * sz(2) * sz(3));
                if L == 1
                    Z = zeros(layer.target_size, layer.target_size, ...
                        sz(3), 'like',X);
                    temp = avgpool(X,[poolsize_d1, poolsize_d2], ...
                        'Stride', [poolsize_d1, poolsize_d2],...
                        'DataFormat','SSCB');
                    Z = temp(1:layer.target_size,1:layer.target_size,:);
                    size(Z)
                else
                    Z = zeros(layer.target_size, layer.target_size, ...
                        sz(3), sz(4), 'like', X);
                    for i = 1: L
                        temp = avgpool(X(:,:,:,i),...
                            [poolsize_d1, poolsize_d2],...
                            'Stride', [poolsize_d1, poolsize_d2],...
                            'DataFormat','SSCB');
                        Z(:,:,:,i) = ...
                            temp(1:layer.target_size, ...
                            1:layer.target_size,:,i);
                    end
                end
            end
        end
    end
end

author: mahmoudnafifi


text:
function [D, matches] = run_icp_fixed(M, N, C_init, max_iters)

verbose = 1;

if verbose
    fprintf('Running ICP...\n');
end

ann_params = struct;
% ann_params.algorithm = 'linear'; %use 'kmeans' for some speedup
% ann_params.trees = 8;
% ann_params.checks = 64;
% ann_params.centers_init = 'kmeanspp';
% ann_params.iterations = -1;

kN = size(N.evecs,1);

% flann_search(target, query)
[matches, dists] = flann_search(C_init*M.evecs', N.evecs', 1, ann_params);
% [matches, dists] = knnsearch(N.S*N.evecs, M.S*M.evecs*C_init');

err = sum(sqrt(dists));
err = err / (kN*size(C_init,1));

if verbose
    fprintf('(0) MSE: %.2e\n', err);
end

if max_iters == 0
    D = C_init;
%     matches = matches';
    return
end

% Start iterations

D_prev = C_init;
err_prev = err;
matches_prev = matches;

% vidx = 2:4;
% figure, plot_cloud([],N.evecs(:,vidx),'b.'); axis equal; hold on; plot_cloud([],M.evecs(:,vidx),'r.');
% figure, plot_cloud([],N.evecs(:,vidx),'b.'); axis equal; hold on; pp = M.evecs*C_init';plot_cloud([],pp(:,vidx),'r.');

% [u,~,v] = svd(C_init);
% D = u*v';

for i=1:max_iters

%     if i>1
        [U,~,V] = svd((M.evecs(matches,:)'*M.S(matches,matches)) * (N.S*N.evecs));
        D = U * V(:,1:size(C_init,2))';
        D = D';
%     end

%     figure, plot_cloud([],N.evecs(:,vidx),'b.'); axis equal; hold on; pp = M.evecs*D';plot_cloud([],pp(:,vidx),'r.');

    %     matches = flann_search(M.evecs', D*N.evecs', 1, ann_params);
    [matches, dists] = flann_search(D*M.evecs', N.evecs', 1, ann_params);
%     [matches, dists] = knnsearch(N.S*N.evecs, M.S*M.evecs*D');
    err = sum(sqrt(dists));
    err = err / (kN*size(C_init,1));

    if verbose
        fprintf('(%d) MSE: %.2e\n', i, err);
    end

    if err > err_prev
        if verbose
            fprintf('Local optimum reached.\n');
        end
        D = D_prev;
        matches = matches_prev;
        break;
    end

    if (err_prev - err) < 5e-6
        if verbose
            fprintf('Local optimum reached.\n');
        end
        break;
    end

    err_prev = err;
    D_prev = D;
    matches_prev = matches;

end

% matches = matches';

end

author: OshriHalimi


text:
% DEMROBOTWIRELESSVARGPLVMDYN1 Run variational GPLVM on robot wireless
% data.
%
% DESC: Run variational GPLVM on robot wireless data with the option to
% also add dynamics to the model.
%
% COPYRIGHT : Michalis K. Titsias, 2009-2011
% COPYRIGHT : Neil D. Lawrence, 2009-2011
% COPYRIGHT : Andreas C. Damianou, 2010-2011
%
% VARGPLVM

%clear;

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'robotWireless';
clear timeStamps; % in case it's left from a previous experiment

%-- Constants
if ~exist('trainModel'), trainModel=1; end
if ~exist('indPoints'), indPoints = 70; end % Default: 50
if ~exist('latentDim'), latentDim = 10; end
if ~exist('dynUsed'),  dynUsed = 1; end
if ~exist('printDiagram'), printDiagram = 1; end
if ~exist('experimentNo'), experimentNo = 404; end
if ~exist('doPredictions'), doPredictions = 0; end
if ~exist('itNo'), itNo = 5000; end
if ~exist('initVardistIters'), initVardistIters = 400; end
if ~exist('dynamicKern')   ,     dynamicKern = {'rbf', 'white', 'bias'}; end
if ~exist('vardistCovarsMult'), vardistCovarsMult = 2.5; end

% load data
%[Y, lbls] = lvmLoadData(dataSetName);
%%% Like lvmLoadData but also parse times.
fprintf(1,'# Preparing the dataset...\n');
baseDir = datasetsDirectory;
dirSep = filesep;
[Ydat, timeStampsdat, wireless_x, wireless_y, storedMacs] = parseWirelessData([baseDir 'uw-floor.txt']);

Ydat = (Ydat + 85)/15;
lbls = 'connect';

if ~exist('dataToKeep')
    dataToKeep = size(Ydat,1); % Default should be: 215
end

Y = Ydat(1:dataToKeep, :);
timeStampsTraining = timeStampsdat(1:dataToKeep);
Ytest = Ydat((dataToKeep+1):end, :);
timeStampsTest = timeStampsdat((dataToKeep+1):end);
%%%

% Set up model
options = vargplvmOptions('dtcvar');
options.kern = 'rbfardjit';%{'rbfard2', 'bias', 'white'};
options.numActive = indPoints; % Default: 50

options.optimiser = 'scg2';
d = size(Y, 2);

capName = dataSetName;
capName(1) = upper(capName(1));
modelType = 'Vargplvm';
saveName = ['dem' capName modelType num2str(experimentNo) '.mat'];

if trainModel
    % demo using the variational inference method for the gplvm model
    fprintf(1,'# Creating the model...\n');
    model = vargplvmCreate(latentDim, d, Y, options);
    %
    model = vargplvmParamInit(model, model.m, model.X);
    model.vardist.covars = 0.5*ones(size(model.vardist.covars)) + 0.001*randn(size(model.vardist.covars));
    if dynUsed
        optionsDyn.type = 'vargpTime';
        optionsDyn.t=timeStampsTraining;
        optionsDyn.inverseWidth=20;
        optionsDyn.kern = dynamicKern;
        optionsDyn.vardistCovars = vardistCovarsMult;

        % Fill in with default values whatever is not already set
        optionsDyn = vargplvmOptionsDyn(optionsDyn);
        model = vargplvmAddDynamics(model, 'vargpTime', optionsDyn, optionsDyn.t, 0, 0,optionsDyn.seq);

        fprintf(1,'# Further calibration of the initial parameters...\n');
        model = vargplvmInitDynamics(model,optionsDyn);
    end
    modelInit = model;

    %     for q=1:model.q, plot(X(:,q)); hold on; plot(model.vardist.means(:,q),'r');
    %         plot(model.X_u(:,q),'+r');  pause(1); hold off;
    %     end
    %

    %model.dynamics.vardist.covars = ones(size(model.dynamics.vardist.covars)); % Good initialization

    % Optimise the model.
    model.initVardist = 1; model.learnSigmaf = 0;
    display = 1;
    fprintf(1,'# Optimising the model (initialising var.distr) for %d iters...\n',initVardistIters);
    model = vargplvmOptimise(model, display, initVardistIters);

    model.initVardist = 0; model.learnSigmaf = 1;

    fprintf(1,'# Optimising the model for %d iters...\n',itNo);
    model = vargplvmOptimise(model, display, itNo);

    fprintf(1,'# Saving the model...(%s)\n',saveName);

    prunedModel = vargplvmPruneModel(model);
    save(saveName, 'prunedModel');

    % order wrt to the inputScales
    mm = vargplvmReduceModel(model,2);
    % plot the two largest twe latent dimensions
    if exist('printDiagram') & printDiagram
        lvmPrintPlot(mm, lbls, capName, experimentNo);
    end

else
    load(saveName);
    model = vargplvmRestorePrunedModel(prunedModel, Y);
end

% fprintf(1,'Plotting inferred latent GPs (intial with blue and final  with red)...\n');
% for q=1:model.q, plot(X(:,q)); hold on; plot(model.vardist.means(:,q),'r');
%     plot(model.X_u(:,q),'+r');  pause(1); hold off;
% end

%%%% Reconstruct training data %%%%%%%%
% mu  = vargplvmPosteriorMeanVar(model, model.vardist.means, model.vardist.covars);
% Varmutr=mu;
% % TODO

if ~(doPredictions && dynUsed)
    return
end

%%%% predictions
Yts = Ytest;
Varmu = [];
Varsigma = [];

fprintf(1, '# Prediction...\n');
t_star = timeStampsTest;
[x varx] = vargplvmPredictPoint(model.dynamics, t_star);

% find the largest dimensions
[max, imax] = sort(model.kern.comp{1}.inputScales,'descend');

N = size(model.X,1);
Nstar = size(Yts,1);
figure;
fprintf(1,'The two larger latent dims after re-training with partial test data. Red are the test points\n');
plot(model.X(1:model.N,imax(1)), model.X(1:model.N,imax(2)),'--rs', 'Color','b');
hold on
plot(x(:,imax(1)),x(:,imax(2)),'--rs', 'Color','r');
title('Visualization of the latent space after re-training with partial test data');
hold off

figure;

% for i=1:size(t_star,1)
%     [mu, sigma] = vargplvmPosteriorMeanVar(model, x(i,:), varx(i,:));
%     Varmu(i,:) = mu;
%     Varsigma(i,:) = sigma;
%     %
% end

[mu, sigma] = vargplvmPosteriorMeanVar(model, x, varx);
Varmu = mu;
Varsigma = sigma;

% Find the square error
errsum = sum((Varmu - Yts).^2);
% Devide by the total number of test points to find a mean
error = errsum / size(Varmu,1);

% Find the mean sq. error for each datapoint, to see how it increases as
% we give points further to the future. This time sum accros dimensions
errSq = sum((Varmu - Yts).^2 ,2);
errPt = errSq ./ size(Varmu,2);

fprintf(1,'*** Mean error: %d\n', mean(error));

% See how the latent space looks like
N = size(model.X,1);
Nstar = size(t_star,1);
newX = zeros( N+Nstar, latentDim );
newX(1:N,:) = model.X(:,:);
newX(N+1:N+Nstar, :) = x(:,:);
figure;plot(errPt);
xlabel('time test datapoints','fontsize',18);
ylabel('Mean error','fontsize',18);
%---

% Now we will plot the latent space INCLUDING the newly predicted points.
model.X = newX;
modelTEMP = vargplvmCreate(latentDim, d, Ydat, options);
model.m = modelTEMP.m;
model.y = modelTEMP.y;
model.N = modelTEMP.N;
model.vardist.means = [model.vardist.means; x];
model.vardist.covars = [model.vardist.covars; varx];
%%%%%%%%% end: predictions

% order wrt to the inputScales
mm = vargplvmReduceModel(model,2);
% plot the two largest twe latent dimensions
if exist('printDiagram') & printDiagram
    lvmPrintPlot(mm, lbls, capName, experimentNo);
    %lvmResultsDynamic(model.type, dataSetName, experimentNo, 'robotWireless','vector')
end

author: SheffieldML

repo: vargplvm