text:
classdef ProbeClass < FileLoadSaveClass

    properties
        wavelengths
        wavelengthsEmission
        sourcePos2D
        detectorPos2D
        landmarkPos2D
        landmarkPos3D
        sourcePos3D
        detectorPos3D
        frequencies
        timeDelays
        timeDelayWidths
        momentOrders
        correlationTimeDelays
        correlationTimeDelayWidths
        sourceLabels
        detectorLabels
        landmarkLabels
    end

    properties (Access = private)
        scaling
    end

    methods

        % -------------------------------------------------------
        function obj = ProbeClass(varargin)
            % Set class properties not part of the SNIRF format
            obj.SetFileFormat('hdf5');
            obj.scaling = 1;

            % Set SNIRF fomat properties
            if nargin>0
                if isstruct(varargin{1})
                    n = NirsClass(varargin{1});
                    SD = n.SD;
                    obj.wavelengths = SD.Lambda;
                    obj.wavelengthsEmission  = [];
                    if isfield(SD,'SrcPos2D') &&  ~isempty(SD.SrcPos2D)
                        obj.sourcePos2D  = SD.SrcPos2D;
                    else
                        obj.sourcePos2D  = SD.SrcPos;
                    end
                    if isfield(SD,'DetPos2D') && ~isempty(SD.DetPos2D)
                        obj.detectorPos2D  = SD.DetPos2D;
                    else
                        obj.detectorPos2D  = SD.DetPos;
                    end
                    if isfield(SD,'SrcPos3D') && ~isempty(SD.SrcPos3D)
                        obj.sourcePos3D  = SD.SrcPos3D;
                    else
                        obj.sourcePos3D  = SD.SrcPos;
                    end
                    if isfield(SD,'DetPos3D') && ~isempty(SD.DetPos3D)
                        obj.detectorPos3D  = SD.DetPos3D;
                    else
                        obj.detectorPos3D  = SD.DetPos;
                    end
                    if isfield(SD,'Landmarks3D')
                        obj.landmarkPos3D = SD.Landmarks3D.pos;
                        obj.landmarkLabels = SD.Landmarks3D.labels;
                        if isfield(SD,'Landmarks2D')
                            obj.landmarkPos2D = SD.Landmarks2D.pos;
                        end
                    else
                        obj.landmarkPos3D = SD.DummyPos;
                        obj.landmarkLabels = SD.Landmarks3D.labels;
                    end
                    obj.frequencies  = 1;
                    obj.timeDelays  = 0;
                    obj.timeDelayWidths  = 0;
                    obj.momentOrders = [];
                    obj.correlationTimeDelays = 0;
                    obj.correlationTimeDelayWidths = 0;
                    for ii=1:size(SD.SrcPos)
                        obj.sourceLabels{ii} = ['S',num2str(ii)];
                    end
                    for ii=1:size(SD.DetPos)
                        obj.detectorLabels{ii} = ['D',num2str(ii)];
                    end
                elseif ischar(varargin{1})
                    obj.SetFilename(varargin{1});
                    obj.Load(varargin{1});
                end
            else
                obj.wavelengths          = [];
                obj.wavelengthsEmission  = [];
                obj.sourcePos2D  = [];
                obj.detectorPos2D  = [];
                obj.sourcePos3D  = [];
                obj.detectorPos3D  = [];
                obj.frequencies  = 1;
                obj.timeDelays  = 0;
                obj.timeDelayWidths  = 0;
                obj.momentOrders = [];
                obj.correlationTimeDelays = 0;
                obj.correlationTimeDelayWidths = 0;
                obj.sourceLabels = {};
                obj.detectorLabels = {};
            end
        end

        % -------------------------------------------------------
        function ForwardCompatibility(obj)
            if size(obj.sourcePos2D,2)<3
                obj.sourcePos2D       = [obj.sourcePos2D, zeros(size(obj.sourcePos2D,1), 1)];
            end
            if size(obj.detectorPos2D,2)<3
                obj.detectorPos2D     = [obj.detectorPos2D, zeros(size(obj.detectorPos2D,1), 1)];
            end
        end

        % -------------------------------------------------------
        function BackwardCompatibility(obj)
            if isempty(obj.sourcePos2D)
                obj.sourcePos2D   = HDF5_DatasetLoad(gid, 'sourcePos', [], '2D');
            end
            if isempty(obj.detectorPos2D)
                obj.detectorPos2D = HDF5_DatasetLoad(gid, 'detectorPos', [], '2D');
            end
            if isempty(obj.landmarkPos2D)
                obj.landmarkPos2D = HDF5_DatasetLoad(gid, 'landmarkPos', [], '2D');
            end
        end

        % ----------------------------------------------
        function isValid = isValidLandmarkLabels(obj)
            isValid = 1;
            refpts_labels = {'T7','T8','Oz','Fpz','Cz','C3','C4','Pz','Fz'};
            for u = 1:length(refpts_labels)
                label = refpts_labels{u};
                idx = ismember(obj.landmarkLabels, label);
                if sum(idx) == 0
                    isValid = 0;
                    return
                end
            end
        end

        % -------------------------------------------------------
        function Project_3D_to_2D(obj)
            if isempty(obj.sourcePos2D) && isempty(obj.detectorPos2D)
                if isempty(obj.landmarkPos3D) || ~obj.isValidLandmarkLabels()
                    nSource = size(obj.sourcePos3D,1);
                    optodePos3D = [obj.sourcePos3D; obj.detectorPos3D];

                    optodePos2D = project_3D_to_2D(optodePos3D);
                    if ~isempty(optodePos2D)
                        if nSource ~= 0
                            obj.sourcePos2D = optodePos2D(1:nSource,:);
                        end
                        if size(optodePos3D,1) > nSource
                            obj.detectorPos2D = optodePos2D(nSource+1:end,:);
                        end
                    end
                else
                     if ~isempty(obj.sourcePos3D) && ~isempty(obj.detectorPos3D)
                        [sphere.label, sphere.theta, sphere.phi, sphere.r, sphere.xc, sphere.yc, sphere.zc] = textread('10-5-System_Mastoids_EGI129.csd','%s %f %f %f %f %f %f','commentstyle','c++');
                        % ref pt labels used for affine transformation
                        refpts_labels = {'T7','T8','Oz','Fpz','Cz','C3','C4','Pz','Fz'};

                        % get positions for refpts_labels from both sphere and probe ref pts
                        for u = 1:length(refpts_labels)
                            label = refpts_labels{u};

                            idx = ismember(obj.landmarkLabels, label);
                            if isempty(idx)
                                return
                            end
                            probe_refps_pos(u,:) = obj.landmarkPos3D(idx,:);
                            idx = ismember(sphere.label, label);
                            sphere_refpts_pos(u,:) = [sphere.xc(idx) sphere.yc(idx) sphere.zc(idx)];
                        end

                        % get affine transformation
                        % probe_refps*T = sphere_refpts
                        probe_refps_pos = [probe_refps_pos(:,1:3) ones(size(probe_refps_pos,1),1)];
                        T = probe_refps_pos\sphere_refpts_pos;

                        % tranform optode positions onto unit sphere.
                        % opt_pos = probe.optpos_reg;
                        % opt_pos = [opt_pos ones(size(opt_pos,1),1)];
                        % sphere_opt_pos = opt_pos*T;
                        % sphere_opt_pos_norm = sqrt(sum(sphere_opt_pos.^2,2));
                        % sphere_opt_pos = sphere_opt_pos./sphere_opt_pos_norm ;
                        %
                        % get 2D circular refpts for current selecetd reference point system
                        probe_refpts_idx =  ismember(sphere.label, obj.landmarkLabels);

                        % refpts_2D.pos = [sphere_xc(probe_refpts_idx) sphere_yc(probe_refpts_idx) sphere_zc(probe_refpts_idx)];
                        refpts_2D.label = sphere.label(probe_refpts_idx);
                        %
                        refpts_theta =  sphere.theta(probe_refpts_idx);
                        refpts_phi = 90 - sphere.phi(probe_refpts_idx); % elevation angle from top axis

                        refpts_theta = (2 * pi * refpts_theta) / 360; % convert to radians
                        refpts_phi = (2 * pi * refpts_phi) / 360;
                        [x,y] = pol2cart(refpts_theta, refpts_phi);      % get plane coordinates
                        xy = [x y];

                        %
                        norm_factor = max(max(xy));
                        xy = xy/norm_factor;               % set maximum to unit length
                        refpts_2D.pos = xy;
                        obj.landmarkPos2D = refpts_2D.pos;

                        %
                        obj.sourcePos2D = convert_optodepos_to_circlular_2D_pos(obj.sourcePos3D, T, norm_factor);
                        obj.detectorPos2D = convert_optodepos_to_circlular_2D_pos(obj.detectorPos3D, T, norm_factor);
                    end
                end
            end
        end

        % -------------------------------------------------------
        function err = LoadHdf5(obj, fileobj, location, LengthUnit)
            err = 0;

            % Arg 1
            if ~exist('fileobj','var') || (ischar(fileobj) && ~exist(fileobj,'file'))
                fileobj = '';
            end

            % Arg 2
            if ~exist('location', 'var') || isempty(location)
                location = '/nirs/probe';
            elseif location(1)~='/'
                location = ['/',location];
            end

            % Arg3
            if exist('LengthUnit','var')
                % Figure out the scaling factor to multiple by to get the coorinates to be in mm units
                if strcmpi(LengthUnit,'m')  % meter units
                    obj.scaling = 1000;
                elseif strcmpi(LengthUnit,'cm')  % centimeter units
                    obj.scaling = 10;
                end
            end

            % Error checking
            if ~isempty(fileobj) && ischar(fileobj)
                obj.SetFilename(fileobj);
            elseif isempty(fileobj)
                fileobj = obj.GetFilename();
            end
            if isempty(fileobj)
               err = -1;
               return;
            end

            try
                % Open group
                [gid, fid] = HDF5_GroupOpen(fileobj, location);

                % Load datasets
                obj.wavelengths               = HDF5_DatasetLoad(gid, 'wavelengths');
                obj.wavelengthsEmission       = HDF5_DatasetLoad(gid, 'wavelengthsEmission');
                obj.sourcePos2D               = HDF5_DatasetLoad(gid, 'sourcePos2D', [], '2D') * obj.scaling;
                obj.detectorPos2D             = HDF5_DatasetLoad(gid, 'detectorPos2D', [], '2D') * obj.scaling;
                obj.landmarkPos2D             = HDF5_DatasetLoad(gid, 'landmarkPos2D', [], '2D') * obj.scaling;
                obj.sourcePos3D               = HDF5_DatasetLoad(gid, 'sourcePos3D', [], '3D') * obj.scaling;
                obj.detectorPos3D             = HDF5_DatasetLoad(gid, 'detectorPos3D', [], '3D') * obj.scaling;
                obj.landmarkPos3D             = HDF5_DatasetLoad(gid, 'landmarkPos3D', [], '2D') * obj.scaling;
                obj.frequencies               = HDF5_DatasetLoad(gid, 'frequencies');
                obj.timeDelays                 = HDF5_DatasetLoad(gid, 'timeDelays');
                obj.timeDelayWidths            = HDF5_DatasetLoad(gid, 'timeDelayWidths');
                obj.momentOrders               = HDF5_DatasetLoad(gid, 'momentOrders');
                obj.correlationTimeDelays      = HDF5_DatasetLoad(gid, 'correlationTimeDelays');
                obj.correlationTimeDelayWidths = HDF5_DatasetLoad(gid, 'correlationTimeDelayWidths');
                obj.sourceLabels              = HDF5_DatasetLoad(gid, 'sourceLabels', obj.sourceLabels);
                obj.detectorLabels            = HDF5_DatasetLoad(gid, 'detectorLabels', obj.detectorLabels);
                obj.landmarkLabels            = HDF5_DatasetLoad(gid, 'landmarkLabels', obj.landmarkLabels);

                obj.Project_3D_to_2D();

                % Close group
                HDF5_GroupClose(fileobj, gid, fid);

                assert(obj.IsValid())

            catch
                err = -1;
                return;
            end

            % Call method to change future current and future versions of
            % SNIRF data to Homer3 compatible structure
            obj.ForwardCompatibility();

            % for Homer3 usage, add 3D positions if they are empty
            if isempty(obj.sourcePos3D)
                obj.sourcePos3D = obj.sourcePos2D;
                if size(obj.sourcePos3D,2)<3
                    obj.sourcePos3D(:,3) = 0;
                end
            end

            if isempty(obj.detectorPos3D)
                obj.detectorPos3D = obj.detectorPos2D;
                if size(obj.detectorPos3D,2)<3
                    obj.detectorPos3D(:,3) = 0;
                end
            end
            obj.SetError(err);
        end

        % -------------------------------------------------------
        function SaveHdf5(obj, fileobj, location)
            % Arg 1
            if ~exist('fileobj', 'var') || isempty(fileobj)
                error('Unable to save file. No file name given.')
            end

            % Arg 2
            if ~exist('location', 'var') || isempty(location)
                location = '/nirs/probe';
            elseif location(1)~='/'
                location = ['/',location];
            end

            if ~exist(fileobj, 'file')
                fid = H5F.create(fileobj, 'H5F_ACC_TRUNC', 'H5P_DEFAULT', 'H5P_DEFAULT');
                H5F.close(fid);
            end

            % Multiple all coordinates by reciprocal of scaling factor when saving to get the
            % coordinates back to original.
            obj.sourcePos2D    = obj.sourcePos2D / obj.scaling;
            obj.detectorPos2D  = obj.detectorPos2D / obj.scaling;
            obj.landmarkPos2D  = obj.landmarkPos2D / obj.scaling;
            obj.sourcePos3D    = obj.sourcePos3D / obj.scaling;
            obj.detectorPos3D  = obj.detectorPos3D / obj.scaling;
            obj.landmarkPos3D  = obj.landmarkPos3D / obj.scaling;

            % Now save
            hdf5write_safe(fileobj, [location, '/wavelengths'], obj.wavelengths, 'array');
            hdf5write_safe(fileobj, [location, '/wavelengthsEmission'], obj.wavelengthsEmission, 'array');
            hdf5write_safe(fileobj, [location, '/sourcePos2D'], obj.sourcePos2D, 'array');
            hdf5write_safe(fileobj, [location, '/detectorPos2D'], obj.detectorPos2D, 'array');
            hdf5write_safe(fileobj, [location, '/landmarkPos2D'], obj.landmarkPos2D, 'array');
            hdf5write_safe(fileobj, [location, '/sourcePos3D'], obj.sourcePos3D, 'array');
            hdf5write_safe(fileobj, [location, '/detectorPos3D'], obj.detectorPos3D, 'array');
            hdf5write_safe(fileobj, [location, '/landmarkPos3D'], obj.landmarkPos3D, 'array');
            hdf5write_safe(fileobj, [location, '/frequencies'], obj.frequencies, 'array');
            hdf5write_safe(fileobj, [location, '/timeDelays'], obj.timeDelays, 'array');
            hdf5write_safe(fileobj, [location, '/timeDelayWidths'], obj.timeDelayWidths, 'array');
            hdf5write_safe(fileobj, [location, '/momentOrders'], obj.momentOrders, 'array');
            hdf5write_safe(fileobj, [location, '/correlationTimeDelays'], obj.correlationTimeDelays, 'array');
            hdf5write_safe(fileobj, [location, '/correlationTimeDelayWidths'], obj.correlationTimeDelayWidths, 'array');
            hdf5write_safe(fileobj, [location, '/sourceLabels'], obj.sourceLabels, 'array');
            hdf5write_safe(fileobj, [location, '/detectorLabels'], obj.detectorLabels, 'array');
            hdf5write_safe(fileobj, [location, '/landmarkLabels'], obj.landmarkLabels, 'array');
        end

        % ---------------------------------------------------------
        function wls = GetWls(obj)
            wls = obj.wavelengths;
        end

        % ---------------------------------------------------------
        function srcpos = GetSrcPos(obj,flag2d)
            if ~exist('flag2d','var')
                flag2d = 0;
            end
            if flag2d==0
                srcpos = obj.sourcePos3D;
            else
                srcpos = obj.sourcePos2D;
            end
        end

        % ---------------------------------------------------------
        function detpos = GetDetPos(obj,flag2d)
            if ~exist('flag2d','var')
                flag2d = 0;
            end
            if flag2d==0
                detpos = obj.detectorPos3D;
            else
                detpos = obj.detectorPos2D;
            end
        end

        % -------------------------------------------------------
        function B = eq(obj, obj2)
            B = false;
            if isempty(obj) && ~isempty(obj2)
                return;
            end
            if isempty(obj) && ~isempty(obj2)
                return;
            end
            if ~isempty(obj) && isempty(obj2)
                return;
            end
            if isempty(obj) && isempty(obj2)
                b = true;
                return;
            end
            if ~all(obj.wavelengths(:)==obj2.wavelengths(:))
                return;
            end
            if ~all(obj.wavelengthsEmission(:)==obj2.wavelengthsEmission(:))
                return;
            end
            if ~all(obj.sourcePos2D(:)==obj2.sourcePos2D(:))
                return;
            end
            if ~all(obj.detectorPos2D(:)==obj2.detectorPos2D(:))
                return;
            end
            if ~all(obj.sourcePos3D(:)==obj2.sourcePos3D(:))
                return;
            end
            if ~all(obj.detectorPos3D(:)==obj2.detectorPos3D(:))
                return;
            end
            if ~all(obj.frequencies(:)==obj2.frequencies(:))
                return;
            end
            if ~all(obj.timeDelays(:)==obj2.timeDelays(:))
                return;
            end
            if ~all(obj.timeDelayWidths(:)==obj2.timeDelayWidths(:))
                return;
            end
            if ~all(obj.momentOrders(:)==obj2.momentOrders(:))
                return;
            end
            if ~all(obj.correlationTimeDelays(:)==obj2.correlationTimeDelays(:))
                return;
            end
            if ~all(obj.correlationTimeDelayWidths(:)==obj2.correlationTimeDelayWidths(:))
                return;
            end
            if length(obj.sourceLabels)~=length(obj2.sourceLabels)
                return;
            end
            for ii=1:length(obj.sourceLabels)
                if ~strcmp(obj.sourceLabels{ii}, obj2.sourceLabels{ii})
                    return;
                end
            end
            if length(obj.detectorLabels)~=length(obj2.detectorLabels)
                return;
            end
            for ii=1:length(obj.detectorLabels)
                if ~strcmp(obj.detectorLabels{ii}, obj2.detectorLabels{ii})
                    return;
                end
            end
            B = true;
        end

        % ----------------------------------------------------------------------------------
        function nbytes = MemoryRequired(obj)
            nbytes = 0;
            fields = properties(obj);
            for ii=1:length(fields)
                nbytes = nbytes + eval(sprintf('sizeof(obj.%s)', fields{ii}));
            end
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmpty(obj)
            b = true;
            if isempty(obj.wavelengths)
                return;
            end
            if isempty(obj.sourcePos2D) && isempty(obj.detectorPos2D) && ...
                    isempty(obj.sourcePos3D) && isempty(obj.detectorPos3D)
                return;
            end
            b = false;
        end

        % ----------------------------------------------------------------------------------
        function b = IsValid(obj)
            b = false;
            if obj.IsEmpty()
                return;
            end
            if iscolumn(obj.sourcePos2D)
                return;
            end
            if length(obj.sourcePos2D)>4
                if size(obj.sourcePos2D,2) > size(obj.sourcePos2D,1)
                    return;
                end
            end
            if iscolumn(obj.detectorPos2D)
                return;
            end
            if length(obj.detectorPos2D)>4
                if size(obj.detectorPos2D,2) > size(obj.detectorPos2D,1)
                    return;
                end
            end
            b = true;
        end

        % ---------------------------------------------------------
        function bbox = GetSdgBbox(obj)
            bbox = [];

            optpos = [obj.GetSrcPos('2D'); obj.GetDetPos('2D')];
            if isempty(optpos)
                return
            end

            xmax = max(optpos(:,1));
            ymax = max(optpos(:,2));

            xmin = min(optpos(:,1));
            ymin = min(optpos(:,2));

            width = xmax-xmin;
            height = ymax-ymin;

            if width==0
                width = 1;
            end
            if height==0
                height = 1;
            end

            px = width * 0.05;
            py = height * 0.05;

            bbox = [xmin-px, xmax+px, ymin-py, ymax+py];
        end

        % ---------------------------------------------------------
        function val = GetScaleFactor(obj)
            val = obj.scaling;
        end

    end

end

author: BUNPC


text:
function imgCompleteFinal = sc_complete(imgFileName)
% SC_COMPLETE: image completion with user-defined mask
%
% Input:
%     - imgFileName: image file name
% Output:
%     - imgCompletion: completed result
%
% Example:
%   imgFileName   = '002_input_hole.png';
%   imgCompletion = sc_complete(imgFileName);
%
% Disclaimer:
%   This is a Matlab re-implementation of the paper:
%
%   Jia-Bin Huang, Sing Bing Kang, Narendra Ahuja, and Johannes Kopf,
%   Image Completion using Planar Structure Guidance,
%   ACM Transactions on Graphics (Proceedings of SIGGRAPH 2014), 33(4), 2014
%
%   It is provided for educational/researrch purpose only. If you find the
%   software useful, please consider cite our paper.
%
% Contact:
%   Jia-Bin Huang
%   University of Illinois, Urbana-Champaign
%   www.jiabinhuang.com
<<<<<<< HEAD
%   jbhuang0604@gmail.com

% Set up required path
startup;

% Load image
imgID = 9;
imgFileName = [num2str(imgID, '%03d'), '_input_hole.png'];
=======
%   jbhuang0604@gmail.com

% % Load image
% imgID = 3;
% imgFileName = [num2str(imgID, '%03d'), '_input_hole.png'];
>>>>>>> 80fe71bd858652004119dc8946439915385cc282

imgFileName = 'test.png';

% Option parameters
[optA, optS] = sc_init_opt;

% =========================================================================
% Planar structure extraction
% =========================================================================
fprintf('- Extract planar structures \n');
tic;
[img, mask, maskD, modelPlane, modelReg] = sc_extract_planar_structure(imgFileName, optA);
tAnalysis = toc;
fprintf('Done in %6.3f seconds.\n\n', tAnalysis);

% =========================================================================
% Guided image completion
% =========================================================================
% Construct image pyramid for coarse-to-fine image completion
fprintf('- Construct image pyramid: \n');
tic;
% Create image pyramid
[imgPyr, maskPyr, scaleImgPyr] = sc_create_pyramid(img, maskD, optS);
% Structure constraints
[modelPlane, modelReg] = sc_planar_structure_pyramid(scaleImgPyr, modelPlane, modelReg);
tImgPyramid = toc;
fprintf('Done in %6.3f seconds.\n\n', tImgPyramid);

% Completion by synthesis
fprintf('- Image completion using planar structure guidance \n');
tic;
imgPyr = sc_synthesis(imgPyr, maskPyr, modelPlane, modelReg, optS);
tSynthesis = toc;
fprintf('Synthesis took %6.3f seconds.\n', tSynthesis);

% Return the top level
imgSyn = imgPyr{optS.topLevel};

% =========================================================================
% Poisson blending
% =========================================================================
tic;
imgCompleteFinal = sc_poisson_blend(img, imgSyn, mask);
tBlend = toc;
fprintf('Blending took %6.3f seconds.\n', tBlend);

imwrite(imgCompleteFinal, fullfile('result', [imgFileName(1:end-4), '_completion.png']));

% [To-Do] Visualizing completion process
% sc_vis_synthesis_process(imgCompletePyr, imgPyrNNF, imgPyrCost);

end

repo: StructCompletion


text:
% this function calculates the required traffic inforamtion for a basic cell in our example the ju area

function traffic()
% here we define oure figure from pushbuttons to edit boxes to axis and etc..
figure('name','traffic','menubar','none','numbertitle','off','units','normalized','color',[1 1 1],'position',[0.25 0.25 0.7 0.6])
uicontrol('style','text','units','normalized','position',[0.05 0.8 0.2 0.1],'string','Number of users');
e1=uicontrol('style','edit','units','normalized','position',[0.25 0.8 0.2 0.1]);
uicontrol('style','text','units','normalized','position',[0.05 0.6 0.2 0.1],'string','Lambda');
e2=uicontrol('style','edit','units','normalized','position',[0.25 0.6 0.2 0.1]);
uicontrol('style','text','units','normalized','position',[0.05 0.5 0.2 0.1],'string','Gamma');
e3=uicontrol('style','edit','units','normalized','position',[0.25 0.5 0.2 0.1]);
uicontrol('style','text','units','normalized','position',[0.05 0.4 0.2 0.1],'string','Totla number of channels');
e4=uicontrol('style','edit','units','normalized','position',[0.25 0.4 0.2 0.1]);
uicontrol('style','text','units','normalized','position',[0.05 0.7 0.2 0.1],'string','percent');
e5=uicontrol('style','edit','units','normalized','position',[0.25 0.7 0.2 0.1]);
uicontrol('style','text','units','normalized','position',[0.05 0.3 0.2 0.1],'string','min C/I');

e6=uicontrol('style','edit','units','normalized','position',[0.25 0.3 0.2 0.1]);
uicontrol('style','text','units','normalized','position',[0.05 0.2 0.2 0.1],'string','Number of Sectors');

e7=uicontrol('style','edit','units','normalized','position',[0.25 0.2 0.2 0.1]);

uicontrol('style','pushbutton','callback',@push,'string','Calculate','units','normalized','position',[0.15 0.01 0.2 0.1]);

tt1=uicontrol('style','text','units','normalized','position',[0.5 0.7 0.4 0.1]);
tt=uicontrol('style','text','units','normalized','position',[0.46 .7 0.03 0.1],'string','A');
tt2=uicontrol('style','text','units','normalized','position',[0.5 0.6 0.4 0.1]);
uicontrol('style','text','units','normalized','position',[0.46 .6 0.03 0.1],'string','Ncell');
tt3=uicontrol('style','text','units','normalized','position',[0.5 0.5 0.4 0.1]);
uicontrol('style','text','units','normalized','position',[0.46 .5 0.03 0.1],'string','#of cell');
mi=	uimenu('label','example','callback',@ex);
	function ex(varargin)
% 		This function loads the already saved example in the directory
		load('example.mat','n','l','g','Nt','pe','ci','nm');
		n;
		l;
		g;
		Nt;
		pe;
		ci;
		nm;
		set(e1,'string',num2str(n));
		set(e2,'string',num2str(l));
		set(e3,'string',num2str(g));
		set(e4,'string',num2str(Nt));
		set(e5,'string',num2str(pe));
		set(e6,'string',num2str(ci));
		set(e7,'string',num2str(nm));
	end

	function push(varargin)
		%calculation of the traffic
		n=str2num(get(e1,'string'));
		l=str2num(get(e2,'string'));
		g=str2num(get(e3,'string'));
		Nt=str2num(get(e4,'string'));
		pe=str2num(get(e5,'string'));
		ci=str2num(get(e6,'string'));
		nm=str2num(get(e7,'string'));
		A=(n*(pe/100)*l)/60;
		set(tt1,'string',num2str(A,6));
		q=((10^(ci/10))*(6/nm))^(1/3);
		er=[0.01;0.15;0.43;0.81;1.26;1.76;2.30;2.87;3.46;4.08;4.71;5.36;6.03;6.71;7.39;8.09;8.80;9.52;10.24;20.97;11.71;12.45;13.21;12.96;14.72;15.49;16.25;17.03;17.81;18.59;19.37;20.16;20.95;21.75;22.55;23.35;24.15;24.96;25.76;26.58;27.39;28.20;29.02;29.84;30.67;31.49;32.31;33.14;33.97;34.81;35.58;36.42;37.27;38.10;38.94;39.78;40.63;41.47;42.32;43.16];
		cc=length(er);
		gg=round(er);
		k=(q^2)/3;
		cd=k;
		i=1;
		j=1;
		k=i^2+j^2+i*j;

		ee=0;
		i=0;
		j=1;
		while ee<cd

			i=i+1;
			ee=i^2+j^2+i*j;

			j=j+1;;
			k=ee;

		end
		if nm == 1
			nc=Nt/(k*nm);
			set(tt2,'string',num2str(nc));
			aaa=er(nc);
			nce=A/aaa;
			nce=ceil(nce);
			set(tt3,'string',num2str(nce));
		else
			Ns=63/(nm*k);
			As=er(round(Ns));
			Acell=3*As;
			set(tt2,'string',num2str(find(gg==ceil(Acell),1)));
			nce=A/Acell;
			nce=ceil(nce);
			set(tt3,'string',num2str(nce))
				save('model.mat','nce')
		end

	end

function push1(varargin)
	%save the number of channels that we calculated to be used on the site
	%configruation.
		save('model.mat','nce')
	end
end

author: Sable

repo: mcbench-benchmarks


text:
function [phi,dphi]=pt2BarycentricTriangCoords(x,v1,v2,v3)
%%PT2BARYCENTRICTIANGCOORDS Convert a point into barycentric points with
%           respect to the vertices of a triangle in 2D. This solves for a
%           specific set of phi such that
%           x=phi(1)*v1+phi(2)*v2+phi(3)*v3
%           and sum(phi)=1. If x is inside the triangle, then phi>=0.
%
%INPUTS: x The 2X1 point whose barycentric coordiantes are desired.
% v1,v2,v3 The three 2X1 vertices of the triangle. It doens't matter if
%          they are in clockwise or counterclockwise order.
%
%OUTPUTS: phi The 3X1 set of barycentric coordinates of x.
%        dPhi A 3X2 set of partial derivatives of phi (rows) with respect
%             to the components of x (columns).
%
%The barycentric coordinate system on a triangle is described in Section 2
%of [1].
%
%EXAMPLE:
%In this example, coordinates inside a triangle are obtained. It is
%verified that the returned derivatives are (within finite precision
%limits) within the bounds of what one gets with finite diferencing. Also,
%the normalization of the weights is verified and the relative error of the
%interpolated values is computed and is also negligible within finite
%precision limits.
% v1=[-2;3];
% v2=[11;0];
% v3=[1;1];
% x=[4;1.5];
% [phi,dphi]=pt2BarycentricTriangCoords(x,v1,v2,v3);
%
% f=@(y)pt2BarycentricTriangCoords(y,v1,v2,v3);
% dPhiNumDiff=numDiff(x,f,3);
%
% RelDerivError=max(max(abs((dPhiNumDiff-dphi)./dPhiNumDiff)))
% normalizationError=sum(phi)-1
%
% interpPoint=phi(1)*v1+phi(2)*v2+phi(3)*v3;
% RelErrVal=max(abs((interpPoint-x)./x))
%
%REFERENCES:
%[1] M. S. Floater, "Generalized barycentric coordinates and applications,"
%    Acta Numerica, vol. 24, pp. 161-214, 1 May 2015.
%
%September 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

Av=triangleArea(v1,v2,v3);

phi=zeros(3,1);
phi(1)=triangleArea(x,v2,v3)/Av;
phi(2)=triangleArea(x,v3,v1)/Av;
phi(3)=triangleArea(x,v1,v2)/Av;

if(nargout>1)
    dphi=zeros(3,2);

    dphi(1,:)=[(v2(2)-v3(2)),(v3(1)-v2(1))]/(2*Av);
    dphi(2,:)=[(v3(2)-v1(2)),(v1(1)-v3(1))]/(2*Av);
    dphi(3,:)=[(v1(2)-v2(2)),(v2(1)-v1(1))]/(2*Av);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%% CUBEPOISSONQ1 Poisson equation in a CUBE domain using trilinear cube element
%
%   cubePoissonQ1 computes trilinear finite element approximations of the
%   Poisson equation in the unit cube on a sequence of hex meshes obtained by
%   uniform refinement. It plots the approximation err vs the number of
%   degree of freedoms.
%
% See also
%   cubePoisson, cubePoissonP2
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all;

%% Parameters
maxIt = 4;
N = zeros(maxIt,1);
h = zeros(maxIt,1);

%% Domain and pde
h0 = 1.0/2.0;
cube = [0 1 0 1 0 1];

pde = sincosdata3;
option.solver = 'amg';
% option.isoparametric = true;

%% Finite Element Method
errL2 = zeros(maxIt,1);
errH1 = zeros(maxIt,1);
erruIuh = zeros(maxIt,1);
for k = 1:maxIt
    [node,elem] = cubehexmesh(cube,h0/(2.0^k));
    bdFlag = setboundary3(node,elem,'Dirichlet','abs(z)>eps','Neumann','abs(z)<eps');
%     bdFlag = setboundary3(node,elem,'Dirichlet');
%     [uh,Du,eqn,info] = Poisson3Q1(node,elem,pde,bdFlag,option);
    [uh,Du,eqn,info] = Poisson3T1(node,elem,pde,bdFlag,option);
    N(k) = size(node,1);
    h(k) = 1./(size(node,1)^(1/3)-1);
    uI = pde.exactu(node);
    erruIuh(k) = sqrt((uh-uI)'*eqn.A*(uh-uI));
    errH1(k) = getH1error3Q1(node,elem,pde.Du,uh);
    errL2(k) = getL2error3Q1(node,elem,pde.exactu,uh);
end

%% Plot convergence rates
figure;
showrateh3(h,errH1,2,'-*', '$|| Du - Du_h ||', ...
           h,errL2,2,'k-+', '|| u - u_h ||', ...
           h,erruIuh,2,'m-+','|| D u_I - D u_h ||');

fprintf('\n');
disp('Table: Error')
colname = {'#Dof','h','|| u-u_h ||','|| Du-Du_h ||','|| Du_I-Du_h ||'};
disptable(colname,N,[],h,'%0.3e',errL2,'%0.5e',errH1,'%0.5e',erruIuh,'%0.5e');

repo: ifem


text:
% Pre-Processing of Code Generation with Simulink Data Object
% Make signal and parameter objects to define variable specifications (RAM/ROM)
% Resolve signal resolution with DISABLEIMPLICITSIGNALRESOLUTION

% make signal objects
sig.name = { ...
	'timer_en1_0', 'timer_en2_0', 'timer_pen_0', ...		% signals
	'operation_mode', 'motor_number', 'pen_idx', ...		% states (data store memories)
	'theta1m_ref_bl', 'theta2m_ref_bl', 'timer_trj', ...
	'ud_backlash1', 'ud_backlash2', ...						% states (unit delay blocks)
	'ud_engage1', 'ud_engage2', ...
	'ud_flag_enter', 'ud_flag_run', 'ud_flag_pen', ...
	'ud_theta1m_ref', 'ud_theta2m_ref' ...
	};
num_sig = length(sig.name);
for n = 1:num_sig
	cmd = strcat( ...
		sig.name{n}, ' = ', 'mpt.Signal;', ...
		sig.name{n}, '.RTWInfo.CustomAttributes.DefinitionFile = ''nxtscara_ram.c'';', ...
		sig.name{n}, '.RTWInfo.CustomAttributes.HeaderFile = ''nxtscara_ram.h'';' ...
		);
	eval(cmd)
end

% resolve signal resolution
sig.handle = disableimplicitsignalresolution(bdroot(gcb));

% make parameter objects
param.name = { ...
	'backlash1', 'backlash2', 'dthetam_bl', ...
	'engage1_iv', 'engage2_iv', 'num_mode', 'num_motor', ...
	'pwm1_adjst', 'pwm2_adjst', 'pwm3_adjst', 'pwm3_pen', ...
	'pwm1_gain', 'pwm1_offset', 'pwm2_gain', 'pwm2_offset', ...
	'theta1m_max', 'theta2m_max', 'theta1m_iv', 'theta2m_iv', ...
	'time_en1', 'time_en2', 'time_finish', ...
	'time_pwm3', 'time_pen', 'time_table_pen', 'time_wait', ...
	'x_ref', 'y_ref',  ...
	};
num_param = length(param.name);
param.value = cell(1, num_param);
for n = 1:num_param
	cmd = strcat( ...
		'param.value{', num2str(n), '} = ', param.name{n}, ';', ...
		param.name{n}, ' = ', 'mpt.Parameter;', ...
		param.name{n}, '.Value = ', 'param.value{', num2str(n), '};', ...
		param.name{n}, '.RTWInfo.CustomAttributes.DefinitionFile = ''nxtscara_rom.c'';', ...
		param.name{n}, '.RTWInfo.CustomAttributes.HeaderFile = ''nxtscara_rom.h'';' ...
		);
	eval(cmd)
end


text:
function [k, sk] = linKernCompute(kern, x, x2)

% LINKERNCOMPUTE Compute the LIN kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the linear
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN k : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the linear
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x : input data matrix in the form of a design matrix.
% RETURN k : the kernel matrix computed at the given points.
%
% SEEALSO : linKernParamInit, kernCompute, kernCreate, linKernDiagCompute
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2009

% KERN

if nargin < 3
  sk = x*x';
else
  sk = x*x2';
end
k = sk*kern.variance;
if issparse(x)
  k = full(k);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout = deletechanGUI(varargin)

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @deletechanGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @deletechanGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before deletechanGUI is made visible.
function deletechanGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to deletechanGUI (see VARARGIN)

% Choose default command line output for deletechanGUI
handles.output = [];

try
typedata = varargin{1};
catch
    typedata = 'EEG';
end

%
% Name & version
%
version = geterplabversion;
set(handles.gui_chassis,'Name', ['ERPLAB ' version '   -   Remove Channel(s) GUI for ' typedata])

%
% Color GUI
%
handles = painterplab(handles);

%
% Set font size
%
handles = setfonterplab(handles);

% Update handles structure
guidata(hObject, handles);

drawnow
uiwait(handles.gui_chassis);

% --- Outputs from this function are returned to the command line.
function varargout = deletechanGUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

% The figure can be deleted now
delete(handles.gui_chassis);
pause(0.1)

%--------------------------------------------------------------------------
function edit_chanindex_Callback(hObject, eventdata, handles)

%--------------------------------------------------------------------------
function edit_chanindex_CreateFcn(hObject, eventdata, handles)

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%--------------------------------------------------------------------------
function pushbutton_cancel_Callback(hObject, eventdata, handles)

handles.output = [];
% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

%--------------------------------------------------------------------------
function pushbutton_ok_Callback(hObject, eventdata, handles)

chanArray = str2num(get(handles.edit_chanindex,'String'));
handles.output = {chanArray}; % sent like a cell string (with formulas)

% Update handles structure
guidata(hObject, handles);
uiresume(handles.gui_chassis);

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function [D,C]=JDDLDR_UDC(X,trls,lambda_a,lambda_b,D,C,Max_iteration)

nClass = max(trls);
Jstep_T = 1e-3;

disp_cycle = floor(nClass/10);
if disp_cycle < 1
    disp_cycle = 1;
end

fprintf('Update of D and C: class ');
for ci = 1:nClass
    cdat = X(:,trls==ci);
    temD = D(ci).M;
    iteration = 1;
    while  iteration < Max_iteration
    % update C
    afa = inv(temD'*temD+lambda_a*eye(size(temD,2)))*(temD'*cdat);
    PinvD = inv(temD'*temD+(lambda_a+lambda_b)*eye(size(temD,2)));
    for tj = 1:10
%         plot(afa(:,1));title(num2str(tj));pause(1);
        avg_afa = repmat(mean(afa,2),[1 size(afa,2)]);
        afa = PinvD*(temD'*cdat+lambda_b*avg_afa);
    end
    C(ci).M = afa;

    % update D
    for i=1:size(temD,2)
        ai        =    afa(i,:);
        Y         =    cdat-temD*afa+temD(:,i)*ai;
        di        =    Y*ai';
        di        =    di./norm(di,2);
        temD(:,i)    =    di;
    end
    D(ci).M  = temD;

    zz            =    cdat-temD*afa;
    zalpha        =    afa(:);
    avg_afa       =    repmat(mean(C(ci).M,2),[1 size(C(ci).M,2)]);
    z_afa         =    afa - avg_afa;
    Jnow          =    zz(:)'*zz(:)+lambda_a*sum(zalpha(:).*zalpha(:))+lambda_b*sum(z_afa(:).*z_afa(:));
    iteration     =    iteration+1;
    end

    if mod(ci, disp_cycle) == 0
        fprintf('%03d ', ci);
    end
end
fprintf('\n');

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
function [ s ] = xml2struct( file )
% Convert xml file into a MATLAB structure
%  [ s ] = xml2struct( file )
%
%  A file containing:
%  <XMLname attrib1="Some value">
%    <Element>Some text</Element>
%    <DifferentElement attrib2="2">Some more text</Element>
%    <DifferentElement attrib3="2" attrib4="1">Even more text</DifferentElement>
%  </XMLname>
%
%  Will produce:
%  s.XMLname.Attributes.attrib1 = "Some value";
%  s.XMLname.Element.Text = "Some text";
%  s.XMLname.DifferentElement{1}.Attributes.attrib2 = "2";
%  s.XMLname.DifferentElement{1}.Text = "Some more text";
%  s.XMLname.DifferentElement{2}.Attributes.attrib3 = "2";
%  s.XMLname.DifferentElement{2}.Attributes.attrib4 = "1";
%  s.XMLname.DifferentElement{2}.Text = "Even more text";
%
%  Please note that the following characters are substituted
%  '-' by '_dash_', ':' by '_colon_' and '.' by '_dot_'
%
%  Written by W. Falkena, ASTI, TUDelft, 21-08-2010
%  Attribute parsing speed increased by 40%  by A. Wanner, 14-6-2011
%  Added CDATA support by I. Smirnov, 20-3-2012
%
%  Modified by X. Mo, University of Wisconsin, 12-5-2012

    if (nargin < 1)
        clc;
        help xml2struct
        return
    end

    if isa(file, 'org.apache.xerces.dom.DeferredDocumentImpl') || isa(file, 'org.apache.xerces.dom.DeferredElementImpl')
        %  input is a java xml object
        xDoc = file;
    else
        % check for existance
        if (exist(file,'file') == 0)
            % Perhaps the xml extension was omitted from the file name. Add the
            % extension and try again.
            if (isempty(strfind(file,'.xml')))
                file = [file '.xml'];
            end

            if (exist(file,'file') == 0)
                error(['The file ' file ' could not be found']);
            end
        end
        % read the xml file
        xDoc = xmlread(file);
    end

    % parse xDoc into a MATLAB structure
    s = parseChildNodes(xDoc);

end

%  ----- Subfunction parseChildNodes -----
function [children,ptext,textflag] = parseChildNodes(theNode)
    %  Recurse over node children.
    children = struct;
    ptext = struct; textflag = 'Text';
    if hasChildNodes(theNode)
        childNodes = getChildNodes(theNode);
        numChildNodes = getLength(childNodes);

        for count = 1:numChildNodes
            theChild = item(childNodes,count-1);
            [text,name,attr,childs,textflag] = getNodeData(theChild);

            if (~strcmp(name,'#text') && ~strcmp(name,'#comment') && ~strcmp(name,'#cdata_dash_section'))
                % XML allows the same elements to be defined multiple times,
                % put each in a different cell
                if (isfield(children,name))
                    if (~iscell(children.(name)))
                        % put existsing element into cell format
                        children.(name) = {children.(name)};
                    end
                    index = length(children.(name))+1;
                    % add new element
                    children.(name){index} = childs;
                    if(~isempty(fieldnames(text)))
                        children.(name){index} = text;
                    end
                    if(~isempty(attr))
                        children.(name){index}.('Attributes') = attr;
                    end
                else
                    % add previously unknown (new) element to the structure
                    children.(name) = childs;
                    if(~isempty(text) && ~isempty(fieldnames(text)))
                        children.(name) = text;
                    end
                    if(~isempty(attr))
                        children.(name).('Attributes') = attr;
                    end
                end
            else
                ptextflag = 'Text';
                if (strcmp(name, '#cdata_dash_section'))
                    ptextflag = 'CDATA';
                elseif (strcmp(name, '#comment'))
                    ptextflag = 'Comment';
                end

                % this is the text in an element (i.e., the parentNode)
                if (~isempty(regexprep(text.(textflag),'[\s]*','')))
                    if (~isfield(ptext,ptextflag) || isempty(ptext.(ptextflag)))
                        ptext.(ptextflag) = text.(textflag);
                    else
                        % what to do when element data is as follows:
                        % <element>Text <!--Comment--> More text</element>

                        % put the text in different cells:
                        %  if (~iscell(ptext)) ptext = {ptext}; end
                        %  ptext{length(ptext)+1} = text;

                        % just append the text
                        ptext.(ptextflag) = [ptext.(ptextflag) text.(textflag)];
                    end
                end
            end

        end
    end
end

%  ----- Subfunction getNodeData -----
function [text,name,attr,childs,textflag] = getNodeData(theNode)
    %  Create structure of node info.

    % make sure name is allowed as structure name
    name = toCharArray(getNodeName(theNode))';
    name = strrep(name, '-', '_dash_');
    name = strrep(name, ':', '_colon_');
    name = strrep(name, '.', '_dot_');

    attr = parseAttributes(theNode);
    if (isempty(fieldnames(attr)))
        attr = [];
    end

    % parse child nodes
    [childs,text,textflag] = parseChildNodes(theNode);

    if (isempty(fieldnames(childs)) && isempty(fieldnames(text)))
        % get the data of any childless nodes
        %  faster than if any(strcmp(methods(theNode), 'getData'))
        %  no need to try-catch (?)
        %  faster than text = char(getData(theNode));
        text.(textflag) = toCharArray(getTextContent(theNode))';
    end

end

%  ----- Subfunction parseAttributes -----
function attributes = parseAttributes(theNode)
    %  Create attributes structure.

    attributes = struct;
    if hasAttributes(theNode)
       theAttributes = getAttributes(theNode);
       numAttributes = getLength(theAttributes);

       for count = 1:numAttributes
            % attrib = item(theAttributes,count-1);
            % attr_name = regexprep(char(getName(attrib)),'[-:.]','_');
            % attributes.(attr_name) = char(getValue(attrib));

            % Suggestion of Adrian Wanner
            str = toCharArray(toString(item(theAttributes,count-1)))';
            k = strfind(str,'=');
            attr_name = str(1:(k(1)-1));
            attr_name = strrep(attr_name, '-', '_dash_');
            attr_name = strrep(attr_name, ':', '_colon_');
            attr_name = strrep(attr_name, '.', '_dot_');
            attributes.(attr_name) = str((k(1)+2):(end-1));
       end
    end
end

author: opencobra

repo: cobratoolbox


text:
%% NLOPT Install for OPTI Toolbox
% Copyright (C) 2012 Jonathan Currie (I2C2)

% This file will help you compile NonLinear OPTimization (NLOPT) for use
% with MATLAB.

% My build platform:
% - Windows 8 SP1 x64
% - Visual Studio 2012

% To recompile you will need to get / do the following:

% 1) Get NLOPT
% NLOPT is available from http://ab-initio.mit.edu/wiki/index.php/NLopt.
% Download the source.

% 2) Compile NLOPT
% The easiest way to compile NLOPT is to use the Visual Studio Project
% Builder included with OPTI. Use the following commands, substituting the
% required path on your computer:
%
% %% Visual Studio Builder Commands
% path = 'full path to NLOPT here'; %e.g. 'C:\Solvers\NLOPT'
% sdir = path;
% name = 'libnlopt';
% opts.exPP = {'_CRT_SECURE_NO_WARNINGS'};
% opts.exclude = {'testfuncs.c','tst.cc','tstc.c','testros.cc','prog.cc','redblack_test.c','DIRparallel.c'};
% opts.exFolder = {'octave','swig','test'};
% VS_WriteProj(sdir,name,[],opts)
% %%
%
% Once complete, you will have a directory called NLOPT\libnlopt. Open the
% Visual Studio 2012 project file, then complete the following steps:
%   a) A default config.h has not been supplied (and I'm guessing only gets
%   made in Linux versions) so I've made one up based on the template. Copy
%   the supplied config.h from:
%       OPTI/Solvers/nlopt/Source/Include
%   To the project directory.
%   b) Build a Win32 or x64 Release to compile the code.
%   c) Copy the generated .lib file to the following folder:
%
%   OPTI/Solvers/nlopt/Source/lib/win32 or win64
%
%   You will also need to copy nlopt.h to the following folder:
%
%   OPTI/Solvers/nlopt/Source/Include
%
%   Note with Global Optimization turned on in VS it will take a long time
%   to link the MEX file.

% 3) NLOPT MEX Interface
% The NLOPT MEX Interface was written by Steven Johnson and is located in
% the octave folder (nlopt_optimize-mex.c) HOWEVER in its original form it
% was not compatible with OPTI (due to the method of adding nonlinear
% constraints cell by cell), as well as the file name caused compile
% problems (no '-' allowed for MEX). Therefore I have modified the
% MEX interface and included is nlopt_optimize_mex.c which you will need to
% use for compatibility with OPTI. You will need to however copy
% nlopt_optimize_usage.h from nlopt-xxx/octave to:
%
%   OPTI/Solvers/nlopt/Source/Include

% 4) Compile the MEX File
% The code below will automatically include all required libraries and
% directories to build the NLOPT MEX file. Once you have completed all
% the above steps, simply run this file to compile NLOPT! You MUST BE in
% the base directory of OPTI!

clear nlopt

% Get Arch Dependent Library Path
libdir = opti_GetLibPath();

fprintf('\n------------------------------------------------\n');
fprintf('NLOPT MEX FILE INSTALL\n\n');

%Get NLOPT Libraries
post = [' -IInclude -L' libdir ' -llibnlopt -llibut -output nlopt'];

%CD to Source Directory
cdir = cd;
cd 'Solvers/nlopt/Source';

%Compile & Move
pre = 'mex -v -largeArrayDims nloptmex.c';
try
    eval([pre post])
    movefile(['nlopt.' mexext],'../','f')
    fprintf('Done!\n');
catch ME
    cd(cdir);
    error('opti:nlopt','Error Compiling NLOPT!\n%s',ME.message);
end
cd(cdir);
fprintf('------------------------------------------------\n');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
if(exist('D:/Datasets/Bosphorus/', 'file'))
    Bosphorus_dir = 'D:\Datasets\Bosphorus/';
elseif(exist('E:/Datasets/Bosphorus/', 'file'))
    Bosphorus_dir = 'E:\Datasets\Bosphorus/';
elseif(exist('/media/tadas/2EBEA130BEA0F20F/datasets/Bosphorus/', 'file'))
    Bosphorus_dir = '/media/tadas/2EBEA130BEA0F20F/datasets/Bosphorus/';
else
    fprintf('Bosphorus dataset location not found (or not defined)\n');
end

hog_data_dir = ['D:\Datasets\face_datasets'];

all_recs = dir([Bosphorus_dir, '/BosphorusDB/BosphorusDB/bs*']);
all_recs_mat = cat(1, all_recs.name);
all_recs = cell(numel(all_recs), 1);

for i=1:size(all_recs_mat,1)

    all_recs{i} = all_recs_mat(i,:);

end

devel_recs = all_recs(1:3:end);
train_recs = setdiff(all_recs, devel_recs);

author: TadasBaltrusaitis

repo: OpenFace


text:
clear all; close all; clc;
addpath('tools');

db_name = '256_ObjectCategories'; % '256_ObjectCategories' as a option
numRetrieval = 10;

% load dataset
if strcmp(db_name, 'faceDataset')
    load feat4096Norml.mat;
    path_imgDB = './faceDataset/';
    addpath(path_imgDB);
    queryID = 1900;
elseif strcmp(db_name, '256_ObjectCategories')
    %load 256feat2048Norml.mat; % 0.495471
    % load 256feat4096Norml.mat; % 0.484413
    load feat4096Norml.mat;
    path_imgDB = './database/';
    addpath(path_imgDB);
    %example
    queryID = 1; %
end

%if not normalize, then do
% featNorm = normalize1(feat);
% save('feat4096Norml.mat','featNorm', 'rgbImgList');

% [pc, ~] = eigs(double(cov(feat)), 128);
% feat = feat*pc;

%virsulazation
retrieval_virsulazation( queryID, numRetrieval, feat_norm, imgNamList);

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function [Z, Z_L, Z_U, T, P, rho, c, g, mu, nu, k, n, n_sum] = atmo(alt,division,units)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Program:    1976 Standard Atmosphere Calculator[0-1000 km]
%   Author:     Brent Lewis(RocketLion@gmail.com)
%               University of Colorado-Boulder
%   History:    Original-1/10/2007
%               Revision-1/12/2007-Corrected for changes in Matlab versions
%               for backward compatability-Many thanks to Rich
%               Rieber(rrieber@gmail.com)
%   Input:      alt:        Final Geometric Altitude[km]
%               division:   Reporting points for output arrays[km]
%                           (.01 km & Divisible by .01 km)
%               units:      1-[Metric]
%                           2-{English}
%   Default:    Values used if no input
%               alt:        1000 km
%               division:   1 km
%               units:      Metric
%   Output:     Each value has a specific region that it is valid in with this model
%               and is only printed out in that region
%               Z:          Total Reporting Altitudes[0<=alt<=1000 km][km]{ft}
%               Z_L:        Lower Atmosphere Reporting Altitudes[0<=alt<=86 km][km]{ft}
%               Z_U:        Upper Atmosphere Reporting Altitudes[86<=alt<=1000 km][km]{ft}
%               T:          Temperature array[0<=alt<=1000 km][K]{R}
%               P:          Pressure array[0<=alt<=1000 km][Pa]{in_Hg}
%               rho:        Density array[0<=alt<=1000 km][kg/m^3]{lb/ft^3}
%               c:          Speed of sound array[0<=alt<=86 km][m/s]{ft/s}
%               g:          Gravity array[0<=alt<=1000 km][m/s^2]{ft/s^2}
%               mu:         Dynamic Viscosity array[0<=alt<=86 km][N*s/m^2]{lb/(ft*s)}
%               nu:         Kinematic Viscosity array[0<=alt<=86 km][m^2/s]{ft^2/s}
%               k:          Coefficient of Thermal Conductivity
%                           array[0<=alt<=86 km][W/(m*K)]{BTU/(ft*s*R)}
%               n:          Number Density of individual gases
%                           (N2 O O2 Ar He H)[86km<=alt<=1000km][1/m^3]{1/ft^3}
%               n_sum:      Number Density of total gases
%                           [86km<=alt<=1000km][1/m^3]{1/ft^3}
%   Acknowledgements:       1976 U.S. Standard Atmosphere
%                           Prof. Adam Norris-Numerical Analysis Class
%                           Steven S. Pietrobon USSA1976 Program
%   Notes:                  Program uses a 5-point Simpson's Rule in 10
%                           meter increments.  Results DO vary by less 1%
%                           compared to tabulated values and is probably
%                           caused by different integration techniques
%   Examples:               atmo() will compute the full atmosphere in 1 km
%                           increments and output in Metric Units
%                           atmo(10) will compute the atmosphere between 0
%                           and 10 km in 1 km increments and output in
%                           Metric Units
%                           atmo(20,.1,2) will compute the atmosphere
%                           between 0 and 20 km in 100 m increments and
%                           output in English Units
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin == 0
    alt = 1000;
    division = 1;
    units = 1;
elseif nargin == 1
    division = 1;
    units = 1;
elseif nargin == 2
    units = 1;
end

%   Error Reporting
if nargin > 3
    error('Too many inputs')
elseif mod(division,.01) ~= 0
    error('Divisions must be multiples of .01 km')
elseif units ~= 1 && units ~= 2
    error('Units Choice Invalid[1-Metric,2-English]')
elseif alt<0 || alt>1000
    error('Program only valid for 0<altitudes<1000 km')
end

%   Matrix Pre-allocation
if alt <= 86
    Z_L = (0:division:alt)';
    Z_U = [];
    n = [];
else
    Z_L = (0:division:86)';
    Z_U = (86:division:alt)';
    if mod(86,division) ~= 0
        Z_L = [Z_L; 86];
    end
    if mod(alt-86,division) ~= 0
        Z_U = [Z_U; alt];
    end
end
T_L = zeros(size(Z_L));
T_M_L = T_L;
T_U = zeros(size(Z_U));

%   Conversion Factor Used in 80<alt<86 km
Z_M = 80:.5:86;
M_M_0 = [1 .999996 .999989 .999971 .999941 .999909 ...
    .999870 .999829 .999786 .999741 .999694 .999641 .999579];

%   Constants
M_0 = 28.9644;
M_i = [28.0134; 15.9994; 31.9988; 39.948; 4.0026; 1.00797];
beta = 1.458e-6;
gamma = 1.4;
g_0 = 9.80665;
R = 8.31432e3;
r_E = 6.356766e3;
S = 110.4;
N_A = 6.022169e26;

%   Temperature
for i = 1 : length(Z_L)
    T_L(i,1) = atmo_temp(Z_L(i));
    T_M_L(i,1) = T_L(i,1);
    if Z_L(i) > 80 && Z_L(i) < 86
        T_L(i,1) = T_L(i)*interp1(Z_M,M_M_0,Z_L(i));
    end
end
for i = 1 : length(Z_U)
    T_U(i,1) = atmo_temp(Z_U(i));
end

%   Number Density
if alt > 86
    n = atmo_compo(alt,division);
    n_sum = sum(n,2);
else
    n = [];
    n_sum = [];
end

%   Pressure
P_L = atmo_p(Z_L);
P_U = atmo_p(Z_U,T_U,n_sum);

%   Density
rho_L = M_0*P_L./(R*T_M_L);
if ~isempty(P_U)
    rho_U = n*M_i/N_A;
else
    rho_U = [];
end

%   Speed of Sound
c = sqrt(gamma*R*T_M_L/M_0);
%   Dynamic Viscosity
mu = beta*T_L.^1.5./(T_L+S);
%   Kinematic Viscosity
nu = mu./rho_L;
%   Thermal Conductivity Coefficient
k = 2.64638e-3*T_L.^1.5./(T_L+245*10.^(-12./T_L));

%   Combine Models
T = [T_L(1:end-1*double(~isempty(T_U)));T_U];
P = [P_L(1:end-1*double(~isempty(T_U)));P_U];
rho = [rho_L(1:end-1*double(~isempty(T_U)));rho_U];
Z = [Z_L(1:end-1*double(~isempty(T_U)));Z_U];

%   Gravity
g = g_0*(r_E./(r_E+Z)).^2;

if units == 2
    unit_c = [3.048e-1 3.048e-1 3.048e-1 5/9 0.0001450377 1.6018463e1...
        3.048e-1 3.048e-1 1.488163944 9.290304e-2 6.226477504e-3...
        3.531466672e2 3.531466672e2];
    Z = Z/unit_c(1);
    Z_L = Z_L/unit_c(2);
    Z_U = Z_U/unit_c(3);
    T = T/unit_c(4);
    P = P/unit_c(5);
    rho = rho/unit_c(6);
    c = c/unit_c(7);
    g = g/unit_c(8);
    mu = mu/unit_c(9);
    nu = nu/unit_c(10);
    k = n/unit_c(11);
    n_sum = n_sum/unit_c(12);
end


text:
function varargout = ensureManifoldMesh(varargin)
%ENSUREMANIFOLDMESH Apply several simplification to obtain a manifold mesh.
%
%   Try to transform an input mesh into a manifold mesh.
%
%   Not all cases of "non-manifoldity" are checked, so please use with
%   care.
%
%   [V2, F2] = ensureManifoldMesh(V, F);
%   [V2, F2] = ensureManifoldMesh(MESH);
%   MESH2 = ensureManifoldMesh(...);
%
%   Example
%   ensureManifoldMesh
%
%   See also
%    meshes3d, isManifoldMesh

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2019-02-01, using Matlab 9.5.0.944444 (R2018b)
% Copyright 2019-2022 INRA - Cepia Software Platform

%% Parse input arguments

[vertices, faces] = parseMeshData(varargin{:});
verbose = true;

%% Pre-processing

% remove duplicate faces if any
if verbose
    disp('remove duplicate faces');
end
faces = removeDuplicateFaces(faces);

%% Iterative processing of multiple edges
% Reduces all edges connected to more than two faces, by collapsing second
% vertex onto the first one.

% iter = 0;
% while ~isManifoldMesh(vertices, faces) && iter < 10
%     iter = iter + 1;
    if verbose
        disp('collapse edges with many faces');
    end

    [vertices, faces] = collapseEdgesWithManyFaces(vertices, faces);
% end

%% Format output

varargout = formatMeshOutput(nargout, vertices, faces);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function m = similarity(sp1, sp2, opts)
% Returns the similarity score of two superpixels using various methods
% As seen with Matlab profiler, this function is responsible for majority
% of calculation for the whole algorithm.

d = zeros(1, length(opts.features));
a = 1;

for fn = 1:length(opts.features)

    % Intersection distance
    %d(fn) = 1 - sum(min([sp1.hist{fn}, sp2.hist{fn}], [], 2));
    % Hellinger's distance (Hellinger's kernel)
    %d(fn) = sqrt(sum((sqrt(sp1.hist{fn}) - sqrt(sp2.hist{fn})).^2));
    % l^1 distance (Intersection kernel)
    %d(fn) = sum(abs(sp1.hist{fn} - sp2.hist{fn}));

    if opts.features(fn) == 6 % size
        d(fn) = sp1.size + sp2.size;
        a = min(sp1.size, sp2.size); % for removing very small superpixels early
    else
        d(fn) = 0.5*sum(((sp1.hist{fn} - sp2.hist{fn}).^2)./(eps + sp1.hist{fn} + sp2.hist{fn})); % \chi^2 distance (\chi^2 kernel)
    end

end

d = 1 - d; % Low distance = high similarity. This is the monotone decreasing mapping we use
m = sum(d.*opts.feature_weights)/sum(opts.feature_weights); % The similarity value. If weights are ones, then this is just mean(d)

% Force small superpixels to merge first
if a < 0.002
    m = 10 + m; % adding constant will preserve the original ordering among overloaded pairs
end

author: Cloud-CV

repo: object-proposals


text:
function [uOutput] = import_Iceland_PT(nFunction, sFilename)

% Filter function switchyard
%%%%     CHANGE THESE LINES %%%%%%%%%%%
if nFunction == FilterOp.getDescription
    uOutput = 'Iceland (with P/T-axes)';
elseif nFunction == FilterOp.getWebpage
  uOutput = 'import_Iceland_PT_doc.html';
 %%%%    DO NOT CHANGE %%%%%%%%%%%

elseif nFunction == FilterOp.importCatalog
    % Read formated data
    mData = textread(sFilename, '%s', 'delimiter', '\n', 'whitespace', '');
    % Create empty catalog
    uOutput = zeros(length(mData), 12);
    % Loop through all lines of catalog and convert them
    mData = char(mData);
    l = find( mData == ' ' );
    mData(l) = '0';
    errc = 1;

        for i = 1:length(mData(:,1))
            if rem(i,100) == 0 ; disp([ num2str(i) ' of ' num2str(length(mData)) ' events processed ']); end
            try
                uOutput(i,:) = readvalues(mData,i,1);
            catch
                errc = errc + 1;
                if errc == 100
                    if stoploop
                        return
                    end
                end
                msg.dbfprintf('Import: Problem in line %d of %s. Line ignored.\n',i, sFilename);
                uOutput(i,:)=nan;
            end
        end

    l = isnan(uOutput(:,1));
    uOutput(l,:) = [];
end

%%%%%%%%%%%%%%%%%%%%%%
%%%%     CHANGE THESE LINES %%%%%%%%%%%

function [uOutput] = readvalues(mData,i,k)

uOutput(k,1) =-str2num(mData(i,29:30)) - (str2num(mData(i,31:35))/60); % Longitude
uOutput(k,2) = str2num(mData(i,20:21)) + (str2num(mData(i,22:26))/60); % Latitude
uOutput(k,3) = str2num(mData(i,2:3)) + 1900; % Year
uOutput(k,4) = str2num(mData(i,4:5));        % Month
uOutput(k,5) = str2num(mData(i,6:7));        % Day
uOutput(k,6) = str2num(mData(i,64:66));      % Magnitude
uOutput(k,7) = str2num(mData(i,37:41));      % Depth
uOutput(k,8) = str2num(mData(i,9:10));       % Hour
uOutput(k,9) = str2num(mData(i,11:12));      % Minute

fPTrend = str2double(mData(i,68:70));
fPDip = str2double(mData(i,72:73));
fTTrend = str2double(mData(i,75:77));
fTDip = str2double(mData(i,79:80));

% DipDirection, Dip, Rake
[uOutput(k,10), uOutput(k,11), uOutput(k,12)] = ex_pt2fm(fPTrend, fPDip, fTTrend, fTDip);

% Convert to decimal years
uOutput(k,3) = decyear([uOutput(k,3) uOutput(k,4) uOutput(k,5) uOutput(k,8) uOutput(k,9)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    DO NOT CHANGE %%%%%%%%%%%

function [mystop] = stoploop()

ButtonName=questdlg('More than 100 lines could not be read. Continue?', ...
    'Interrupt?', ...
    'Yes','No','Nope');

switch ButtonName
case 'Yes'
    disp('going on');
    mystop = 0;
case 'No'
    mystop = 1;
end % switch

author: CelsoReyes


text:
function [ball] = pwrbal1(pp,pw,ee)
%PWRBAL1 compute the ortho-balls associated with a 1-simplex
%triangulation embedded in R^2 or R^3.
%   [BB] = PWRBAL1(PP,PW,TT) returns the set of power balls
%   associated with the edges in [PP,PW,EE], such that BB =
%   [XC,YC,RC.^2]. PW is a vector of vertex weights.

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 02/05/2018

%---------------------------------------------- basic checks
    if ( ~isnumeric(pp) || ...
         ~isnumeric(pw) || ...
         ~isnumeric(ee) )
        error('pwrbal1:incorrectInputClass' , ...
            'Incorrect input class.');
    end

%---------------------------------------------- basic checks
    if (ndims(pp) ~= +2 || ...
        ndims(pw) ~= +2 || ...
        ndims(ee) ~= +2 )
        error('pwrbal1:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    if (size(pp,2) < +2 || ...
            size(pp,1)~= size(pw,1) || ...
                size(ee,2) < +2 )
        error('pwrbal1:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    switch (size(pp,2))
        case  +2
    %-------------------------------------------- lin offset
        pp12 = pp(ee(:,1),:) - ...
               pp(ee(:,2),:) ;

        ww12 = pw(ee(:,1),1) - ...
               pw(ee(:,2),1) ;

        dp12 = sum(pp12.*pp12,2) ;

        tpwr = +.5 * (ww12+dp12)./dp12 ;

        ball = zeros(size(ee,1),3) ;
        ball(:,1:2) = ...
            pp(ee(:,1),:) - tpwr.*pp12 ;

        vsq1 = ...
            pp(ee(:,1),:) - ball(:,1:2);
        vsq2 = ...
            pp(ee(:,2),:) - ball(:,1:2);

    %-------------------------------------------- mean radii
        rsq1 = sum(vsq1 .^ 2,2) ;
        rsq2 = sum(vsq2 .^ 2,2) ;

        rsq1 = rsq1-pw(ee(:,1)) ;
        rsq2 = rsq2-pw(ee(:,2)) ;

        ball(:,3) = (rsq1 + rsq2) / 2. ;

        case  +3
    %-------------------------------------------- lin offset
        pp12 = pp(ee(:,1),:) - ...
               pp(ee(:,2),:) ;

        ww12 = pw(ee(:,1),1) - ...
               pw(ee(:,2),1) ;

        dp12 = sum(pp12.*pp12,2) ;

        tpwr = +.5 * (ww12+dp12)./dp12 ;

        ball = zeros(size(ee,1),4) ;
        ball(:,1:3) = ...
            pp(ee(:,1),:) - tpwr.*pp12 ;

        vsq1 = ...
            pp(ee(:,1),:) - ball(:,1:3);
        vsq2 = ...
            pp(ee(:,2),:) - ball(:,1:3);

    %-------------------------------------------- mean radii
        rsq1 = sum(vsq1 .^ 2,2) ;
        rsq2 = sum(vsq2 .^ 2,2) ;

        rsq1 = rsq1-pw(ee(:,1)) ;
        rsq2 = rsq2-pw(ee(:,2)) ;

        ball(:,4) = (rsq1 + rsq2) / 2. ;

    otherwise

    error('pwrbal2:unsupportedDimension' , ...
            'Dimension not supported.');

    end

end

author: dengwirda


text:
% Minimum volume ellipsoid covering a finite set
% Section 8.4.1, Boyd & Vandenberghe "Convex Optimization"
% Almir Mutapcic - 10/05
% (a figure is generated)
%
% Given a finite set of points x_i in R^2, we find the minimum volume
% ellipsoid (described by matrix A and vector b) that covers all of
% the points by solving the optimization problem:
%
%           maximize     log det A
%           subject to   || A x_i + b || <= 1   for all i
%
% CVX cannot yet handle the logdet function, but this problem can be
% represented in an equivalent way as follows:
%
%           maximize     det(A)^(1/n)
%           subject to   || A x_i + b || <= 1   for all i
%
% The expression det(A)^(1/n) is SDP-representable, and is implemented
% by the MATLAB function det_rootn().

% Generate data
x = [ 0.55  0.0;
      0.25  0.35
     -0.2   0.2
     -0.25 -0.1
     -0.0  -0.3
      0.4  -0.2 ]';
[n,m] = size(x);

% Create and solve the model
cvx_begin
    variable A(n,n) symmetric
    variable b(n)
    maximize( det_rootn( A ) )
    subject to
        norms( A * x + b * ones( 1, m ), 2 ) <= 1;
cvx_end

% Plot the results
clf
noangles = 200;
angles   = linspace( 0, 2 * pi, noangles );
ellipse  = A \ [ cos(angles) - b(1) ; sin(angles) - b(2) ];
plot( x(1,:), x(2,:), 'ro', ellipse(1,:), ellipse(2,:), 'b-' );
axis off

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function jac = p21_jac ( neqn, t, y )

%*****************************************************************************80
%
%% P21_JAC evaluates the jacobian for problem p21.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the jacobian.
%
%    Output, real JAC(NEQN,NEQN), the jacobian matrix.
%
  jac = zeros ( neqn, neqn );

  jac(1,1) = 0.0;
  jac(1,2) = 1.0;

  jac(2,1) = - 1.0 + 0.25 / ( t + 1.0 ).^2;
  jac(2,2) = - 1.0 / ( t + 1.0 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Gamma,Gammasum,Xi,LL] = hsinference(data,T,hmm,residuals,options,XX,Gamma0,Xi0,cv)
%
% inference engine for HMMs.
%
% INPUT
%
% data      Observations - a struct with X (time series) and C (classes)
% T         Number of time points for each time series
% hmm       hmm data structure
% residuals in case we train on residuals, the value of those.
% XX        optionally, XX, as computed by setxx.m, can be supplied
% Gamma0    initial Gamma (only used for ehmm model)
% cv        is this called from a CV routine?
%
% OUTPUT
%
% Gamma     Probability of hidden state given the data
% Gammasum  sum of Gamma over t
% Xi        joint Prob. of child and parent states given the data
% LL        Log-likelihood, summed across time points
%
% Author: Diego Vidaurre, OHBA, University of Oxford

if length(hmm.train.embeddedlags_batched) > 1
    [Gamma,Gammasum,Xi,LL] = hsinference_batched(T,hmm,residuals,XX);
    return
end

N = length(T);
K = hmm.K;
p = hmm.train.lowrank; do_HMM_pca = (p > 0);
if nargin < 7, Gamma0 = []; end
if nargin < 8, Xi0 = []; end
if nargin < 9, cv = 0; end
if ~isfield(hmm,'train')
    if nargin<5 || isempty(options)
        error('You must specify the field options if hmm.train is missing');
    end
    hmm.train = checkoptions(options,data.X,T,0);
end
mixture_model = isfield(hmm.train,'id_mixture') && hmm.train.id_mixture;
if isfield(hmm.train,'episodic'), episodic = hmm.train.episodic;
else, episodic = 0;
end
if episodic, rangeK = 1:K+1; else, rangeK = 1:K; end
order = hmm.train.maxorder;

if iscell(data)
    data = cell2mat(data);
end
if ~isstruct(data)
    data = struct('X',data);
    data.C = NaN(size(data.X,1)-order*length(T),K);
end

if do_HMM_pca
    ndim = size(data.X,2);
elseif nargin<4 || isempty(residuals)
    ndim = size(data.X,2);
    if ~isfield(hmm.train,'Sind')
        orders = formorders(hmm.train.order,hmm.train.orderoffset,hmm.train.timelag,hmm.train.exptimelag);
        hmm.train.Sind = formindexes(orders,hmm.train.S) == 1;
        if ~hmm.train.zeromean, hmm.train.Sind = [true(1,ndim); hmm.train.Sind]; end
    end
    residuals =  getresiduals(data.X,T,hmm.train.S,hmm.train.maxorder,hmm.train.order,...
        hmm.train.orderoffset,hmm.train.timelag,hmm.train.exptimelag,hmm.train.zeromean);
else
    ndim = size(residuals,2);
end

if (episodic && ~isfield(hmm.state(1),'P')) || (~episodic && ~isfield(hmm,'P'))
    hmm = hmmhsinit(hmm);
end

if nargin<6 || isempty(XX)
    setxx;
end

if ~hmm.train.acrosstrial_constrained
    Gamma = cell(N,1);
    LL = cell(N,1);
end
Gammasum = zeros(N,K);
if ~mixture_model
    Xi = cell(N,1);
else
    Xi = [];
end

setstateoptions;
% Cache shared results for use in obslike
for k = rangeK
    %hmm.cache = struct();
    hmm.cache.train{k} = train;
    %hmm.cache.order{k} = order;
    %hmm.cache.orders{k} = orders;
    %hmm.cache.Sind{k} = Sind;
    %hmm.cache.S{k} = S;
    if do_HMM_pca
        W = hmm.state(k).W.Mu_W;
        v = hmm.Omega.Gam_rate / hmm.Omega.Gam_shape;
        C = W * W' + v * eye(ndim);
        ldetWishB = 0.5*logdet(C); PsiWish_alphasum = 0;
    elseif k == 1 && strcmp(train.covtype,'uniquediag') || strcmp(train.covtype,'shareddiag')
        ldetWishB = 0;
        PsiWish_alphasum = 0;
        for n = 1:ndim
            if ~regressed(n), continue; end
            ldetWishB = ldetWishB+0.5*log(hmm.Omega.Gam_rate(n));
            PsiWish_alphasum = PsiWish_alphasum+0.5*psi(hmm.Omega.Gam_shape);
        end
        C = hmm.Omega.Gam_shape ./ hmm.Omega.Gam_rate;
    elseif k == 1 && strcmp(train.covtype,'uniquefull') || strcmp(train.covtype,'sharedfull')
        ldetWishB = 0.5*logdet(hmm.Omega.Gam_rate(regressed,regressed));
        PsiWish_alphasum = 0;
        for n = 1:sum(regressed)
            PsiWish_alphasum = PsiWish_alphasum+psi(hmm.Omega.Gam_shape/2+0.5-n/2);
        end
        PsiWish_alphasum=PsiWish_alphasum*0.5;
        C = hmm.Omega.Gam_shape * hmm.Omega.Gam_irate;
        for kk = rangeK
            [hmm.cache.WishTrace{kk},hmm.cache.codevals] = computeWishTrace(hmm,regressed,XX,C,kk);
        end
    elseif strcmp(train.covtype,'diag')
        ldetWishB=0;
        PsiWish_alphasum = 0;
        for n = 1:ndim
            if ~regressed(n), continue; end
            ldetWishB = ldetWishB+0.5*log(hmm.state(k).Omega.Gam_rate(n));
            PsiWish_alphasum = PsiWish_alphasum+0.5*psi(hmm.state(k).Omega.Gam_shape);
        end
        C = hmm.state(k).Omega.Gam_shape ./ hmm.state(k).Omega.Gam_rate;
        if episodic, iC = hmm.state(k).Omega.Gam_rate / hmm.state(k).Omega.Gam_shape; end
    elseif strcmp(train.covtype,'full')
        ldetWishB = 0.5*logdet(hmm.state(k).Omega.Gam_rate(regressed,regressed));
        PsiWish_alphasum = 0;
        for n = 1:sum(regressed)
            PsiWish_alphasum = PsiWish_alphasum+0.5*psi(hmm.state(k).Omega.Gam_shape/2+0.5-n/2);
        end
        C = hmm.state(k).Omega.Gam_shape * hmm.state(k).Omega.Gam_irate;
        if episodic, iC = hmm.state(k).Omega.Gam_rate / hmm.state(k).Omega.Gam_shape; end
        [hmm.cache.WishTrace{k},hmm.cache.codevals] = computeWishTrace(hmm,regressed,XX,C,k);
    end
    % Set up cache
    if ~isfield(train,'distribution') || strcmp(train.distribution,'Gaussian')
        hmm.cache.ldetWishB{k} = ldetWishB;
        hmm.cache.PsiWish_alphasum{k} = PsiWish_alphasum;
        hmm.cache.C{k} = C;
        if episodic && (strcmp(train.covtype,'diag') || strcmp(train.covtype,'full'))
            hmm.cache.iC{k} = iC;
        end
    end
end

if hmm.train.acrosstrial_constrained % state probabilities are shared across trials

    ttrial = T(1)-order; L = zeros(N*ttrial,K);
    if episodic
        [Gammastar,Xistar,LL] = fb_Gamma_inference_ehmm(XX,hmm,residuals,Gamma0,Xi0,T);
    else
        for j = 1:N
            t = (1:ttrial) + (j-1)*ttrial;
            t2 = (1:ttrial+order) + (j-1)*(ttrial+order);
            try
                L(t2,:) = obslike([],hmm,residuals(t,:),XX(t,:),hmm.cache,cv);
            catch
                error('obslike function is giving trouble - Out of precision?')
            end
        end
        L = squeeze(sum(reshape(L,[ttrial+order,N,K]),2));
        if mixture_model
            Gammastar = id_Gamma_inference(L,hmm.Pi,order); Xistar = [];
        else
            [Gammastar,Xistar,LL] = fb_Gamma_inference_sub(L,hmm.P,hmm.Pi,T(1),...
                order,hmm.train.Gamma_constraint);
        end
    end
    Gamma = zeros(ttrial,N,K);
    if episodic
        Xi = zeros(ttrial-1,N,K,4);
    elseif mixture_model
        Xi = zeros(ttrial-1,N,K,K);
    end
    for t = 1:ttrial
        Gamma(t,:,:) = repmat(Gammastar(t,:),N,1);
        if t==ttrial, break; end
        if episodic
            Xi(t,:,:,:) = reshape(repmat(Xistar(t,:,:),N,1),[1,N,K,4]);
        elseif mixture_model
            Xi(t,:,:,:) = reshape(repmat(Xistar(t,:,:),N,1),[1,N,K,K]);
        end
    end
    Gamma = reshape(Gamma,ttrial*N,K);
    if episodic
        Xi = reshape(Xi,(ttrial-1)*N,K,2,2);
    elseif mixture_model
        Xi = reshape(Xi,(ttrial-1)*N,K,K);
    end

elseif hmm.train.useParallel==1 && N>1

    % to duplicate this code is really ugly but there doesn't seem to be
    % any other way - more Matlab's fault than mine

    % to avoid the entire data being copied entirely in each parfor loop
    XX_copy = cell(N,1); residuals_copy = cell(N,1);
    Gamma0_copy = cell(N,1); Xi0_copy = cell(N,1);
    C_copy = cell(N,1);
    for j = 1:N
        t0 = sum(T(1:j-1)); s0 = t0 - order*(j-1); s02 = t0 - (order+1)*(j-1);
        XX_copy{j} = XX(s0+1:s0+T(j)-order,:);
        if ~do_HMM_pca
            residuals_copy{j} = residuals(s0+1:s0+T(j)-order,:);
        end
        if isfield(data,'C')
            C_copy{j} = data.C(s0+1:s0+T(j)-order,:);
        end
        if ~isempty(Gamma0)
            Gamma0_copy{j} = Gamma0(s0+1:s0+T(j)-order,:);
        end
        if ~isempty(Xi0)
            Xi0_copy{j} = Xi(s02+1:s02+T(j)-order-1,:,:,:);
        end
    end; clear data; clear Gamma0; clear Xi0; clear XX; clear residuals

    parfor j = 1:N
        xit = []; llt = [];
        if order>0
            R = [zeros(order,size(residuals_copy{j},2)); residuals_copy{j}];
            if ~isempty(C_copy{j})
                C = [zeros(order,K); C_copy{j}];
            else
                C = NaN(size(R,1),K);
            end
        else
            if do_HMM_pca
                R = XX_copy{j};
            else
                R = residuals_copy{j};
            end
            if ~isempty(C_copy{j})
                C = C_copy{j};
            else
                C = NaN(size(R,1),K);
            end
        end
        % we jump over the fixed parts of the chain
        t = order+1;
        nsegments = computeNoSegments(T(j),t,C);
        xi = cell(nsegments,1); gamma = cell(nsegments,1);
        ll = cell(nsegments,1); ns = 1;
        while t <= T(j)
            if isnan(C(t,1)), no_c = find(~isnan(C(t:T(j),1)));
            else, no_c = find(isnan(C(t:T(j),1)));
            end
            if t>order+1
                if isempty(no_c), slicer = (t-1):T(j); %slice = (t-order-1):T(in);
                else, slicer = (t-1):(no_c(1)+t-2); %slice = (t-order-1):(no_c(1)+t-2);
                end
            else
                if isempty(no_c), slicer = t:T(j); %slice = (t-order):T(in);
                else, slicer = t:(no_c(1)+t-2); %slice = (t-order):(no_c(1)+t-2);
                end
            end
            slicepoints = slicer - order;
            XXt = XX_copy{j}(slicepoints,:);
            if isnan(C(t,1))
                if episodic
                    slicer2 = slicer(1:end-1);
                    Gamma0t = Gamma0_copy{j}(slicer,:);
                    Xi0t = []; if ~isempty(Xi0_copy{j}), Xi0t = Xi0_copy{j}(slicer2,:); end
                    [gammat,xit,llt] = fb_Gamma_inference_ehmm(XXt,hmm,R(slicer,:),Gamma0t,Xi0t);
                else
                    [gammat,xit,llt] = ...
                        fb_Gamma_inference(XXt,hmm,R(slicer,:),slicepoints,hmm.train.Gamma_constraint,cv);
                end
            else
                gammat = zeros(length(slicer),K);
                if t==order+1, gammat(1,:) = C(slicer(1),:); end
                if ~mixture_model
                    if episodic, xit = zeros(length(slicer)-1, K, 4);
                    else, xit = zeros(length(slicer)-1, K^2);
                    end
                end
                for i = 2:length(slicer)
                    gammat(i,:) = C(slicer(i),:);
                    if ~mixture_model
                        if episodic
                            for k = 1:K
                                gg = [gammat(i-1,k) (1-gammat(i-1,k))];
                                xitr = gg' * gg;
                                xit(i-1,k,:) = xitr(:)';
                            end
                        else
                            xitr = gammat(i-1,:)' * gammat(i,:) ;
                            xit(i-1,:) = xitr(:)';
                        end
                    end
                end
            end
            if t>order+1
                gammat = gammat(2:end,:);
            end
            if ~mixture_model
                xi{ns} = xit;
            end
            gamma{ns} = gammat;
            ll{ns} = llt;
            ns = ns + 1;
            if isempty(no_c), break;
            else, t = no_c(1)+t-1;
            end
        end
        Gamma{j} = cell2mat(gamma);
        Gammasum(j,:) = sum(Gamma{j});
        LL{j} = cell2mat(ll);
        if ~mixture_model
            if episodic
                Xi{j} = reshape(cell2mat(xi),T(j)-order-1,K,2,2);
            else
                Xi{j} = reshape(cell2mat(xi),T(j)-order-1,K,K);
            end
        end
    end

else

    for j = 1:N % this is exactly the same than the code above but changing parfor by for

        t0 = sum(T(1:j-1)); s0 = t0 - order*(j-1); s02 = t0 - (order+1)*(j-1);
        if order > 0
            R = [zeros(order,size(residuals,2)); residuals(s0+1:s0+T(j)-order,:)];
            if isfield(data,'C')
                C = [zeros(order,K); data.C(s0+1:s0+T(j)-order,:)];
            else
                C = NaN(size(R,1),K);
            end
        else %s0==t0
            if do_HMM_pca
                R = XX(s0+1:s0+T(j),:);
            else
                R = residuals(s0+1:s0+T(j),:);
            end
            if isfield(data,'C')
                C = data.C(s0+1:s0+T(j),:);
            else
                C = NaN(size(R,1),K);
            end
        end
        % we jump over the fixed parts of the chain
        t = order+1;
        nsegments = computeNoSegments(T(j),t,C);
        xi = cell(nsegments,1); gamma = cell(nsegments,1);
        ll = cell(nsegments,1); ns = 1;
        while t <= T(j)
            if isnan(C(t,1)), no_c = find(~isnan(C(t:T(j),1)));
            else, no_c = find(isnan(C(t:T(j),1)));
            end
            if t > order+1
                if isempty(no_c), slicer = (t-1):T(j); %slice = (t-order-1):T(in);
                else, slicer = (t-1):(no_c(1)+t-2); %slice = (t-order-1):(no_c(1)+t-2);
                end
            else
                if isempty(no_c), slicer = t:T(j); %slice = (t-order):T(in);
                else, slicer = t:(no_c(1)+t-2); %slice = (t-order):(no_c(1)+t-2);
                end
            end
            slicepoints = slicer + s0 - order;
            XXt = XX(slicepoints,:);
            if isnan(C(t,1))
                if episodic
                    slicepoints2 = slicer(1:end-1) + s02 - order;
                    Gamma0t = Gamma0(slicepoints,:);
                    Xi0t = []; if ~isempty(Xi0), Xi0t = Xi0(slicepoints2,:,:,:); end
                    [gammat,xit,llt] = ...
                        fb_Gamma_inference_ehmm(XXt,hmm,R(slicer,:),Gamma0t,Xi0t);
                else
                    [gammat,xit,llt] = ...
                        fb_Gamma_inference(XXt,hmm,R(slicer,:),slicepoints,hmm.train.Gamma_constraint,cv);
                end
                if any(isnan(gammat(:))) % this will never come up - we treat it within fb_Gamma_inference
                    error(['State time course inference returned NaN (out of precision). ' ...
                        'There are probably extreme events in the data'])
                end
            else
                gammat = zeros(length(slicer),K);
                llt = NaN(length(slicer),1);
                if t==order+1, gammat(1,:) = C(slicer(1),:); end
                if ~mixture_model
                    if episodic, xit = zeros(length(slicer)-1, K, 4);
                    else, xit = zeros(length(slicer)-1, K^2);
                    end
                end
                for i = 2:length(slicer)
                    gammat(i,:) = C(slicer(i),:);
                    if ~mixture_model
                        if episodic
                            for k = 1:K
                                gg = [gammat(i-1,k) (1-gammat(i-1,k))];
                                xitr = gg' * gg;
                                xit(i-1,k,:) = xitr(:)';
                            end
                        else
                            xitr = gammat(i-1,:)' * gammat(i,:) ;
                            xit(i-1,:) = xitr(:)';
                        end
                    end
                end
            end
            if t > order+1
                gammat = gammat(2:end,:);
            end
            if ~mixture_model
                xi{ns} = xit;
            end
            gamma{ns} = gammat;
            ll{ns} = llt;
            ns = ns + 1;
            if isempty(no_c), break;
            else, t = no_c(1)+t-1;
            end
        end
        Gamma{j} = cell2mat(gamma);
        Gammasum(j,:) = sum(Gamma{j});
        LL{j} = cell2mat(ll);
        if ~mixture_model
            if episodic
                Xi{j} = reshape(cell2mat(xi),T(j)-order-1,K,2,2);
            else
                Xi{j} = reshape(cell2mat(xi),T(j)-order-1,K,K);
            end
        end
    end

end

% join
if ~hmm.train.acrosstrial_constrained
    Gamma = cell2mat(Gamma);
    LL = cell2mat(LL);
    if ~mixture_model, Xi = cell2mat(Xi); end
end

% orthogonalise = 1;
% Gamma0 = Gamma;
% if orthogonalise && hmm.train.tuda
%     T = length(Gamma) / length(T) * ones(length(T),1);
%     Gamma = reshape(Gamma,[T(1) length(T) size(Gamma,2) ]);
%     Y = reshape(residuals(:,end),[T(1) length(T)]);
%     for t = 1:T(1)
%        y = Y(t,:)';
%        for k = 1:size(Gamma,3)
%            x = Gamma(t,:,k)';
%            b = y \ x;
%            Gamma(t,:,k) = Gamma(t,:,k) - (y * b)';
%        end
%     end
%     Gamma = reshape(Gamma,[T(1)*length(T) size(Gamma,3) ]);
%     Gamma = rdiv(Gamma,sum(Gamma,2));
%     Gamma = Gamma - min(Gamma(:));
%     Gamma = rdiv(Gamma,sum(Gamma,2));
% end

end

function nsegments = computeNoSegments(Tj,t0,C)
t = t0; nsegments = 0;
while t <= Tj
    if isnan(C(t,1)), no_c = find(~isnan(C(t:Tj,1)));
    else, no_c = find(isnan(C(t:Tj,1)));
    end
    nsegments = nsegments + 1;
    if isempty(no_c), break;
    else, t = no_c(1)+t-1;
    end
end
end

author: OHBA-analysis

repo: HMM-MAR


text:
function res = format(varargin)
%FINISH       Called when exiting Matlab: set rounding mode to nearest
%

% written  05/09/09     S.M. Rump
%

  setround(0)

repo: OpenMAS


text:
%% Test
disp('Accuracy per subject:')
disp(acc)
disp('Average accuracy:')
disp(mean(acc))

confusionchart(total_conf);

author: PatternRecognition

repo: OpenBMI


text:
classdef(Abstract) AbstractThrottleCurve < matlab.mixin.SetGet & matlab.mixin.Copyable
    %AbstractThrottleCurve Summary of this class goes here
    %   Detailed explanation goes here

    properties
        elems AbstractCurveElement

        constValue(1,1) double = NaN
    end

    properties(Transient)
        curve
    end

    methods
        function addElement(obj, newElement)
            obj.elems(end+1) = newElement;
            obj.sortElems();
            obj.generateCurve();
        end

        function removeElement(obj, element)
            obj.elems(obj.elems == element) = [];
            obj.sortElems();
            obj.generateCurve();
        end

        function generateCurve(obj)
            obj.sortElems();

            if(length(obj.elems) > 2)
                x = [obj.elems.indepVar];
                y = [obj.elems.depVar];

%                 xc = [max(x)];
%                 yc = [0];
%                 cc = [0;1];
%                 con = struct('xc',xc,'cc',cc,'yc',yc);
%                 obj.curve = splinefit(x,y,x,2); %I set things to use a piecewise linear interpolation now. 2020/12/18
                obj.curve = griddedInterpolant(x,y,'linear','nearest');

                if(all(not(diff(y))))
                    obj.constValue = y(1);
                else
                    obj.constValue = NaN;
                end
            elseif(length(obj.elems) == 2)
                x = [obj.elems.indepVar];
                y = [obj.elems.depVar];

%                 obj.curve = splinefit(x,y,1,2);
                obj.curve = griddedInterpolant(x,y,'linear','nearest');

                if(y(1) == y(2))
                    obj.constValue = y(1);
                else
                    obj.constValue = NaN;
                end
            else
                error('Cannot generate throttle curve: the number of elements in the curve must be greater than or equal to 2.');
            end
        end

        function yq = evalCurve(obj, xq)
            if(isnan(obj.constValue))
%                 if(isempty(obj.curve))
%                     obj.generateCurve();
%                 end

                yq = obj.curve(xq);
            else
                yq = ones(size(xq)) * obj.constValue;
            end
        end

        function sortElems(obj)
            [~,I] = sort([obj.elems.indepVar], 'ascend');
            obj.elems = obj.elems(I);
        end

        function [x, y] = getPlotablePoints(obj)
            obj.sortElems();

            x = [obj.elems.indepVar];
            y = [obj.elems.depVar];
        end
    end

    methods(Static)
        function obj = loadobj(obj)
            obj.generateCurve();
        end
    end

    methods(Abstract)
        [listBoxStr, elemArr] = getListboxStr(obj)

        curveName = getCurveName(obj)

        indepVarName = getIndepVarName(obj)
        indepVarUnit = getIndepVarUnit(obj)
        depVarName = getDepVarName(obj)
        depVarUnit = getDepVarUnit(obj)

        newElem = createNewElement(obj)

        listBoxTooltipStr = getListboxTooltipStr(obj)
    end

    methods(Abstract, Access = protected)
        newObj = copyElement(obj)
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ fea, out ] = ex_navierstokes17( varargin )
%EX_NAVIERSTOKES17 2D Example for incompressible turbulent flow in a channel.
%
%   [ FEA, OUT ] = EX_NAVIERSTOKES17( VARARGIN ) Sets up and solves
%   stationary turbulent flow at Re = 42800 in a between two flat
%   parallel plates. The inflow profile is constant and the outflow
%   should assume a fully developed turbulent profile (Ref. Laufer,
%   J. 1950). Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       igrid       scalar {3}             Grid refinement level
%                                          (>0=quadrilaterals, <0=triangles)
%       sf_u        string {sflag1}        Shape function for velocity
%       sf_p        string {sflag1}        Shape function for pressure
%       iplot       scalar 0/{1}           Plot solution and error (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { ...
            'igrid',    3;
            'solver',   '';
            'sf_u',     'sflag1';
            'sf_p',     'sflag1';
            'iplot',    1;
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Model parameters.
Re        = 42800;
rho       = 1;           % Density.
miu       = 1/Re;        % Molecular/dynamic viscosity.
uin       = 1;
% Geometry and grid parameters.
h         = 1;           % Height of rectangular domain.
l         = 5;           % Length of rectangular domain.
% Discretization parameters.
sf_u      = opt.sf_u;    % FEM shape function type for velocity.
sf_p      = opt.sf_p;    % FEM shape function type for pressure.

% Geometry definition.
gobj = gobj_rectangle( 0, l, -h/2, h/2 );
fea.geom.objects = { gobj };
fea.sdim = { 'x', 'y' };   % Coordinate names.

% Grid generation.
fea.grid = rectgrid( linspace(0,5,20), [-0.5 -0.45 -0.35 -0.2 0 0.2 0.35 0.45 0.5] );
for i=1:abs(opt.igrid)-1
  fea.grid = gridrefine( fea.grid, fid );
end
if( opt.igrid<0 )
  fea.grid = quad2tri( fea.grid );
end

% Problem definition.
fea = addphys( fea, @navierstokes );
fea.phys.ns.eqn.coef{1,end} = { rho };
fea.phys.ns.eqn.coef{2,end} = { miu };
fea.phys.ns.eqn.coef{5,end} = { uin };
if( strcmp(opt.solver,'openfoam') )
  fea.phys.ns.sfun = { 'sflag1', 'sflag1', 'sflag1' };
else
  fea.phys.ns.sfun = { sf_u sf_u sf_p };           % Set shape functions.
end

% Boundary conditions.
i_inflow  = 4;
i_outflow = 2;
fea.phys.ns.bdr.sel(i_inflow)  = 2;
fea.phys.ns.bdr.sel(i_outflow) = 3;
fea.phys.ns.bdr.coef{2,end}{1,i_inflow} = uin;     % Set inflow velocity.

fea.phys.ns.prop.turb.model = 'algebraic';
if( strcmp(opt.solver,'openfoam') )
  fea.phys.ns.prop.turb.model = 'SpalartAllmaras';
end
fea = parsephys(fea);   % Check and parse physics modes.

% Parse and solve problem.
fea = parseprob(fea);             % Check and parse problem struct.
if( strcmp(opt.solver,'openfoam') )
  logfid = fid; if( ~got.fid ), fid = []; end
  fea.sol.u = openfoam( fea, 'fid', fid, 'logfid', logfid );
  fid = logfid;
else
  mxnit = 25;
  nlrlx = 0.5;
  told  = 1e-3;
  tolc  = 1e-4;
  fea.sol.u = solvestat( fea, 'fid', fid, 'maxnit', mxnit, 'nlrlx', nlrlx, 'tolchg', tolc, 'toldef', told );
end

% Postprocessing.
y_ref = [0,0.09830508474576277,0.2101694915254238,0.31186440677966104,0.3661016949152544,0.4508474576271188,0.5491525423728814,0.6271186440677967,0.7016949152542376,0.7864406779661017,0.8440677966101697,0.898305084745763,0.9593220338983053,0.9796610169491529,0.9864406779661019,1]/2;
u_ref = [0.9972834919897843,0.9886928256326912,0.9780357557464595,0.9674135128859997,0.9487346180636177,0.9278964476433715,0.9029022521476668,0.8718133271418624,0.8469003947062924,0.7972951009983751,0.7457278848386352,0.6962270722080338,0.5357441374506622,0.32608544230322756,0.19455537497097786,0]';
n = length(y_ref);
x = 0.96*l*ones(1,n);
u = evalexpr( 'u', [x;y_ref], fea );
if( opt.iplot>0 )
  figure
  subplot(2,1,1)
  postplot( fea, 'surfexpr', 'sqrt(u^2+v^2)', 'isoexpr', 'sqrt(u^2+v^2)' )
  title('Velocity field')
  subplot(2,1,2)
  postplot( fea, 'surfexpr', 'p' )
  title('Pressure')

  figure
  plot( u/max(u), y_ref, 'b.-' )
  hold on
  plot( u_ref, y_ref, 'ro--' )
  ylabel('y')
  xlabel('u/u_{max}')
  grid on
  legend('Computed','Laufer (1950)','location','southwest')
end

% Error checking.
err = sqrt(sum((u-u_ref).^2)/sum(u_ref.^2));

if( ~isempty(fid) )
  fprintf(fid,'\nL2 Error: %f\n',err)
  fprintf(fid,'\n\n')
end

out.err  = err;
out.pass = err<0.26;
if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function c8vec_unity_test ( )

%*****************************************************************************80
%
%% C8VEC_UNITY_TEST tests C8VEC_UNITY;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2015
%
%  Author:
%
%    John Burkardt
%
  n = 12;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8VEC_UNITY_TEST\n' );
  fprintf ( 1, '  C8VEC_UNITY sets A to the N roots of unity\n' );

  a = c8vec_unity ( n );

  c8vec_print ( n, a, '  The N roots of unity:' );

  return
end


text:
% A helper function to check that all HOWTOs are working

HOWTO_spectra;
HOWTO_features;
HOWTO_sinusoidal;
HOWTO_egg;
HOWTO_glottalsource;
HOWTO_VAD;
HOWTO_envelope;
HOWTO_formant;
HOWTO_hmpd;
HOWTO_filterbanks;

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [Lmk, Frm, Fac] = makeMeasFactor(Lmk, Obs, Frm, Fac)

% MAKEMEASFACTOR Make measurement factor.
%   [Lmk, Frm, Fac] = MAKEMEASFACTOR(Lmk, Obs, Frm, Fac) creates a new
%   measurement factor Fac linking frame Frm to landmark Lmk, using the
%   observation informaiton in Obs.

%   Copyright 2015-     Joan Sola @ IRI-UPC-CSIC

Fac.used = true; % Factor is being used ?
Fac.id = newId; % Factor unique ID

Fac.type = 'measurement'; % {'motion','measurement','absolute'}
Fac.rob = Frm.rob;
Fac.sen = Obs.sen;
Fac.lmk = Obs.lmk;
Fac.frames = Frm.frm;

% Ranges
Fac.state.r1 = Frm.state.r;
Fac.state.r2 = Lmk.state.r;

% Measurement is the straight data
Fac.meas.y = Obs.meas.y;
Fac.meas.R = Obs.meas.R;
Fac.meas.W = Obs.meas.R^-1; % measurement information matrix

% Expectation has zero covariance -- and info in not defined
Fac.exp.e = Fac.meas.y; % expectation
Fac.exp.E = zeros(size(Obs.meas.R)); % expectation cov
%     Fac.exp.W = Fac.meas.W; % expectation information matrix

% Error is zero at this stage, and takes covariance and info from measurement
Fac.err.z     = zeros(size(Fac.meas.y)); % error or innovation (we call it error because we are on graph SLAM)
Fac.err.Z     = Fac.meas.R;              % error cov matrix
Fac.err.W     = Fac.meas.W;              % error information matrix
Fac.err.Wsqrt = chol(Fac.err.W);

% Jacobians are zero at this stage. Just make size correct.
Fac.err.J1 = Obs.Jac.E_r; % Jac. of error wrt. node 1 - robot pose
Fac.err.J2 = Obs.Jac.E_l; % Jac. of error wrt. node 2 - lmk state

Fac.err.size  = numel(Fac.err.z);

% Append factor to Frame's and Lmk's factors lists.
Frm.factors = [Frm.factors Fac.fac];
Lmk.factors = [Lmk.factors Fac.fac];

% ========== 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