text:
% std_precomp() - Precompute measures (ERP, spectrum, ERSP, ITC) for channels in a study.
%                 If channels are interpolated before computing the measures, the updated
%                 EEG datasets are also saved to disk. Called by pop_precomp(). Follow with
%                 pop_plotstudy(). See Example below.
% Usage:
% >> [ALLEEG,STUDY] = std_precomp(STUDY, ALLEEG, chanorcomp, 'key', 'val', ...);
%
% Required inputs:
%   STUDY        - an EEGLAB STUDY set of loaded EEG structures
%   ALLEEG       - ALLEEG vector of one or more loaded EEG dataset structures
%   chanorcomp   - ['components'|'channels'| or channel cell array] The string
%                  'components' forces the program to precompute all selected
%                  measures for components. The string 'channels' forces the
%                  program to compute all measures for all channels.
%                  A channel cell array containing channel labels will precompute
%                  the selected measures. Note that the name of the channel is
%                  not case-sensitive.
% Optional inputs:
%  'design'   - [integer] use specific study index design to compute
%  measure.
%  'erp'      - ['on'|'off'] pre-compute ERPs for each dataset.
%  'spec'     - ['on'|'off'] pre-compute spectrum for each dataset.
%               Use 'specparams' to set spectrum parameters.
%  'ersp'     - ['on'|'off'] pre-compute ERSP for each dataset.
%               Use 'erspparams' to set time/frequency parameters.
%  'itc'      - ['on'|'off'] pre-compute ITC for each dataset.
%               Use 'erspparams' to set time/frequency parameters.
%  'scalp'    - ['on'|'off'] pre-compute scalp maps for components.
%  'allcomps' - ['on'|'off'] compute ERSP/ITC for all components ('off'
%               only use pre-selected components in the pop_study interface).
%  'specparams' - [cell array] Parameters for the spectopo function are given as
%              optional arguments:
%                   'freqrange' = [min max] frequency range to calculate. Changes x-axis limits {default:
%                                1 Hz for the min and Nyquist (srate/2) for the max. If specified
%                                power distribution maps are plotted, the highest mapped frequency
%                                determines the max freq}. Note that it is better here to compute
%                                spectrum over a wide range of frequencies (it will then
%                                be possible to select another subrange for plotting).
%                   'freqfac'  = [integer] ntimes to oversample -> frequency resolution {default: 2}
%                   'nfft'     = [integer] length to zero-pad data to. Overwrites 'freqfac' above.
%                   'winsize'  = [integer] window size in data points
%                   {default: from data}
%                   'overlap'  = [integer] window overlap in data points {default: 0}
%                   'percent'  = [float 0 to 100] percent of the data to sample for computing the
%                                spectra. Values < 100 speed up the computation. {default: 100}.
%                   'mapnorm'  = [float vector] If 'data' contain the activity of an independant
%                                component, this parameter should contain its scalp map. In this case
%                                the spectrum amplitude will be scaled to component RMS scalp power.
%                                Useful for comparing component strengths {default: none}
%                   'rmdc'     = ['on'|'off'] 'on' -> remove DC {default: 'off'}
%
%              Note that it is advised to compute spectrum
%              over all frequencies since plotting function can always reduce
%              the range of plotted frequencies.
%  'erspparams' - [cell array] Optional arguments are 'cycles', 'freqrange',
%              'padratio', 'winsize', 'alpha' (see newtimef()). Note that it
%              is adivised to select the largest frequency range and time window
%              as plotting function are capable of plotting subranges of
%              these. An important optional parameter that is
%                    'savetrials' = ['on'|'off'] save single-trials ERSP.
%                                   Requires a lot of disk space (dataset
%                                   space on disk times 10) but allow for
%                                   refined single-trial statistics.
%  'recompute' - ['on'|'off'] force recomputing ERP file even if it is
%                already on disk.
%
% Outputs:
%   ALLEEG       - the input ALLEEG vector of EEG dataset structures, modified by adding preprocessing
%                  data as pointers to Matlab files that hold the pre-clustering component measures.
%   STUDY        - the input STUDY set with pre-clustering data added, for use by pop_clust()
%
% Example:
%   >> [ALLEEG STUDY] = std_precomp(STUDY, ALLEEG, { 'cz' 'oz' }, 'interpolate', 'on', 'erp', 'on', ...
%          'spec', 'on', 'ersp', 'on', 'erspparams', { 'cycles' [ 3 0.5 ], 'alpha', 0.01, 'padratio' 1 });
%
%           % This prepares, channels 'cz' and 'oz' in the STUDY datasets.
%           % If a data channel is missing in one dataset, it will be
%           % interpolated (see eeg_interp()). The ERP, spectrum, ERSP, and
%           % ITC for each dataset is then computed.
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, 2006-

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, 2006, arno@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [ STUDY, ALLEEG ] = std_precomp(STUDY, ALLEEG, chanlist, varargin)

    if nargin < 2
        help std_precomp;
        return;
    end;

    if nargin == 2
        chanlist = 'channels'; % default to clustering the whole STUDY
    end
    Ncond = length(STUDY.condition);
    if Ncond == 0
        Ncond = 1;
    end

    g = finputcheck(varargin, { 'erp'         'string'  { 'on' 'off' }     'off';
                                'interp'      'string'  { 'on' 'off' }     'off';
                                'ersp'        'string'  { 'on' 'off' }     'off';
                                'recompute'   'string'  { 'on' 'off' }     'off';
                                'spec'        'string'  { 'on' 'off' }     'off';
                                'scalp'       'string'  { 'on' 'off' }     'off';
                                'allcomps'    'string'  { 'on' 'off' }     'off';
                                'itc'         'string'  { 'on' 'off' }     'off';
                                'savetrials'  'string'  { 'on' 'off' }     'off';
                                'rmicacomps'  'string'  { 'on' 'off' }     'off';
                                'design'      'integer' []                 STUDY.currentdesign;
                                'rmclust'     'integer' []                 [];
                                'rmbase'      'integer' []                 [];
                                'specparams'        'cell'    {}                 {};
                                'erpparams'         'cell'    {}                 {};
                                'erspparams'        'cell'    {}                 {}}, 'std_precomp');
    if isstr(g), error(g); end;
    if ~isempty(g.rmbase), g.erpparams = { g.erpparams{:} 'rmbase' g.rmbase }; end;

    % union of all channel structures
    % -------------------------------
    computewhat = 'channels';
    if isstr(chanlist)
        if strcmpi(chanlist, 'channels')
            chanlist = [];
        else % components
            computewhat = 'components';
            if strcmpi(g.allcomps, 'on')
                chanlist = {};
                for index = 1:length(STUDY.datasetinfo)
                    chanlist = { chanlist{:} [1:size(ALLEEG(STUDY.datasetinfo(index).index).icaweights,1)] };
                end;
            else
                chanlist = { STUDY.datasetinfo.comps };
            end;
        end;
    end;
    if isempty(chanlist)
        alllocs = eeg_mergelocs(ALLEEG.chanlocs);
        chanlist = { alllocs.labels };
    elseif ~isnumeric(chanlist{1})
        alllocs = eeg_mergelocs(ALLEEG.chanlocs);
        [tmp c1 c2] = intersect( lower({ alllocs.labels }), lower(chanlist));
        [tmp c2] = sort(c2);
        alllocs = alllocs(c1(c2));
    end;

    % test if interp and reconstruct channel list
    % -------------------------------------------
    if strcmpi(computewhat, 'channels')
        if strcmpi(g.interp, 'on')
            STUDY.changrp = [];
            STUDY = std_changroup(STUDY, ALLEEG, chanlist, 'interp');
            g.interplocs = alllocs;
        else
            STUDY.changrp = [];
            STUDY = std_changroup(STUDY, ALLEEG, chanlist);
            g.interplocs = struct([]);
        end;
    end;

    % components or channels
    % ----------------------
    if strcmpi(computewhat, 'channels')
         curstruct = STUDY.changrp;
    else curstruct = STUDY.cluster;
    end;

    % compute ERPs
    % ------------
    if strcmpi(g.erp, 'on')
        % check dataset consistency
        % -------------------------
        if length(unique([ALLEEG([STUDY.design(g.design).cell.dataset]).pnts])) > 1
            error([ 'Cannot compute ERPs because datasets' 10 'do not have the same number of data points' ])
        end;

        for index = 1:length(STUDY.design(g.design).cell)
            desset = STUDY.design(g.design).cell(index);
            addopts = { 'savetrials', g.savetrials, 'recompute', g.recompute, 'fileout', desset.filebase, 'trialindices', desset.trials };
            if strcmpi(computewhat, 'channels')
                [tmpchanlist opts] = getchansandopts(STUDY, ALLEEG, chanlist, desset.dataset, g);
                std_erp(ALLEEG(desset.dataset), 'channels', tmpchanlist, opts{:}, addopts{:}, g.erpparams{:});
            else
                if length(desset.dataset)>1 && ~isequal(chanlist{desset.dataset})
                    error(['ICA decompositions must be identical if' 10 'several datasets are concatenated to build' 10 'the design, abording' ]);
                end;
                std_erp(ALLEEG(desset.dataset), 'components', chanlist{desset.dataset(1)}, addopts{:}, g.erpparams{:});
            end;
        end;
        if isfield(curstruct, 'erpdata')
            curstruct = rmfield(curstruct, 'erpdata');
            curstruct = rmfield(curstruct, 'erptimes');
        end;
    end;

    % compute spectrum
    % ----------------
    if strcmpi(g.spec, 'on')
        % check dataset consistency
        % -------------------------
        for index = 1:length(STUDY.design(g.design).cell)
            desset = STUDY.design(g.design).cell(index);
            addopts = { 'savetrials', g.savetrials, 'recompute', g.recompute, 'fileout', desset.filebase, 'trialindices', desset.trials };
            if strcmpi(computewhat, 'channels')
                [tmpchanlist opts] = getchansandopts(STUDY, ALLEEG, chanlist, desset.dataset, g);
                std_spec(ALLEEG(desset.dataset), 'channels', tmpchanlist, opts{:}, addopts{:}, g.specparams{:});
            else
                if length(desset.dataset)>1 && ~isequal(chanlist{desset.dataset})
                    error(['ICA decompositions must be identical if' 10 'several datasets are concatenated to build' 10 'the design, abording' ]);
                end;
                std_spec(ALLEEG(desset.dataset), 'components', chanlist{desset.dataset(1)}, addopts{:}, g.specparams{:});
            end;
        end;
        if isfield(curstruct, 'specdata')
            curstruct = rmfield(curstruct, 'specdata');
            curstruct = rmfield(curstruct, 'specfreqs');
        end;
    end;

    % compute component scalp maps
    % ----------------------------
    if strcmpi(g.scalp, 'on')
        for index = 1:length(STUDY.datasetinfo)

            % find duplicate
            % --------------
            found = [];
            ind1 = STUDY.datasetinfo(index).index;
            inds = strmatch(STUDY.datasetinfo(index).subject, { STUDY.datasetinfo(1:index-1).subject });
            for index2 = inds'
                ind2 = STUDY.datasetinfo(index2).index;
                if isequal(ALLEEG(ind1).icawinv, ALLEEG(ind2).icawinv)
                    found = ind2;
                end;
            end;

            % make link if duplicate
            % ----------------------
            fprintf('Computing/checking topo file for dataset %d\n', ind1);
            if ~isempty(found)
                tmpfile1 = fullfile( ALLEEG(index).filepath, [ ALLEEG(index).filename(1:end-3) 'icatopo' ]);
                tmp.file = fullfile( ALLEEG(found).filepath, [ ALLEEG(found).filename(1:end-3) 'icatopo' ]);
                std_savedat(tmpfile1, tmp);
            else
                std_topo(ALLEEG(index), chanlist{index}, 'none', 'recompute', g.recompute);
            end;
        end;
        if isfield(curstruct, 'topo')
            curstruct = rmfield(curstruct, 'topo');
            curstruct = rmfield(curstruct, 'topox');
            curstruct = rmfield(curstruct, 'topoy');
            curstruct = rmfield(curstruct, 'topoall');
            curstruct = rmfield(curstruct, 'topopol');
        end;
    end;

    % compute ERSP and ITC
    % --------------------
    if strcmpi(g.ersp, 'on') | strcmpi(g.itc, 'on')
        % check dataset consistency
        % -------------------------
        if length(unique([ALLEEG([STUDY.design(g.design).cell.dataset]).pnts])) > 1
            error([ 'Cannot compute ERSPs/ITCs because datasets' 10 'do not have the same number of data points' ])
        end;

        if strcmpi(g.ersp, 'on') & strcmpi(g.itc, 'on'), type = 'both';
        elseif strcmpi(g.ersp, 'on')                   , type = 'ersp';
        else                                             type = 'itc';
        end;

        % check for existing files
        % ------------------------
        guimode = 'guion';
        [ tmpX tmpt tmpf g.erspparams ] = std_ersp(ALLEEG(1), 'channels', 1, 'type', type, 'recompute', 'on', 'getparams', 'on', 'savetrials', g.savetrials, g.erspparams{:});
        if strcmpi(g.recompute, 'off')
            for index = 1:length(STUDY.design(g.design).cell)
                desset = STUDY.design(g.design).cell(index);
                if strcmpi(computewhat, 'channels')
                     filename = [ desset.filebase '.datersp'];
                else filename = [ desset.filebase '.icaersp'];
                end;
                [guimode, g.erspparams] = std_filecheck(filename, g.erspparams, guimode, { 'plotitc' 'plotersp' 'plotphase' });
                if strcmpi(guimode, 'cancel'), return; end;
            end;
            if strcmpi(guimode, 'usedisk') | strcmpi(guimode, 'same'), g.recompute = 'off';
            else                                                       g.recompute = 'on';
            end;
        end;

        % check for existing files
        % ------------------------
        if isempty(g.erspparams),
            tmpparams = {};
        elseif iscell(g.erspparams),
            tmpparams = g.erspparams;
        else
            tmpparams      = fieldnames(g.erspparams); tmpparams = tmpparams';
            tmpparams(2,:) = struct2cell(g.erspparams);
        end;
        tmpparams = { tmpparams{:} 'recompute' g.recompute };
        for index = 1:length(STUDY.design(g.design).cell)
            desset = STUDY.design(g.design).cell(index);
            if strcmpi(computewhat, 'channels')
                [tmpchanlist opts] = getchansandopts(STUDY, ALLEEG, chanlist, desset.dataset, g);
                std_ersp(ALLEEG(desset.dataset), 'channels', tmpchanlist, 'type', type, 'fileout', desset.filebase, 'trialindices', desset.trials, opts{:}, tmpparams{:});
            else
                if length(desset.dataset)>1 && ~isequal(chanlist{desset.dataset})
                    error(['ICA decompositions must be identical if' 10 'several datasets are concatenated to build' 10 'the design, abording' ]);
                end;
                std_ersp(ALLEEG(desset.dataset), 'components', chanlist{desset.dataset(1)}, 'type', type, 'fileout', desset.filebase, 'trialindices', desset.trials, tmpparams{:});
            end;
        end;
        if isfield(curstruct, 'erspdata')
            curstruct = rmfield(curstruct, 'erspdata');
            curstruct = rmfield(curstruct, 'ersptimes');
            curstruct = rmfield(curstruct, 'erspfreqs');
        end;
        if isfield(curstruct, 'itcdata')
            curstruct = rmfield(curstruct, 'itcdata');
            curstruct = rmfield(curstruct, 'itctimes');
            curstruct = rmfield(curstruct, 'itcfreqs');
        end;
    end;

    % components or channels
    % ----------------------
    if strcmpi(computewhat, 'channels')
         STUDY.changrp = curstruct;
    else STUDY.cluster = curstruct;
    end;

    return;

    % find components in cluster for specific dataset
    % -----------------------------------------------
    function rmcomps = getclustcomps(STUDY, rmclust, settmpind);

        rmcomps   = cell(1,length(settmpind));
        for idat = 1:length(settmpind) % scan dataset for which to find component clusters
            currentdat = settmpind(idat);
            for rmi = 1:length(rmclust) % scan clusters
                compinds = STUDY.cluster(rmclust(rmi)).allinds;

                for ind = 1:length(STUDY.cluster(rmclust(rmi)).setinds(:)) % scan datasets of cluster
                    setinds  = [STUDY.design(STUDY.currentdesign).cell(STUDY.cluster(rmclust(rmi)).setinds{ind}).dataset ];
                    indmatch = find(setinds == currentdat);
                    if ~isempty(indmatch)
                        rmcomps{idat} = union( rmcomps{idat}, compinds{ind}(indmatch));
                    end;
                end;
            end;
        end;

    % make option array and channel list (which depend on interp) for any type of measure
    % ----------------------------------------------------------------------
    function [tmpchanlist, opts] = getchansandopts(STUDY, ALLEEG, chanlist, idat, g);

        opts = { };
        if ~isempty(g.rmclust)
            opts = { opts{:} 'rmcomps' getclustcomps(STUDY, g.rmclust, idat) };
        elseif strcmpi(g.rmicacomps, 'on')
            for ind = 1:length(idat)
                rmcomps{ind} = find(ALLEEG(idat(1)).reject.gcompreject);
            end;
            opts = { opts{:} 'rmcomps' rmcomps };
        end;
        if strcmpi(g.interp, 'on')
            tmpchanlist = chanlist;
            allocs = eeg_mergelocs(ALLEEG.chanlocs);
            [tmp ind1 ind2] = intersect({allocs.labels}, chanlist);
            opts = { opts{:} 'interp' allocs(ind1) };
        else
            newchanlist = [];
            tmpchanlocs = ALLEEG(idat(1)).chanlocs;
            chanlocs = { tmpchanlocs.labels };
            for i=1:length(chanlist)
                newchanlist = [ newchanlist strmatch(chanlist(i), chanlocs, 'exact') ];
            end;
            tmpchanlocs =  ALLEEG(idat(1)).chanlocs;
            tmpchanlist = { tmpchanlocs(newchanlist).labels };
        end;

author: PatternRecognition

repo: OpenBMI


text:
classdef prtBrvDiscreteStickBreaking < prtBrv & prtBrvVbOnline

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Properties required by prtAction
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties (SetAccess = private)
        name = 'Discrete Stick Breaking Bayesian Random Variable';
        nameAbbreviation = 'BRVSB';
    end

    properties (SetAccess = protected)
        isSupervised = false;
        isCrossValidateValid = true;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrv
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        function self = estimateParameters(self, x)
            self = conjugateUpdate(self, self, x);
        end

        function y = predictivePdf(self, x)
            y = exp(predictiveLogPdf(self, x));
        end

        function y = predictiveLogPdf(self, x)
            %%%% FIXME
            % The true predictive here is a product of beta-binomials
            % Since that isn't implemented yet we use the average
            % variational loglikelihood

            y = conjugateVariationalAverageLogLikelihood(self, x);
        end

        function val = getNumDimensions(self)
            val = size(self.model.beta,1);
        end

        function self = initialize(self, x)
            x = self.parseInputData(x);
            if ~self.model.isValid
                self.model = self.model.defaultParameters(size(x,2));
            end
        end

        % Optional methods
        %------------------------------------------------------------------
        function kld = conjugateKld(obj, priorObj)
            kld = obj.model.kld(priorObj.model);
        end

        function s = posteriorMeanStruct(obj)
            s.probabilities = exp(obj.model.expectedValueLogProbabilities);
        end

        function val = plotLimits(self)
            val = [0 length(self(1).model.lambda)+1 0 length(self(1).model.lambda)+1];
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrvVb
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods
        function [self, training] = vbBatch(self,x)
            % Since we are purely conjugate we actually don't need vbBatch
            % However we must implement it.
            self = conjugateUpdate(self,x);
            training = struct([]);
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrvVbOnline
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods
        function self = vbOnlineInitialize(self, x) %#ok<INUSD>
            randDraw = rand(1,self.nDimensions);
            randDraw = randDraw./sum(randDraw);

            self = self.conjugateUpdate(self, randDraw);
        end

        function [self, training] = vbOnlineUpdate(self, priorSelf, x, training, prevSelf, learningRate, D) %#ok<INUSL>
            x = self.parseInputData(x);

            [self.model, training] = self.model.vbOnlineWeightedUpdate(priorObj.model, sum(x,1), [], lambda, D, prevObj.model);
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrvMembershipModel
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods
        % Don't actually inherit from prtBrvMembershipModel but these two
        % methods are abstracted there
        function self = conjugateUpdate(self, prior, x)
            x = parseInputData(self,x);
            self.model = self.model.conjugateUpdate(prior.model,x);
        end

        function obj = weightedConjugateUpdate(obj, priorObj, x, weights)
            x = parseInputData(self,x);
            obj.model = obj.model.conjugateUpdate(priorObj.model,bsxfun(@times,x,weights));
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrvVbOnlineMembershipModel
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods
        function [obj, training] = vbOnlineWeightedUpdate(obj, priorObj, x, weights, lambda, D, prevObj)
            x = obj.parseInputData(x);
            if ~isempty(weights)
                x = bsxfun(@times,x,weights);
            end

            [obj.model, training] = obj.model.vbOnlineWeightedUpdate(priorObj.model, sum(x,1), [], lambda, D, prevObj.model);
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Properties for prtBrvDiscreteStickBreaking use
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties
        model = prtBrvDiscreteStickBreakingHierarchy;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods for prtBrvDiscreteStickBreaking use
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods
        function obj = prtBrvDiscreteStickBreaking(varargin)
            if nargin < 1
                return
            end
            obj = obj.constructorInputParse(varargin{:});
        end

        function val = expectedLogMean(obj)
            val = obj.model.expectedValueLogProbabilities(:)';
        end

        function model = modelDraw(obj)
            model.probabilities = draw(obj.model);
        end
    end

    methods (Hidden)
        function x = parseInputData(self,x) %#ok<MANU>
            if isnumeric(x)
                return
            elseif prtUtilIsSubClass(class(x),'prtDataSetBase')
                x = x.getObservations();
            else
                error('prt:prtBrvDiscreteStickBreaking:parseInputData','prtBrvDiscreteStickBreaking requires a prtDataSet or a numeric 2-D matrix');
            end
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function rot = power(rot,n)
% r.^n
%
% Syntax
%   rot = rot^(-1) % inverse rotation
%   rot = rot.^2
%   rot = rot.^[0,1,2,3]
%
%
% Input
%  rot - @rotation
%
% Output
%  rot - @rotation
%
% See also
% rotation/log

rot = power@quaternion(rot,n);

% change inversion
rot.i = (1-(1-2*rot.i).^n)./2;

author: mtex-toolbox

repo: mtex

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


text:
function dz = dynamics(z,u,param)
% dz = dynamics(z,u,param)
%
% Computes the first-order form of the dynamics for the combined chain
% integrator and pendulum system
%
% INPUTS:
%   z = [x;v1;v2;a2;j2];
%   u = [u1;u2];
%
% OUTPUTS:
%   dz = dz/dt
%

x = z(1,:);
v1 = z(2,:);
% v2 = z(3,:);   %Unused
a2 = z(4,:);
j2 = z(5,:);
u1 = u(1,:);
u2 = u(2,:);

% Pendulum physics
dv1 = pendulum(x,v1,u1,param);

% Integrator chain physics:
dx = v1;
dv2 = a2;
da2 = j2;
dj2 = u2;

% Combine:
dz = [dx;dv1;dv2;da2;dj2];

end

author: MatthewPeterKelly

repo: OptimTraj


text:
%By Manas ranjan Das%
function varargout = untitled(varargin)
% UNTITLED M-file for untitled.fig
%      UNTITLED, by itself, creates a new UNTITLED or raises the existing
%      singleton*.
%
%      H = UNTITLED returns the handle to a new UNTITLED or the handle to
%      the existing singleton*.
%
%      UNTITLED('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in UNTITLED.M with the given input arguments.
%
%      UNTITLED('Property','Value',...) creates a new UNTITLED or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before untitled_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to untitled_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help untitled

% Last Modified by GUIDE v7.5 01-Dec-2010 02:35:10

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @untitled_OpeningFcn, ...
                   'gui_OutputFcn',  @untitled_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 untitled is made visible.
function untitled_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 untitled (see VARARGIN)

% Choose default command line output for untitled
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes untitled wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = untitled_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% 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 edit2_Callback(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit2 as text
%        str2double(get(hObject,'String')) returns contents of edit2 as a double

% --- Executes during object creation, after setting all properties.
function edit2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% 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 edit3_Callback(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit3 as text
%        str2double(get(hObject,'String')) returns contents of edit3 as a double

% --- Executes during object creation, after setting all properties.
function edit3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% 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

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
freq=str2double(get(handles.edit1,'String'));
er=str2double(get(handles.edit2,'String'));
h=str2double(get(handles.edit3,'String'));
height = h/1000*2.54;
Width=30.0/(2.0*freq)*sqrt(2.0/(er+1.0))*10;
ereff=(er+1.0)/2.0+(er-1)/(2.0*sqrt(1.0+12.0*height/Width));
dl=0.412*height*((ereff+0.3)*(Width/height+0.264))/((ereff-0.258)*(Width/...
    height+0.8));
lamda=30.0/(freq*sqrt(ereff));
lamda0=30.0/freq;
Leff=lamda/2;
L=Leff-2*dl;
Length = L*10;
k0=2*pi*freq/30;
phi1=0:360;
theta1=0:180;
phi=phi1./180.*pi;
theta=theta1./180.*pi;

Rr=120*lamda0/(1-k0^2*height^2/24);
Radiation_resistance = Rr;

set(handles.text4,'String',Length);
set(handles.text5,'String',Width);
set(handles.text11,'String',Radiation_resistance);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Radiation patterns
% Plotting of radiation patterns on polar and rectangular co-ordinate  system

% Calculating the normalized field values for polar plot
Etheta=sin(k0*height/2.*cos(theta)).*cos(k0*L/2.*cos(theta))/k0/height*...
    2./cos(theta);
Ethetamax=max(Etheta);
Ethetanor=Etheta./Ethetamax;

Ephi=sin(k0*Width/2.*cos(phi)).*sin(phi)/k0/Width*2./cos(phi);
Ephimax=max(Ephi);
Ephinor=Ephi./Ephimax;
axes(handles.axes1);
cla;

polar(theta,Ethetanor,'-r')
hold on;
polar(phi,Ephinor,'-b')
title('Radiation plot of E and H plane patterns');
% title('E- and H-plane Patterns of Rectangular Microstrip Antenna',...
%     'fontsize',[12]);
legend('E plane','H plane');

% Space wave, surface wave and gain calculations

k0d=k0*height;

%Space wave power
Psp = 377*k0^2*k0d^2/3/pi*(1-1/ereff+2/5/ereff^2);

%Surface wave power
s=sqrt(ereff-1);
k0ds=k0d*s;
alpha0=s*tan(k0ds);
alpha1=-[tan(k0ds)+k0ds/(cos(k0ds))^2];
x0d=ereff^2-alpha1^2;
x0n=-ereff^2+alpha0*alpha1+ereff*sqrt(ereff^2-2*alpha0*alpha1+alpha0^2);
x0 = 1+x0n/x0d;

Psur = 377*k0^2/4*ereff*(x0^2-1)/[ereff*[1/sqrt(x0^2-1)+sqrt(x0^2-1)/(...
    ereff-x0^2)]+k0d*[1+ereff^2*(x0^2-1)/(ereff-x0^2)]];

Efficiency=Psp/(Psp+Psur);

Gain=10*log10(Efficiency*2*pi*Leff*Width/lamda0);

set(handles.text14,'String',Efficiency);
set(handles.text15,'String',Gain);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% recessed feedline width
% characteristic impedance

x=30*pi/sqrt(er)/50-0.441;

if sqrt(er)*50 <= 120
    W0=height*x;
else
    W0=height*(0.85-sqrt(0.6-x));
end
% Width_of_microstrip_feed=W0;

if W0/height <= 1
    Z0=60/sqrt(ereff)*log(8*height/W0+W0/4/height);
else
    Z0=120*pi/sqrt(ereff)*inv(W0/height+1.393+.667*log(W0/height+1.444));
end
% characteristic_Impedance = Z0
set(handles.text17,'String',Z0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes during object creation, after setting all properties.
function text1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to text1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

author: Sable

repo: mcbench-benchmarks


text:
function [n_signif,index_signif] = fdr(pvals,qlevel,method,adjustment_method,adjustment_args)
% Determine significance based on a false discovery rate (FDR) method
%
%  [n_signif,index_signif] = ...
%     fdr(pvals,qlevel,method,adjustment_method,adjustment_args)
%
% This is the main function designed for general usage for determining
% significance based on the FDR approach. The specific method used depends
% on the calling parameters.
%
% Arguments:
%
%   pvals:  a vector of pvals on which to conduct the multiple testing
%
%   qlevel: the proportion of false positives desired
%
%   method: method for performing the testing.  'original' follows
%   Benjamini & Hochberg (1995); 'general' is much more
%   conservative, requiring no assumptions on the p-values (see
%   Benjamini & Yekutieli (2001)).  We recommend using 'original',
%   and if desired, using 'adjustment.method="mean" ' to increase
%   power.
%
%   adjustment_method: method for increasing the power of the
%   procedure by estimating the proportion of alternative p-values,
%   one of 'mean', the modified Storey estimator that we suggest in
%   Ventura et al. (2004), 'storey', the method of Storey (2002),
%   or 'two-stage', the iterative approach of Benjamini et
%   al. (2001)
%
%   adjustment_args: arguments to adjustment.method; see prop_alt()
%   for description, but note that for 'two-stage', qlevel and
%   fdr_method are taken from the qlevel and method arguments to
%   fdr()
%
% Value:
%
%   index_signif:  a vector of the indices of the significant tests
%   or NA if no significant tests
%
%   n_signif: number of significant tests
%
% Examples:
%
%   signif <- fdr(pvals,method='original',adjustment.method='mean')
%
%  References:
%             Ventura, V., C.J. Paciorek, and J.S. Risbey.  2004.
%             Controlling the proportion of falsely-rejected
%             hypotheses when conducting multiple tests with
%             climatological data.  Journal of Climate, in press.
%             Also Carnegie Mellon University, Department of
%             Statistics technical report 775
%             (www.stat.cmu.edu/tr/tr775/tr775.html).
%
%             Benjamini, Y, and Y. Hochberg. 1995. Controlling the
%             false discovery rate: a practical and powerful
%             approach to multiple testing.  JRSSB 57:289-300.
%
%             Benjamini, Y. and D. Yekutieli.  2001.  The control
%             of the false discovery rate in multiple testing under
%             dependency. Annals of Statistics 29:1165-1188.
%
%             Benjamini, Y., A. Krieger, and D. Yekutieli.  2001.
%             Two staged linear step up FDR controlling procedure.
%             Technical Report, Department of Statistics and
%             Operations Research, Tel Aviv University.  URL:
%             http://www.math.tau.ac.il/~ybenja/Papers.html
%
%             Storey, J. 2002.  A direct approach to false discovery rates.  JRSSB 64: 479-498.
%
% Attribution:
%
%    Version 0.1.3;  May 10, 2004
%
%    License: GPL version 2 or later
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
%
%    Author: Chris Paciorek - please contact the author with bug
%    reports: paciorek AT alumni.cmu.edu
%

if exist('isfield') % matlab
    endline='';
else % octave
    endline='\n';
end

n = length(pvals);

a = 0;   % initialize proportion of alternative hypotheses
if nargin<2
    error(['pvals and qlevel are required arguments',endline]);
end
if nargin==2
    method='original';
end
if nargin>=4 && ~isempty(adjustment_method)
    if strcmp(adjustment_method,'two-stage')  % set things up for the "two-stage" estimator
        qlevel = qlevel/(1+qlevel);   % see Benjamini et al. (2001) for proof that this controls the FDR at level qlevel
        adjustment_args.qlevel = qlevel;
        adjustment_args.fdr_method = method;
        %fprintf(['Adjusting cutoff using two-stage approach with fdr ', ...
        %    'method %s and qlevel %.3f\n'],adjustment_args.fdr_method,...
        %    adjustment_args.qlevel);
    end
    if strcmp(adjustment_method,'mean') & nargin==4  % default arguments for "mean" method of Ventura et al. (2004)
        adjustment_args.edf_lower=0.8;
        adjustment_args.num_steps=20;
        %fprintf(['Adjusting cutoff using mean approach with edf_lower=0.8 ', ...
        %    'and num_steps=20\n']);
    end
    if strcmp(adjustment_method,'storey') & nargin==4
        error(['adjustment_args.edf_quantile not specified for Storey ', ...
            'adjustment',endline]);
    end
    a = prop_alt(pvals,adjustment_method,adjustment_args);
end
if a==1    % all hypotheses are estimated to be alternatives
    index_signif=1:n;
    n_signif=n;
else      % adjust for estimate of a; default is 0
    qlevel = qlevel/(1-a);
    [n_signif,index_signif]=fdr_master(pvals,qlevel,method);
end

function [n_signif,index_signif] = fdr_master(pvals,qlevel,method)
%
% Description:
%
%    This is an internal function that performs various versions of
%    the FDR procedure, but without the modification described in
%    section 4 of Ventura et al. (2004)
%
% Arguments:
%
%   pvals (required):  a vector of pvals on which to conduct the multiple testing
%
%   qlevel: the proportion of false positives desired
%
%   method: one of 'original', the original method of Benjamini &
%   Hochberg (1995), or 'general', the method of Benjamini &
%   Yekutieli (2001), which requires no assumptions about the
%   p-values, but which is much more conservative.  We recommend
%   'original' for climatological data, and suspect it works well
%   generally for spatial data.
%
% Value:
%
%   index_signif:  a vector of the indices of the significant tests
%   or NA if no significant tests
%
%   n_signif: number of significant tests
%
%
if nargin==2
    method='original';
end
n = length(pvals);
if strcmp(method,'general')
    qlevel = qlevel/sum(1./(1:n));  % This requires fewer assumptions but is much more conservative
else
    if(~strcmp(method,'original'))
        error('No method of type: %s\n',method);
    end
end
[n_signif,index_signif] = fdr_basic(pvals,qlevel);

function [n_signif,index_signif] = fdr_basic(pvals,qlevel)
%
% Description:
%
%    This is an internal function that performs the basic FDR of Benjamini & Hochberg (1995).
%
% Arguments:
%
%   pvals:  a vector of pvals on which to conduct the multiple testing
%
%   qlevel: the proportion of false positives desired
%
% Value:
%
%   index_signif:  a vector of the indices of the significant tests
%   or NA if no significant tests
%
%   n_signif: number of significant tests
%
n = length(pvals);
if n>1 & size(pvals,2)==1
    pvals=pvals';
end
[sorted_pvals,sort_index] = sort(pvals);
indices = (1:n).*(sorted_pvals<=qlevel*(1:n)/n);
n_signif = max(indices);
if(n_signif)
    indices = 1:n_signif;
    index_signif=sort(sort_index(indices));
else
    index_signif=nan;
end

function [a] = storey(edf_quantile,pvals)
%
% Description:
%
%    This is an internal function that calculates the basic Storey
%    (2002) estimator of a, the proportion of alternative
%    hypotheses.
%
% Arguments:
%
%   edf_quantile (required):  the quantile of the empirical
%   distribution function at which to estimate a
%
%   pvals (required):  a vector of pvals on which to estimate a
%
% Value:
%
%   a: estimate of a, the number of alternative hypotheses
%
%
if exist('isfield') % matlab
    endline='';
else % octave
    endline='\n';
end
if edf_quantile >= 1 | edf_quantile <=0
    error(['edf_quantile should be between 0 and 1',endline]);
end
a = (mean(pvals<=edf_quantile)-edf_quantile)/(1-edf_quantile);
if(a<0)
    a=0;
end

function [a] = prop_alt(pvals,adjustment_method,adjustment_args)
%
% Description:
%
%    This is an internal function that calculates an estimate of a,
%    the proportion of alternative hypotheses, using one of several
%    methods.
%
% Arguments:
%
%   pvals (required):  a vector of pvals from which to estimate a
%
%   adjustment_method: method for  estimating the proportion of
%   alternative p-values, one of "mean", the modified Storey
%   estimator suggested in Ventura et al. (2004); "storey", the
%   method of Storey (2002); or "two-stage", the iterative approach
%   of Benjamini et al. (2001)
%
%   adjustment_args: arguments to adjustment.method;
%
%      for 'mean', specify edf_lower, the smallest quantile at
%      which to estimate a, and num_steps, the number of quantiles
%      to use - the approach uses the average of the Storey (2002)
%      estimator for the num.steps quantiles starting at edf.lower
%      and finishing just less than 1
%
%      for 'storey', specify edf_quantile, the quantile at which to
%      calculate the estimator
%
%      for 'two-stage', the method uses a standard FDR approach to
%      estimate which p-values are significant; this number is the
%      estimate of a; therefore the method requires specification
%      of qlevel, the proportion of false positives and fdr.method
%      ('original' or 'general'), the FDR method to be used.  We do
%      not recommend 'general' as this is very conservative and
%      will underestimate a.

%
% Value:
%
%   a:  estimate of a, the number of alternative hypotheses
%
% Examples:
%
%   a = prop_alt(pvals,'mean')
%
if exist('isfield') % matlab
    endline='';
else % octave
    endline='\n';
end

stop=0;
n = length(pvals);
if strcmp(adjustment_method,'two-stage')
    if ~exist('adjustment_args')
        stop=1;
    else
        if exist('isfield')  % check if in Matlab
            if ~(isfield(adjustment_args,'qlevel') & ...
                    isfield(adjustment_args,'fdr_method'))
                stop=1;
            end
        else  % in octave
            if ~(struct_contains(adjustment_args,'qlevel') & ...
                    struct_contains(adjustment_args,'fdr_method'))
                stop=1;
            end
        end
    end
    if stop
        error(['adjustment_args.qlevel or adjustment_args.fdr_method ', ...
            'not specified.  Two-stage estimation of the number ',...
            'of alternative hypotheses requires specification of ', ...
            'the  FDR threshold and FDR method',endline]);
    end
    [n_signif,index_signif]=fdr_master(pvals,adjustment_args.qlevel,...
        adjustment_args.fdr_method);
    a=n_signif/n;
end
if strcmp(adjustment_method,'storey')
    if ~(exist('adjustment_args'))
        stop=1;
    else
        if exist('isfield')  % check if in Matlab
            if ~(isfield(adjustment_args,'edf_quantile'))
                stop=1;
            end
        else  % octave
            if ~(struct_contains(adjustment_args,'edf_quantile'))
                stop=1;
            end
        end
    end
    if stop
        error(['adjustment_args.edf_quantile not specified. Using method of Storey for estimating  the number of alternative hypotheses requires specification of the argument of the p-value EDF at which to do the estimation (a number close to one is recommended)',endline]);
    end
    a=storey(adjustment_args.edf_quantile,pvals);
end

if strcmp(adjustment_method,'mean')
    if ~(exist('adjustment_args'))
        stop=1;
    else
        if exist('isfield')  % check if in Matlab
            if ~(isfield(adjustment_args,'edf_lower') & ...
                    isfield(adjustment_args,'num_steps'))
                stop=1;
            end
        else % octave
            if ~(struct_contains(adjustment_args,'edf_lower') & ...
                    struct_contains(adjustment_args,'num_steps'))
                stop=1;
            end
        end
    end
    if stop
        error(['adjustment_args.edf_lower or adjustment_args.num_steps ' ...
            'is not specified. Using the method of Ventura et al. (2004) for estimating  the number of alternative hypotheses requires specification of the lowest quantile of the p-value EDF at which to do the estimation (a number close to one is recommended) and the number of steps between edf_lower and 1, starting at edf_lower, at which to do the estimation',endline])
    end
    if adjustment_args.edf_lower >=1 | adjustment_args.edf_lower<=0
        error(['adjustment_args.edf_lower must be between 0 and 1',endline]);
    end
    if adjustment_args.num_steps<1 | ~(rem(adjustment_args.num_steps,1)==0)
        error(['adjustment_args.num_steps must be an integer greater than 0',endline]);
    end
    stepsize = (1-adjustment_args.edf_lower)/adjustment_args.num_steps;
    edf_quantile=adjustment_args.edf_lower;
    a_vec=storey(edf_quantile,pvals);
    for step=1:adjustment_args.num_steps
        a_vec(step)=storey(edf_quantile,pvals);
        edf_quantile=edf_quantile+stepsize;
    end
    a=mean(a_vec);
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
fname=input('Input File (in ASCII format)? ','s');
hash_foutname=input('Output File for md5 Hash? ','s');
time1=clock;
%Open the input file and get the first line of data
fid=fopen(fname);
M = fread(fid);
fclose(fid);

pwd='sri';
ini=reshape(dec2bin(pwd,8),1,24);

for ii = 2:length(M)
    ini=cat(2,ini,dec2bin(M(ii),8));
end

s2=length(ini)/8;

block_temp = [ ini, ...
                 '1', ...
                 num2str(zeros(mod(448-1-s2*8,512),1))' ...
                dec2bin(s2*8,64) ];
            nb=length(block_temp)/512;
block = reshape(block_temp,512,nb)';

shi=[7,12,17,22,7,12,17,22,7,12,17,22,7,12,17,22,5,9,14,20,5,9,14,20,5,9,14,20,5,9,14,20,4,11,16,23,4,11,16,23,4,11,16,23,4,11,16,23,6,10,15,21,6,10,15,21,6,10,15,21,6,10,15,21];

mp=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,2,7,12,1,6,11,16,5,10,15,4,9,14,3,8,13,6,9,12,15,2,5,8,11,14,1,4,7,10,13,16,3,1,8,15,6,13,4,11,2,9,16,7,14,5,12,3,10];

a=dec2bin(hex2dec('01234567'),32);a1=a;
b=dec2bin(hex2dec('89abcdef'),32);b1=b;
c=dec2bin(hex2dec('fedcba98'),32);c1=c;
d=dec2bin(hex2dec('76543210'),32);d1=d;
for b=1:nb

inp=reshape(block(b,:),32,16)';

%disp('round1');
%Round1

for i=1:64

    a=d;
    d=c;
    c=b;
if i>48
t1= bin2dec2(xor(c,(b&~d)));
elseif i>32
t1= bin2dec2(xor(b,xor(c,d)));
elseif i>16
t1= bin2dec2((b&d)|(c&~d));
else
t1= bin2dec2((b&c)|(~b&d));
end

t2= floor(2^32*abs(sin(i)));
t3= bin2dec2(inp(mp(i),:));
t4=bin2dec2(b);
anst=dec2bin(mod(t1+t2+t3,2^32),32);
ans=cls(anst,shi(i));
b=dec2bin(mod(bin2dec2(ans)+t4,2^32),32);

end

a=dec2bin(mod(bin2dec2(a)+bin2dec2(a1),2^32),32);
b=dec2bin(mod(bin2dec2(b)+bin2dec2(b1),2^32),32);
c=dec2bin(mod(bin2dec2(c)+bin2dec2(c1),2^32),32);
d=dec2bin(mod(bin2dec2(d)+bin2dec2(d1),2^32),32);

end
a=dec2hex(bin2dec(a),8);disp(a);
b=dec2hex(bin2dec(b),8);disp(b);
c=dec2hex(bin2dec(c),8);disp(c);
d=dec2hex(bin2dec(d),8);disp(d);

time2=clock;

disp(etime(time2,time1));

fid=fopen(hash_foutname,'w+');
fprintf(fid,'%s',a);
fprintf(fid,'%s',b);
fprintf(fid,'%s',c);
fprintf(fid,'%s',d);

fclose(fid);


text:
%% DEMO_febio_0002_beam_force
% Below is a demonstration for:
%
% * Building geometry for a beam with hexahedral elements
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement results

%% Keywords
%
% * febio_spec version 4.0
% * febio, FEBio
% * beam force loading
% * force control boundary condition
% * hexahedral elements, hex8, hex20
% * beam, rectangular
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=20;
faceAlpha1=0.8;
markerSize=40;
markerSize2=20;
lineWidth=3;

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[febioFebFileNamePart,'.txt']; %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement

%Specifying dimensions and number of elements
beamWidth=10;
sampleWidth=beamWidth; %Width
sampleThickness=4*beamWidth; %Thickness
sampleHeight=beamWidth; %Height
pointSpacings=3*ones(1,3); %Desired point spacing between nodes
numElementsWidth=round(sampleWidth/pointSpacings(1)); %Number of elemens in dir 1
numElementsThickness=round(sampleThickness/pointSpacings(2)); %Number of elemens in dir 2
numElementsHeight=round(sampleHeight/pointSpacings(3)); %Number of elemens in dir 3

elementType='hex20'; %'hex8' or 'hex20'

%Define applied force
appliedForce=[0 0 -2e-3];

%Material parameter set
c1=1e-3; %Shear-modulus-like parameter
m1=8; %Material parameter setting degree of non-linearity
k_factor=1e2; %Bulk modulus factor
k=c1*k_factor; %Bulk modulus

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=6; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size

runMode='external';% 'internal' or 'external'

%% Creating model geometry and mesh
% A box is created with tri-linear hexahedral (hex8) elements using the
% |hexMeshBox| function. The function offers the boundary faces with
% seperate labels for the top, bottom, left, right, front, and back sides.
% As such these can be used to define boundary conditions on the exterior.

% Create a box with hexahedral elements
beamDimensions=[sampleWidth sampleThickness sampleHeight]; %Dimensions
beamElementNumbers=[numElementsWidth numElementsThickness numElementsHeight]; %Number of elements
outputStructType=2; %A structure compatible with mesh view
[meshStruct]=hexMeshBox(beamDimensions,beamElementNumbers,outputStructType);

%Access elements, nodes, and faces from the structure
E=meshStruct.elements; %The elements
V=meshStruct.nodes; %The nodes (vertices)
Fb=meshStruct.facesBoundary; %The boundary faces
Cb=meshStruct.boundaryMarker; %The "colors" or labels for the boundary faces
elementMaterialIndices=ones(size(E,1),1); %Element material indices

if strcmp(elementType,'hex20')
    [E,V,~,Fb]=hex8_hex20(E,V,{},Fb);
    meshStruct.elements=E;
    meshStruct.nodes=V;
    meshStruct.Fb=Fb;
end

%%
% Plotting model boundary surfaces and a cut view

hFig=cFigure;

subplot(1,2,1); hold on;
title('Model boundary surfaces and labels','FontSize',fontSize);
hp=gpatch(Fb,V,Cb,'k',faceAlpha1);
hp.Marker='.';
hp.MarkerSize=markerSize2;

colormap(gjet(6)); icolorbar;
axisGeom(gca,fontSize);

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
meshView(meshStruct,optionStruct);
axisGeom(gca,fontSize);

drawnow;

%% Defining the boundary conditions
% The visualization of the model boundary shows colors for each side of the
% cube. These labels can be used to define boundary conditions.

%Define supported node set
bcSupportList=unique(Fb(Cb==4,:)); %Node set part of selected face

%Prescribed force nodes
bcPrescribeList=unique(Fb(Cb==3,:)); %Node set part of selected face

%%
% Visualizing boundary conditions. Markers plotted on the semi-transparent
% model denote the nodes in the various boundary condition lists.

hf=cFigure;
title('Boundary conditions','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fb,V,'kw','k',0.5);

hl(1)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize);
hl(2)=plotV(V(bcPrescribeList,:),'r.','MarkerSize',markerSize);

legend(hl,{'BC support','BC prescribe'});

axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='4.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.qn_method.max_ups=max_ups;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1;
febio_spec.Material.material{1}.m1=m1;
febio_spec.Material.material{1}.c2=c1;
febio_spec.Material.material{1}.m2=-m1;
febio_spec.Material.material{1}.k=k;

%Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type=elementType; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E; %The element matrix

% -> NodeSets
nodeSetName1='bcSupportList';
febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.VAL=bcSupportList(:)';

nodeSetName2='bcPrescribeList';
febio_spec.Mesh.NodeSet{2}.ATTR.name=nodeSetName2;
febio_spec.Mesh.NodeSet{2}.VAL=bcPrescribeList(:)';

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.name='zero_displacement_xyz';
febio_spec.Boundary.bc{1}.ATTR.type='zero displacement';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.x_dof=1;
febio_spec.Boundary.bc{1}.y_dof=1;
febio_spec.Boundary.bc{1}.z_dof=1;

%Loads section
% -> Prescribed nodal forces
nodalLoadType='force';
switch nodalLoadType
    case 'load' %Apply same force to each node
        febio_spec.Loads.nodal_load{1}.ATTR.name='PrescribedForceX';
        febio_spec.Loads.nodal_load{1}.ATTR.type='nodal_load';
        febio_spec.Loads.nodal_load{1}.ATTR.node_set=nodeSetName2;
        febio_spec.Loads.nodal_load{1}.dof='x';
        febio_spec.Loads.nodal_load{1}.scale.ATTR.lc=1;
        febio_spec.Loads.nodal_load{1}.scale.VAL=appliedForce(1)/numel(bcPrescribeList);

        febio_spec.Loads.nodal_load{2}.ATTR.name='PrescribedForceY';
        febio_spec.Loads.nodal_load{2}.ATTR.type='nodal_load';
        febio_spec.Loads.nodal_load{2}.ATTR.node_set=nodeSetName2;
        febio_spec.Loads.nodal_load{2}.dof='y';
        febio_spec.Loads.nodal_load{2}.scale.ATTR.lc=1;
        febio_spec.Loads.nodal_load{2}.scale.VAL=appliedForce(2)/numel(bcPrescribeList);

        febio_spec.Loads.nodal_load{3}.ATTR.name='PrescribedForceZ';
        febio_spec.Loads.nodal_load{3}.ATTR.type='nodal_load';
        febio_spec.Loads.nodal_load{3}.ATTR.node_set=nodeSetName2;
        febio_spec.Loads.nodal_load{3}.dof='z';
        febio_spec.Loads.nodal_load{3}.scale.ATTR.lc=1;
        febio_spec.Loads.nodal_load{3}.scale.VAL=appliedForce(3)/numel(bcPrescribeList);
    case 'force' %Apply a force vector to a collection of nodes (destributed)
        febio_spec.Loads.nodal_load{1}.ATTR.name='PrescribedForceX';
        febio_spec.Loads.nodal_load{1}.ATTR.type='nodal_force';
        febio_spec.Loads.nodal_load{1}.ATTR.node_set=nodeSetName2;
        febio_spec.Loads.nodal_load{1}.value.ATTR.lc=1;
        febio_spec.Loads.nodal_load{1}.value.VAL=appliedForce/numel(bcPrescribeList);
end

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.name='LC_1';
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
%febio_spec.LoadData.load_controller{1}.extend='CONSTANT';
febio_spec.LoadData.load_controller{1}.points.pt.VAL=[0 0; 1 1];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode=runMode;

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    %%
    % Importing nodal displacements from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),1,1);

    %Access data
    N_disp_mat=dataStruct.data; %Displacement
    timeVec=dataStruct.time; %Time

    %Create deformed coordinate set
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    DN_magnitude=sqrt(sum(N_disp_mat(:,:,end).^2,2)); %Current displacement magnitude

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('Displacement magnitude [mm]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),DN_magnitude,'k',1); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    gpatch(Fb,V,0.5*ones(1,3),'k',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([0 max(DN_magnitude)]);
    axis(axisLim(V_DEF)); %Set axis limits statically
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN_magnitude=sqrt(sum(N_disp_mat(:,:,qt).^2,2)); %Current displacement magnitude

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),DN_magnitude}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

end

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

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [A,R,w_jk,mu_jk,sigma_jk,extra] = gmm_hu_ex(I, endmembers, options)
%GMM_EX Summary of this function goes here
%   Detailed explanation goes here
if nargin < 3
    options = [];
end

if max(I(:)) > 10 || max(endmembers{1}(:)) > 10
    disp('Warning! I or endmember spectra are not in the range 0 - 1.');
end

if class(I) ~= class(endmembers{1})
    disp('Error! I and endmember spectra do not have the same data type');
end

options = insert_param_when_absent(options, 'reduced_dim', 10);
options = insert_param_when_absent(options, 'max_num_comp', 4);
options = insert_param_when_absent(options, 'show_fig', 0);

options = insert_param_when_absent(options, 'fix_w_k', 1);
options = insert_param_when_absent(options, 'fix_mu_jk', 1);
options = insert_param_when_absent(options, 'fix_sigma_jk', 1);

use_predefined_projection_components = parse_param(options, ...
    'use_predefined_projection_components', 0);

[Y_ori,~,rows,cols] = reshape_hsi(I,[]);
M = length(endmembers);

if use_predefined_projection_components
    mapping = options.project_mapping;
    Y = gmm_project(Y_ori, mapping);
    components = options.predefined_components;
    K = components.K;
    w_jk = components.w_jk;
    mu_jk = components.mu_jk;
    sigma_jk = components.sigma_jk;
else
    % project_mode = 'image' or 'endmembers'
    project_mode = parse_param(options,'project_mode','endmembers');

    if strcmp(project_mode, 'image')
        [Y, mapping] = pca(Y_ori, options.reduced_dim);
    elseif strcmp(project_mode, 'endmembers')
        mapping = calc_projection_from_library(endmembers, options.reduced_dim);
        Y = gmm_project(Y_ori, mapping);
    elseif strcmp(project_mode, 'combined')
        % Gives pretty bad results, should figure out a new implemenation
        mapping = calc_projection_from_image_library(Y_ori, endmembers, options.reduced_dim);
        Y = gmm_project(Y_ori, mapping);
    elseif strcmp(project_mode, 'custom')
        mapping = options.project_mapping;
        Y = gmm_project(Y_ori, mapping);
    else
        disp('ERROR: Unknown projection mode');
    end
    options = insert_param_when_absent(options, 'project_mapping', mapping);

    % Y = gmm_project(Y_ori, mapping);
    [K,w_jk,mu_jk,sigma_jk] = estimate_components(endmembers, mapping, options);
end

if options.show_fig
    names = options.names;
    endmember_scatter_plot_end_var(Y,w_jk,mu_jk,sigma_jk,names);
end

options.w_jk = w_jk;
options.mu_jk = mu_jk;
options.sigma_jk = sigma_jk;
options.K = K;

options = insert_param_when_absent(options, 'beta1', 0);
options = insert_param_when_absent(options, 'beta2', 0);
% options = insert_param_when_absent(options, 'convergence_thresh', 0.002);

% no need to project back to the original dimension, gmm_hu will do it.
[A,R,w_jk,mu_jk,sigma_jk,extra] = gmm_hu(I,M,options);

function mapping = calc_projection_from_image_library(Y_ori, endmembers, reduced_dim)
spectra_train = cell2mat(endmembers');
if size(spectra_train,1) > size(Y_ori,1)
    Y_ori = repmat(Y_ori, round(size(spectra_train,1) / size(Y_ori,1)), 1);
end
spectra_train = cat(1, spectra_train, Y_ori);
[~,mapping] = pca(spectra_train, reduced_dim);

function [beta1,beta2] = automatic_params(I,M,options)
init_beta1 = 10;
init_beta2 = 10;

options.fix_w_k = 0;
options.convergence_thresh = 0.1;

w_jk = options.w_jk;

beta1s = 10.^(-5:5)';
beta2s = 10.^(-5:5)';

errs = Inf(length(beta1s), length(beta2s));
% A_errs = Inf(length(beta1s), length(beta2s));

[~,ind1] = min(abs(beta1s - init_beta1));
[~,ind2] = min(abs(beta2s - init_beta2));
while 1
    ind1s = [ind1-1, ind1, ind1+1];
    ind2s = [ind2-1, ind2, ind2+1];
%     ind2s = ind2;
    for i = ind1s
        for j = ind2s
            % only the inner errors will be updated so the current point
            % will not move to the boundary
            if i > 1 && i < length(beta1s) && j > 1 && ...
                    j < length(beta2s) && errs(i,j) == Inf
                options.beta1 = beta1s(i);
                options.beta2 = beta2s(j);
                [A,R,w_jk1,mu_jk,sigma_jk] = gmm_hu(I,M,options);
                title = ['beta1: ',num2str(options.beta1),', beta2: ',num2str(options.beta2)];
                show_abundances(A, rows, cols, title);
%                 A_errs(i,j) = calc_abundance_error(A_gt,A);
                errs(i,j) = mdiff(w_jk,w_jk1);
            end
        end
    end
    if errs(ind1,ind2) <= errs(ind1-1,ind2) && errs(ind1,ind2) <= errs(ind1+1,ind2) ...
            && errs(ind1,ind2) <= errs(ind1,ind2-1) && errs(ind1,ind2) <= errs(ind1,ind2+1)
        break;
    else
        errs1 = errs(ind1-1:ind1+1,ind2-1:ind2+1);
        [~,ind] = min(errs1(:));
        [y,x] = ind2sub([3 3], ind);
        ind1 = ind1 + y - 2;
        ind2 = ind2 + x - 2;
    end
end

beta1 = beta1s(ind1);
beta2 = beta2s(ind2);

author: zhouyuanzxcv

repo: Hyperspectral


text:
function x = box_behnken ( dim_num, x_num, range )

%*****************************************************************************80
%
%% BOX_BEHNKEN returns a Box-Behnken design for the given number of factors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    George Box, Donald Behnken,
%    Some new three level designs for the study of quantitative variables,
%    Technometrics,
%    Volume 2, pages 455-475, 1960.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer X_NUM, the number of elements of the design.
%    X_NUM should be equal to DIM_NUM * 2**(DIM_NUM-1) + 1.
%
%    Input, real RANGE(DIM_NUM,2), the minimum and maximum
%    value for each component.
%
%    Output, real X(DIM_NUM,X_NUM), the elements of the design.
%

%
%  Ensure that the range is legal.
%
  if ( any ( range(1:dim_num,2) <= range(1:dim_num,1) ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BOX_BEHNKEN - Fatal error!\n' );
    fprintf ( 1, '  For some index I,\n' );
    fprintf ( 1, '  RANGE(I,2) <= RANGE(I,1).\n' );
    error ( 'BOX_BEHNKEN - Fatal error!' );
  end
%
%  The first point is the center.
%
  j = 1;

  x(1:dim_num,j) = ( range(1:dim_num,1) + range(1:dim_num,2) ) / 2.0;
%
%  For subsequent elements, one entry is fixed at the middle of the range.
%  The others are set to either extreme.
%
  for i = 1 : dim_num

    j = j + 1;

    x(1:dim_num,j) = range(1:dim_num,1);
    x(i,j) = ( range(i,1) + range(i,2) ) / 2.0;
%
%  The next element is made by finding the last low value, making it
%  high, and all subsequent high values low.
%
    while ( 1 )

      last_low = -1;

      for i2 = 1 : dim_num
        if ( x(i2,j) == range(i2,1) )
          last_low = i2;
        end
      end

      if ( last_low == -1 )
        break
      end

      j = j + 1;
      x(1:dim_num,j) = x(1:dim_num,j-1);
      x(last_low,j) = range(last_low,2);

      for i2 = last_low + 1 : dim_num
        if ( x(i2,j) == range(i2,2) )
          x(i2,j) = range(i2,1);
        end
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Image super-resolution using sparse representation
% Example code
%
% Nov. 2, 2007. Jianchao Yang
% IFP @ UIUC
%
% Revised version. April, 2009.
%
% Reference
% Jianchao Yang, John Wright, Thomas Huang and Yi Ma. Image superresolution
% via sparse representation of raw image patches. IEEE Computer Society
% Conference on Computer Vision and Pattern Recognition (CVPR), 2008.
%
% For any questions, email me by jyang29@illinois.edu

clear all;
clc;

addpath('Solver');
addpath('Sparse coding');

% =====================================================================
% specify the parameter settings

patch_size = 3; % patch size for the low resolution input image
overlap = 1; % overlap between adjacent patches
lambda = 0.1; % sparsity parameter
zooming = 3; % zooming factor, if you change this, the dictionary needs to be retrained.

tr_dir = 'Data/training'; % path for training images
skip_smp_training = true; % sample training patches
skip_dictionary_training = true; % train the coupled dictionary
num_patch = 50000; % number of patches to sample as the dictionary
codebook_size = 1024; % size of the dictionary

regres = 'L1'; % 'L1' or 'L2', use the sparse representation directly, or use the supports for L2 regression
% =====================================================================
% training coupled dictionaries for super-resolution

if ~skip_smp_training,
    disp('Sampling image patches...');
    [Xh, Xl] = rnd_smp_dictionary(tr_dir, patch_size, zooming, num_patch);
    save('Data/Dictionary/smp_patches.mat', 'Xh', 'Xl');
    skip_dictionary_training = false;
end;

if ~skip_dictionary_training,
    load('Data/Dictionary/smp_patches.mat');
    [Dh, Dl] = coupled_dic_train(Xh, Xl, codebook_size, lambda);
    save('Data/Dictionary/Dictionary.mat', 'Dh', 'Dl');
else
    load('Data/Dictionary/Dictionary.mat');
end;

% =====================================================================
% Process the test image

fname = 'Data/Test/1.bmp';
testIm = imread(fname); % testIm is a high resolution image, we downsample it and do super-resolution

if rem(size(testIm,1),zooming) ~=0,
    nrow = floor(size(testIm,1)/zooming)*zooming;
    testIm = testIm(1:nrow,:,:);
end;
if rem(size(testIm,2),zooming) ~=0,
    ncol = floor(size(testIm,2)/zooming)*zooming;
    testIm = testIm(:,1:ncol,:);
end;

imwrite(testIm, 'Data/Test/high.bmp', 'BMP');

lowIm = imresize(testIm,1/zooming, 'bicubic');
imwrite(lowIm,'Data/Test/low.bmp','BMP');

interpIm = imresize(lowIm,zooming,'bicubic');
imwrite(uint8(interpIm),'Data/Test/bb.bmp','BMP');

% work with the illuminance domain only
lowIm2 = rgb2ycbcr(lowIm);
lImy = double(lowIm2(:,:,1));

% bicubic interpolation for the other two channels
interpIm2 = rgb2ycbcr(interpIm);
hImcb = interpIm2(:,:,2);
hImcr = interpIm2(:,:,3);

% ======================================================================
% Super-resolution using sparse representation

disp('Start superresolution...');

[hImy] = L1SR(lImy, zooming, patch_size, overlap, Dh, Dl, lambda, regres);

ReconIm(:,:,1) = uint8(hImy);
ReconIm(:,:,2) = hImcb;
ReconIm(:,:,3) = hImcr;

nnIm = imresize(lowIm, zooming, 'nearest');
figure, imshow(nnIm);
title('Input image');
pause(1);
figure, imshow(interpIm);
title('Bicubic interpolation');
pause(1)

ReconIm = ycbcr2rgb(ReconIm);
figure,imshow(ReconIm,[]);
title('Our method');
imwrite(uint8(ReconIm),'Data/Test/L1SR.bmp','BMP');

author: thomas-koehler

repo: SupER


text:
function yn = strcontains(C1, C2)
%STRCONTAINS Does a cell array contain a particular string?
%
% YN = STRCONTAINS(C1, C2).  C1 and C2 should be a cell array of strings
% and a single string, in either order.  If the single string is one of
% the strings in the cell array, the function returns true, otherwise
% false.
%
% $Author: rcotton $  $Date: 2010/10/01 13:58:36 $	$Revision: 1.1 $
% $Copyright: Health and Safety Laboratory 2010 $

% Basic input checking
if (~iscellstr(C1) || ~ischar(C2)) && (~iscellstr(C2) || ~ischar(C1))
   warning('strcontains:BadInputs', ...
      'Inputs should be a cell array of strings and a single string');
   yn = false;
   return;
end

yn = any(strcmp(C1, C2));
end


text:
function [tests,pass,perf]=test_rdsamp(varargin)

inputs={'verbose'};
verbose=0;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

%Test the examples
test1_str=['[signal,Fs,tm]=rdsamp(''mghdb/mgh001'', [1 3 5],1000,[]);'];
test2_str=['[signal,Fs,tm]=rdsamp(''mghdb/mgh001'', [1 3 5],1000,[],2);'];
test3_str=['[signal,Fs,tm]=rdsamp(''challenge/2013/set-a/a01'',[],1000);'...
    'plot(tm,signal(:,1));close all'];
test4_str=['[sig,Fs,tm] = rdsamp(''drivedb/drive02'',[1],[],[],[],1);'];

%test3_str=['[tm2,sig2]=rdsamp(''mimic2wdb/30/3003521/3003521_0001'',[2 4 5], 30499638,30484638);'];
%test4_str=['[tm2,sig2]=rdsamp(''mimic2wdb/30/3003521/3003521_0001'',[2 4 5],[],[],2);'];
%For now avoid these test conditions from above, which still need to be
%fixed. This looks like it is an issue with reading large multi-record
%signals with N and N0 defined.

test_string={test1_str,test2_str,test3_str,test4_str};

clean_up={[''],[''],[''],['']};
[tests,pass,perf]=test_wrapper(test_string,clean_up,verbose);

author: ikarosilva

repo: wfdb-app-toolbox


text:
function Out   = RiemannExpMap(P,X)

[U Delta] = eig(P);
G = U*sqrt(Delta);
Y = inv(G)*X*inv(G)';
[V Sigma] = eig(Y);
Out = (G*V)*diag(exp(diag(Sigma)))*(G*V)';

author: alexandrebarachant

repo: covariancetoolbox


text:
function name = p26_name ( )

%*****************************************************************************80
%
%% P26_NAME returns the name of problem 26.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string NAME, the name of the problem.
%
  name = 'Patterson #1';

  return
end


text:
function [results] = trackerMain(p, im, bg_area, fg_area, area_resize_factor)
%TRACKERMAIN contains the main loop of the tracker, P contains all the parameters set in runTracker
    %% INITIALIZATION
    num_frames = numel(p.img_files);
    % used for OTB-13 benchmark
    OTB_rect_positions = zeros(num_frames, 4);
	pos = p.init_pos;
    target_sz = p.target_sz;
	num_frames = numel(p.img_files);
    % patch of the target + padding
    patch_padded = getSubwindow(im, pos, p.norm_bg_area, bg_area);
    % initialize hist model
    new_pwp_model = true;
    [bg_hist, fg_hist] = updateHistModel(new_pwp_model, patch_padded, bg_area, fg_area, target_sz, p.norm_bg_area, p.n_bins, p.grayscale_sequence);
    new_pwp_model = false;
    % Hann (cosine) window
    if isToolboxAvailable('Signal Processing Toolbox')
        hann_window = single(hann(p.cf_response_size(1)) * hann(p.cf_response_size(2))');
    else
        hann_window = single(myHann(p.cf_response_size(1)) * myHann(p.cf_response_size(2))');
    end
    % gaussian-shaped desired response, centred in (1,1)
    % bandwidth proportional to target size
    output_sigma = sqrt(prod(p.norm_target_sz)) * p.output_sigma_factor / p.hog_cell_size;
    y = gaussianResponse(p.cf_response_size, output_sigma);
    yf = fft2(y);

    %Context patches
    offset = [-target_sz(1) 0; 0 -target_sz(2); target_sz(1) 0; 0 target_sz(2)];
    elapsed_time =0;

    %% SCALE ADAPTATION INITIALIZATION
    if p.scale_adaptation
        % Code from DSST
        scale_factor = 1;
        base_target_sz = target_sz;
        scale_sigma = sqrt(p.num_scales) * p.scale_sigma_factor;
        ss = (1:p.num_scales) - ceil(p.num_scales/2);
        ys = exp(-0.5 * (ss.^2) / scale_sigma^2);
        ysf = single(fft(ys));
        if mod(p.num_scales,2) == 0
            scale_window = single(hann(p.num_scales+1));
            scale_window = scale_window(2:end);
        else
            scale_window = single(hann(p.num_scales));
        end;

        ss = 1:p.num_scales;
        scale_factors = p.scale_step.^(ceil(p.num_scales/2) - ss);

        if p.scale_model_factor^2 * prod(p.norm_target_sz) > p.scale_model_max_area
            p.scale_model_factor = sqrt(p.scale_model_max_area/prod(p.norm_target_sz));
        end

        scale_model_sz = floor(p.norm_target_sz * p.scale_model_factor);
        % find maximum and minimum scales
        min_scale_factor = p.scale_step ^ ceil(log(max(5 ./ bg_area)) / log(p.scale_step));
        max_scale_factor = p.scale_step ^ floor(log(min([size(im,1) size(im,2)] ./ target_sz)) / log(p.scale_step));
    end

    %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%
    %% MAIN LOOP
    for frame = 1:num_frames

        im = imread(fullfile(p.img_path, p.img_files{frame}));
        tic()

        if frame>1

	    %% TESTING step
            % extract patch of size bg_area and resize to norm_bg_area
            im_patch_cf = getSubwindow(im, pos, p.norm_bg_area, bg_area);
            pwp_search_area = round(p.norm_pwp_search_area / area_resize_factor);
            % extract patch of size pwp_search_area and resize to norm_pwp_search_area
            im_patch_pwp = getSubwindow(im, pos, p.norm_pwp_search_area, pwp_search_area);
            % compute feature map
            xt = getFeatureMap(im_patch_cf, p.feature_type, p.cf_response_size, p.hog_cell_size);
            % apply Hann window
            xt_windowed = bsxfun(@times, hann_window, xt);
            % compute FFT
            xtf = fft2(xt_windowed);
            % Correlation between filter and test patch gives the response
            % Solve diagonal system per pixel.
            if p.den_per_channel
% 				hf = hf_num ./ (hf_den + p.lambda);
                hf = hf_num ./ hf_den;
            else
				%hf = bsxfun(@rdivide, hf_num, sum(hf_den, 3)+p.lambda);
                hf = bsxfun(@rdivide, hf_num, sum(hf_den, 3));
            end
            %response_cf = ensure_real(ifft2(sum(conj(hf) .* xtf, 3)));
            response_cf = ensure_real(ifft2(sum(hf .* xtf, 3)));

            % Crop square search region (in feature pixels).
            response_cf = cropFilterResponse(response_cf, ...
                floor_odd(p.norm_delta_area / p.hog_cell_size));
            if p.hog_cell_size > 1
                % Scale up to match center likelihood resolution.
                response_cf = mexResize(response_cf, p.norm_delta_area,'auto');
                %response_cf = imresize(response_cf, p.norm_delta_area,'nearest');
            end

            [likelihood_map] = getColourMap(im_patch_pwp, bg_hist, fg_hist, p.n_bins, p.grayscale_sequence);
            % (TODO) in theory it should be at 0.5 (unseen colors shoud have max entropy)
            likelihood_map(isnan(likelihood_map)) = 0;

            % each pixel of response_pwp loosely represents the likelihood that
            % the target (of size norm_target_sz) is centred on it
            response_pwp = getCenterLikelihood(likelihood_map, p.norm_target_sz);

            %% ESTIMATION
            response = mergeResponses(response_cf, response_pwp, p.merge_factor, p.merge_method);
            [row, col] = find(response == max(response(:)), 1);
            center = (1+p.norm_delta_area) / 2;
            pos = pos + ([row, col] - center) / area_resize_factor;
            rect_position = [pos([2,1]) - target_sz([2,1])/2, target_sz([2,1])];

            %% SCALE SPACE SEARCH
            if p.scale_adaptation
                im_patch_scale = getScaleSubwindow(im, pos, base_target_sz, scale_factor * scale_factors, scale_window, scale_model_sz, p.hog_scale_cell_size);
                xsf = fft(im_patch_scale,[],2);
                scale_response = real(ifft(sum(sf_num .* xsf, 1) ./ (sf_den + p.lambda) ));
                recovered_scale = ind2sub(size(scale_response),find(scale_response == max(scale_response(:)), 1));
                %set the scale
                scale_factor = scale_factor * scale_factors(recovered_scale);

                if scale_factor < min_scale_factor
                    scale_factor = min_scale_factor;
                elseif scale_factor > max_scale_factor
                    scale_factor = max_scale_factor;
                end
                % use new scale to update bboxes for target, filter, bg and fg models
                target_sz = round(base_target_sz * scale_factor);
                avg_dim = sum(target_sz)/2;
                bg_area = round(target_sz + avg_dim);
                if(bg_area(2)>size(im,2)),  bg_area(2)=size(im,2)-1;    end
                if(bg_area(1)>size(im,1)),  bg_area(1)=size(im,1)-1;    end

                bg_area = bg_area - mod(bg_area - target_sz, 2);
                fg_area = round(target_sz - avg_dim * p.inner_padding);
                fg_area = fg_area + mod(bg_area - fg_area, 2);
                % Compute the rectangle with (or close to) params.fixed_area and
                % same aspect ratio as the target bboxgetScaleSubwindow
                area_resize_factor = sqrt(p.fixed_area/prod(bg_area));
            end

            if p.visualization_dbg==1
                mySubplot(2,1,5,1,im_patch_cf,'FG+BG','gray');
                mySubplot(2,1,5,2,likelihood_map,'obj.likelihood','parula');
                mySubplot(2,1,5,3,response_cf,'CF response','parula');
                mySubplot(2,1,5,4,response_pwp,'center likelihood','parula');
                mySubplot(2,1,5,5,response,'merged response','parula');
                drawnow
            end
        end

        %% TRAINING
        % extract patch of size bg_area and resize to norm_bg_area
        im_patch_bg = getSubwindow(im, pos, p.norm_bg_area, bg_area);
        % compute feature map, of cf_response_size
        xt = getFeatureMap(im_patch_bg, p.feature_type, p.cf_response_size, p.hog_cell_size);
        % apply Hann window
        xt = bsxfun(@times, hann_window, xt);
        % compute FFT
        xtf = fft2(xt);
        kf = (conj(xtf) .* xtf); %/ prod(p.cf_response_size);

        kfn = zeros([size(xtf) length(offset)]);
        for j=1:length(offset)
        %obtain a subwindow close to target for regression to 0
		im_patch_bgn = getSubwindow(im, pos+offset(j,:), p.norm_bg_area, bg_area);
        % compute feature map, of cf_response_size
        xtn = getFeatureMap(im_patch_bgn, p.feature_type, p.cf_response_size, p.hog_cell_size);
        % apply Hann window
        xtn = bsxfun(@times, hann_window, xtn);
        % compute FFT
        xtfn = fft2(xtn);
        kfn(:,:,:,j) = conj(xtfn) .*xtfn; %/ prod(p.cf_response_size);
        end

        %% FILTER UPDATE
        % Compute expectations over circular shifts,
        % therefore divide by number of pixels.
		new_hf_num = bsxfun(@times, conj(xtf),yf); %/ prod(p.cf_response_size);
		new_hf_den = kf + p.lambda + p.lambda2.*sum(kfn,4);

        if frame == 1
            % first frame, train with a single image
		    hf_den = new_hf_den;
		    hf_num = new_hf_num;
		else
		    % subsequent frames, update the model by linear interpolation
        	hf_den = (1 - p.learning_rate_cf) * hf_den + p.learning_rate_cf * new_hf_den;
	   	 	hf_num = (1 - p.learning_rate_cf) * hf_num + p.learning_rate_cf * new_hf_num;

            %% BG/FG MODEL UPDATE
            % patch of the target + padding
            [bg_hist, fg_hist] = updateHistModel(new_pwp_model, im_patch_bg, bg_area, fg_area, target_sz, p.norm_bg_area, p.n_bins, p.grayscale_sequence, bg_hist, fg_hist, p.learning_rate_pwp);
        end

        %% SCALE UPDATE
        if p.scale_adaptation
            im_patch_scale = getScaleSubwindow(im, pos, base_target_sz, scale_factor*scale_factors, scale_window, scale_model_sz, p.hog_scale_cell_size);
            xsf = fft(im_patch_scale,[],2);
            new_sf_num = bsxfun(@times, ysf, conj(xsf));
            new_sf_den = sum(xsf .* conj(xsf), 1);
            if frame == 1,
                sf_den = new_sf_den;
                sf_num = new_sf_num;
            else
                sf_den = (1 - p.learning_rate_scale) * sf_den + p.learning_rate_scale * new_sf_den;
                sf_num = (1 - p.learning_rate_scale) * sf_num + p.learning_rate_scale * new_sf_num;
            end
        end

        % update bbox position
        if frame==1, rect_position = [pos([2,1]) - target_sz([2,1])/2, target_sz([2,1])]; end

        rect_position_padded = [pos([2,1]) - bg_area([2,1])/2, bg_area([2,1])];

        OTB_rect_positions(frame,:) = rect_position;

        if p.fout > 0,  fprintf(p.fout,'%.2f,%.2f,%.2f,%.2f\n', rect_position(1),rect_position(2),rect_position(3),rect_position(4));   end

        elapsed_time = elapsed_time + toc();

        %% VISUALIZATION
        if p.visualization == 1
            if isToolboxAvailable('Computer Vision System Toolbox')
                im = insertShape(im, 'Rectangle', rect_position, 'LineWidth', 4, 'Color', 'black');
                im = insertShape(im, 'Rectangle', rect_position_padded, 'LineWidth', 4, 'Color', 'yellow');
                % Display the annotated video frame using the video player object.
                step(p.videoPlayer, im);
            else
                figure(1)
                imshow(im)
                rectangle('Position',rect_position, 'LineWidth',2, 'EdgeColor','g');
                rectangle('Position',rect_position_padded, 'LineWidth',2, 'LineStyle','--', 'EdgeColor','b');
                drawnow
            end
        end
    end

    % save data for OTB-13 benchmark
    results.type = 'rect';
    results.res = OTB_rect_positions;
    results.fps = num_frames/(elapsed_time );
end

% Reimplementation of Hann window (in case signal processing toolbox is missing)
function H = myHann(X)
    H = .5*(1 - cos(2*pi*(0:X-1)'/(X-1)));
end

% We want odd regions so that the central pixel can be exact
function y = floor_odd(x)
    y = 2*floor((x-1) / 2) + 1;
end

function y = ensure_real(x)
    assert(norm(imag(x(:))) <= 1e-5 * norm(real(x(:))));
    y = real(x);
end

repo: Context-Aware-CF-Tracking


text: any(cellfun(@exist,{'lsl_append_child', 'lsl_append_child_value', 'lsl_append_copy', 'lsl_child', 'lsl_child_value', 'lsl_child_value_n', 'lsl_close_stream', 'lsl_create_inlet', 'lsl_create_outlet', 'lsl_create_streaminfo', 'lsl_destroy_inlet', 'lsl_destroy_outlet', 'lsl_destroy_streaminfo', 'lsl_empty', 'lsl_first_child', 'lsl_freelib_', 'lsl_get_channel_count', 'lsl_get_channel_format', 'lsl_get_created_at', 'lsl_get_desc', 'lsl_get_fullinfo', 'lsl_get_hostname', 'lsl_get_info', 'lsl_get_name', 'lsl_get_nominal_srate', 'lsl_get_session_id', 'lsl_get_source_id', 'lsl_get_type', 'lsl_get_uid', 'lsl_get_version', 'lsl_get_xml', 'lsl_have_consumers', 'lsl_is_text', 'lsl_last_child', 'lsl_loadlib_', 'lsl_local_clock', 'lsl_name', 'lsl_next_sibling', 'lsl_next_sibling_n', 'lsl_open_stream', 'lsl_parent', 'lsl_prepend_child', 'lsl_prepend_child_value', 'lsl_prepend_copy', 'lsl_previous_sibling', 'lsl_previous_sibling_n', 'lsl_protocol_version', 'lsl_pull_chunk_d', 'lsl_pull_sample_d', 'lsl_pull_sample_str', 'lsl_push_chunk', 'lsl_push_sample', 'lsl_remove_child', 'lsl_remove_child_n', 'lsl_resolve_all_', 'lsl_resolve_bypred_', 'lsl_resolve_byprop_', 'lsl_set_child_value', 'lsl_set_name', 'lsl_set_value', 'lsl_time_correction', 'lsl_value', 'lsl_wait_for_consumers'})==3)


text:
function setcurrentfile(filename, lboxhandle, currfilehandle)

set(lboxhandle, 'UserData', filename );

[path,name,ext] = fileparts(filename);

set(currfilehandle, 'String', name);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
classdef GemCheckbox < GemUserInterfaceObject
    % GemCheckbox GEM class for a check box control
    %
    %
    %     Licence
    %     -------
    %     Part of GEM. https://github.com/tomdoel/gem
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        FontSize
        FontColour
        BackgroundColour
        HorizontalAlignment
    end

    properties (Access = private)
        Text
        Tag
        ToolTip
        Checked
    end

    events
        CheckChanged
    end

    methods
        function obj = GemCheckbox(parent, text, tooltip, tag)
            obj = obj@GemUserInterfaceObject(parent);
            obj.Text = text;
            obj.ToolTip = tooltip;
            obj.Tag = tag;
            obj.FontSize = 11;
            obj.HorizontalAlignment = 'left';
            obj.Checked = false;
            obj.FontColour = obj.StyleSheet.TextPrimaryColour;
            obj.BackgroundColour = obj.StyleSheet.BackgroundColour;
        end

        function CreateGuiComponent(obj, position)
            obj.GraphicalComponentHandle = uicontrol('Style', 'checkbox', 'Parent', obj.Parent.GetContainerHandle, ...
                'Units', 'pixels', 'Position', position, 'Tag', obj.Tag, 'ToolTipString', obj.ToolTip, ...
                'ForegroundColor', obj.FontColour, 'BackgroundColor', obj.BackgroundColour, ...
                'FontUnits', 'pixels', 'FontSize', obj.FontSize, 'FontAngle', 'normal', ...
                'HorizontalAlignment', obj.HorizontalAlignment, 'String', obj.Text, ...
                'Callback', @obj.CheckboxCallback, 'Value', obj.Checked);
        end

        function ChangeChecked(obj, checked)
            obj.Checked = checked;
            if obj.ComponentHasBeenCreated && ishandle(obj.GraphicalComponentHandle)
                set(obj.GraphicalComponentHandle, 'Value', checked);
            end
        end
    end

    methods (Access = protected)
        function CheckboxCallback(obj, hObject, ~, ~)
            checked = get(hObject,'Value');
            obj.Checked = checked;
            notify(obj, 'CheckChanged', CoreEventData(checked));
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
classdef TotalStateOfChargeOptimizationVariable < AbstractOptimizationVariable
    %TotalStateOfChargeOptimizationVariable Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) TotalVehicleEpsStateOfChargeTermCondition = TotalVehicleEpsStateOfChargeTermCondition(0);

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        useTf(1,1) = false;
    end

    methods
        function obj = TotalStateOfChargeOptimizationVariable(varObj)
            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = [];

            if(obj.useTf)
                x = obj.varObj.perialt;
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            lb = obj.lb(obj.useTf);
            ub = obj.ub(obj.useTf);
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lwrBnd;
            ub = obj.uprBnd;
        end

        function setBndsForVariable(obj, lb, ub)
            obj.lb = lb;
            obj.ub = ub;
        end

        function useTf = getUseTfForVariable(obj)
            useTf = obj.useTf;
        end

        function setUseTfForVariable(obj, useTf)
            obj.useTf = useTf;
        end

        function updateObjWithVarValue(obj, x)
            obj.varObj.totalStateOfCharge = x;
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum)
            nameStrs = {sprintf('Event %i Power Total State of Charge Termination Condition', evtNum)};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout=minDist(varargin)

% function [D1,minIND]=minDist(V1,V2,maxVarSize,selfAvoid,numFreeBytes)

%% Parse input
if nargin<2
    error('Insufficient input arguments');
end

V1=varargin{1};
V2=varargin{2};
switch nargin
    case 2
        maxVarSize=[]; %Empty will force calcucation below
        selfAvoid=0;
        numFreeBytes=[];
    case 3
        maxVarSize=varargin{3};
        selfAvoid=0;
        numFreeBytes=[];
    case 4
        maxVarSize=varargin{3};
        selfAvoid=varargin{4};
        numFreeBytes=[];
    case 5
        maxVarSize=varargin{3};
        selfAvoid=varargin{4};
        numFreeBytes=varargin{5};
end

%Get free memory
if isempty(numFreeBytes)
    [numFreeBytes]=freeMemory;
end

%Get max variable size available
if isempty(maxVarSize)
    maxVarSize=numFreeBytes/2;
end

if isnan(maxVarSize)
    numSteps=1;
else
    %Derive class dependent variable size
    [~,b1]=maxnumel(V1(1),numFreeBytes);
    [~,b2]=maxnumel(V2(1),numFreeBytes);
    b=max([b1 b2]);
    numelVar=numel(V1)*numel(V2);
    varSize=numelVar*b;

    numSteps=ceil(varSize/maxVarSize);
    indSteps=round(linspace(0,size(V1,1),numSteps));
    indSteps=sort(unique(indSteps));
    numSteps=numel(indSteps);
end

if numSteps>1 %In steps
    D1=zeros(size(V1,1),1);
    minIND=zeros(size(V1,1),1);
    for q=1:1:numSteps-1
        v1=V1(indSteps(q)+1:indSteps(q+1),:);
        try
            d=dist(v1,V2'); %dist from Neural network toolbox
        catch
            d=distND(v1,V2); %GIBBON's dist function
        end
        if selfAvoid
            %Set "diagonal" to something too large so self is avoided in
            %minimum (could use NaN and nanmin but the latter is a toolbox
            %function)
            I=1:size(v1,1);
            J=indSteps(q)+1:indSteps(q+1);
            ind=sub2ind(size(d),I,J); %Indices of selfies
            d(ind)=1+max(d(:)); %Overide selfies
        end

        [min_d,min_ind]=min(d,[],2);
        D1(indSteps(q)+1:indSteps(q+1))=min_d;
        minIND(indSteps(q)+1:indSteps(q+1))=min_ind;
    end
else %In one go
    try
        D=dist(V1,V2'); %dist from Neural network toolbox
    catch
        D=distND(V1,V2); %GIBBON's dist function
    end
    if selfAvoid
        %Set "diagonal" to something too large so self is avoided in
        %minimum (could use NaN and nanmin but the latter is a toolbox
        %function)
        L=eye(size(D))>0;
        D(L)=1+max(D(:));
    end
    [D1,minIND]=min(D,[],2);
    D1=D1(:);
    minIND=minIND(:);
end

switch nargout
    case 1
        varargout{1}=D1;
    case 2
        varargout{1}=D1;
        varargout{2}=minIND;
    otherwise
        error('wrong number of output arguments');
end

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

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC