text:
% create metabolite database

metabolite_structure = struct();
if 1
    % this file seems to contain a lot of errors
    % [NUM,TXT,RAW]=xlsread('VMH_Metabolites.xlsx');

    % use instead as starting point as it comes from rBioNet
    [NUM,TXT,RAW]=xlsread('MetabolitesAGORA2_1.xlsx');
    for i = 2 : size(RAW,1)
        Ori = RAW{i,1};
        RAW{i,1} = regexprep(RAW{i,1},'-','_minus_');
        RAW{i,1} = regexprep(RAW{i,1},'(','_parentO_');
        RAW{i,1} = regexprep(RAW{i,1},')','_parentC_');
        metabolite_structure.(strcat('VMH_',RAW{i,1})) = struct();

        for j = 2:size(RAW,2)
            if j ~= 3
                metabolite_structure.(strcat('VMH_',RAW{i,1})).(RAW{1,j}) = RAW{i,j};
                metabolite_structure.(strcat('VMH_',RAW{i,1})).VMHId = Ori;
            end
        end
    end
    [metabolite_structure] = cleanUpMetabolite_structure(metabolite_structure)
    removeErrors = 1;
    [metabolite_structure,errorFlag] = sanityCheckMetIds(metabolite_structure,removeErrors);
    metabolite_structureStart = metabolite_structure;

end
if 1
    % add fields
    field2Add={
        'phenolExplorer'
        'neutralFormula'
        'chembl'
        'miriam'
        'hepatonetId'
        'iupac'
        'echa_id'
        'fda_id'
        'iuphar_id'
        'mesh_id'
        'chodb_id'
        'gtopdb'
        'inchiKey'
        'avgmolweight'
        'monoisotopicweight'
        'keggId'
        'inchiString'
        'pubChemId'
        'cheBIId'
        'hmdb'
        'pdmapName'
        'reconMap'
        'reconMap3'
        'food_db'
        'chemspider'
        'biocyc'
        'biggId'
        'wikipedia'
        'drugbank'
        'seed'
        'metanetx'
        'knapsack'
        'metlin'
        'casRegistry'
        'epa_id'
        'inchiKey'
        'smile'
        'lipidmaps'
        'reactome'
        'GNPS'
        'Recon3D'
        'Agora2'
        'massbank'
        'MoNa'
        'bindingdb'
        'metabolights'
        'rhea'
        'swisslipids'};
  %  field2Add = unique(field2Add);
    Mets = fieldnames(metabolite_structure);
    fields = fieldnames(metabolite_structure.(Mets{1}));

    [missingfields,map] = setdiff(field2Add,fields);
    for i = 1 : length(Mets)
        for j = 1 : length(missingfields)
            metabolite_structure.(Mets{i}).(missingfields{j}) = NaN;
        end
    end
end
cd FileDumps/
%[IDs,IDcountA,Table] = getStats(metabolite_structure);
%[metabolite_structure] = getMetIdsFromUniChemOffline(metabolite_structure);
assignMetaboliteIDs;
[IDs,IDcountB,Table] = getStats(metabolite_structure);
[metabolite_structure] = getMetIdsFromUniChemOffline(metabolite_structure);
[metabolite_structure,IDsAdded ] = parseChebiIdWebPage(metabolite_structure);
save vmh_met_anno
%assignMetaboliteIDs;
%[IDs,IDcountC,Table] = getStats(metabolite_structure);
%save vmh_met_anno
inchiKeyCheck = 1;
inchiStringCheck = 1;
inchiKeyAltCheck = 0;
cheBIIdCheck = 0;
drugBankCheck = 0;
pubChemIdCheck = 0;
keggIdCheck = 0;
hmdbCheck = 0;
vmhIdCheck = 0;
[metabolite_structure,IDsAdded] = getMetIdsFromUniChem(metabolite_structure,vmhIdCheck,cheBIIdCheck,drugBankCheck, pubChemIdCheck,keggIdCheck,hmdbCheck, inchiKeyCheck, inchiStringCheck,inchiKeyAltCheck );
[IDs,IDcountD,Table] = getStats(metabolite_structure);
save vmh_met_anno
translation = {'keggId' 'kegg'
    'inchiKey'  'inchikey'
    'hmdb'  'Human%20Metabolome%20Database'
    'cheBIId'   'chebi'
    'pubChemId' 'PubChem%20CID'
    'inchiString'   'inchi%20code'
    'bindingdb' 'bindingdb'
    'drugbank'  'drugbank'
    'biocyc'    'biocyc'
    'chemspider'    'chemspider'
    'casRegistry'   'cas'
    'lipidmaps' 'lipidmaps'
    'epa_id'    'epa%20dsstox'
    };
for i = 1 : 5%size(translation,1)
    for j = 1:size(translation,1)
        if i~=j
            i
            j
            [metabolite_structure,IDsAdded] = getIDsfromFiehnLab(metabolite_structure, translation{i,1},translation{j,1});
            save vmh_met_anno
        end
    end
    done = i
end
save vmh_met_anno
inchiKeyCheck = 1;
inchiStringCheck = 1;
inchiKeyAltCheck = 1;
[metabolite_structure,IDsAdded] = getMetIdsFromInchiKeys(metabolite_structure, inchiKeyCheck, inchiStringCheck,inchiKeyAltCheck);
save vmh_met_anno
assignMetaboliteIDs;
[IDs,IDcountD,Table] = getStats(metabolite_structure);
save vmh_met_anno

writecell(Table);

author: opencobra

repo: cobratoolbox


text:
function x = daub12_transform_inverse ( n, y )

%*****************************************************************************80
%
%% DAUB12_TRANSFORM_INVERSE inverts the DAUB12 transform of a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%    N must be a power of 2 and at least 4.
%
%    Input, real Y(N), the transformed vector.
%
%    Output, real X(N), the original vector.
%
  c = [ ...
    0.1115407433501095; ...
    0.4946238903984533; ...
    0.7511339080210959; ...
    0.3152503517091982; ...
   -0.2262646939654400; ...
   -0.1297668675672625; ...
    0.0975016055873225; ...
    0.0275228655303053; ...
   -0.0315820393174862; ...
    0.0005538422011614; ...
    0.0047772575109455; ...
   -0.0010773010853085 ];
  p = 11;
  x(1:n,1) = y(1:n);
  m = 4;
  q = floor ( ( p - 1 ) / 2 );

  while ( m <= n )

    z(1:m,1) = 0.0;

    j = 1;

    mh = floor ( m / 2 );

    for i = - q + 1 : mh - q

      for k = 0 : 2 : p - 1
        i0 = i4_wrap ( i      + k / 2,     1,      mh );
        i1 = i4_wrap ( i + mh + k / 2,     mh + 1, m  );
        z(j,1)   = z(j,1)   + c(p-k)   * x(i0) + c(k+2) * x(i1);
        z(j+1,1) = z(j+1,1) + c(p-k+1) * x(i0) - c(k+1) * x(i1);
      end

      j = j + 2;

    end

    x(1:m,1) = z(1:m);

    m = m * 2;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
echo off
%load seed.mat
rand('seed',0);
echo on
% This script shows how to use the ga. You should see the demos for
% more information as well. gademo1, gademo2, gademo3
global bounds

% Crossover Operators
xFns = 'arithXover heuristicXover simpleXover';
xOpts = [2 0; 2 3; 2 0];

% Mutation Operators
mFns = 'boundaryMutation multiNonUnifMutation nonUnifMutation unifMutation';
mOpts = [4 0 0;6 10 3;4 10 3;4 0 0];

% Termination Operators
termFns = 'maxGenTerm';
termOps = [10];

% Selection Function
selectFn = 'normGeomSelect';
selectOps = [0.06];

% Evaluation Function takes two options
% prob to use gradient, prob to perform Lamarkian evolution
evalFn = 'gaZBGradEval';
evalOps = [1.00 1.00];

% Bounds on the variables
bounds = [-3 12.1; 4.1 5.8];

% GA Options [epsilon float/binar display]
gaOpts=[1e-6 1 1];

% Generate an intialize population of size 80
startPop = initializega(80,bounds,evalFn,evalOps,[1e-6 1]);
evalOps = [1.00 0.00]; % 1 - Peform learning 0-Do not update

% Lets run the GA using Baldwinian Evolution

[x endPop bestPop trace]=ga(bounds,evalFn,evalOps,startPop,gaOpts,...
    termFns,termOps,selectFn,selectOps,xFns,xOpts,mFns,mOpts);

% x is the best solution found
x
pause

% endPop is the ending population
endPop
pause

% bestPop is the best solution tracked over generations
bestPop
pause

% trace is a trace of the best value and average value of generations
trace
pause

% Plot the best over time
clf
plot(trace(:,1),trace(:,2));
pause

% Add the average to the graph
hold on
plot(trace(:,1),trace(:,3));
pause

% Lets run the GA using Lamarkian Evolution
evalOps = [1.00 1.00];

[x endPop bestPop trace]=ga(bounds,evalFn,evalOps,startPop,gaOpts,...
    termFns,termOps,selectFn,selectOps,xFns,xOpts,mFns,mOpts);

% x is the best solution found
x
pause

% endPop is the ending population
endPop
pause

% bestPop is the best solution tracked over generations
bestPop
pause

% trace is a trace of the best value and average value of generations
trace
pause

% Plot the best over time
clf
plot(trace(:,1),trace(:,2));
pause

% Add the average to the graph
hold on
plot(trace(:,1),trace(:,3));
pause

author: HuangCongQing


text:
function Y = squeeze(X)
%SQUEEZE Remove singleton dimensions from a tensor.
%
%   Y = SQUEEZE(X) returns a tensor Y with the same elements as
%   X but with all the singleton dimensions removed.  A singleton
%   is a dimension such that size(X,dim)==1.
%
%   If X has *only* singleton dimensions, then Y is a scalar.
%
%   Examples
%   squeeze( tenrand([2,1,3]) ) %<-- returns a 2-by-3 tensor
%   squeeze( tenrand([1 1]) ) %<-- returns a scalar
%
%   See also TENSOR, TENSOR/RESHAPE, TENSOR/PERMUTE.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if all(X.size > 1)
  % No singleton dimensions to squeeze
  Y = X;
else
  idx = find(X.size > 1);
  if numel(idx) == 0
    % Scalar case - only singleton dimensions
    Y = X.data;
  else
    siz = X.size(idx);
    Y = tensor(squeeze(X.data),siz);
  end
end

return;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
classdef Violin < handle
    % Violin creates violin plots for some data
    %   A violin plot is an easy to read substitute for a box plot
    %   that replaces the box shape with a kernel density estimate of
    %   the data, and optionally overlays the data points itself.
    %   It is also possible to provide two sets of data which are supposed
    %   to be compared by plotting each column of the two datasets together
    %   on each side of the violin.
    %
    %   Additional constructor parameters include the width of the
    %   plot, the bandwidth of the kernel density estimation, the
    %   X-axis position of the violin plot, and the categories.
    %
    %   Use violinplot for a
    %   boxplot-like wrapper for
    %   interactive plotting.
    %
    %   See for more information on Violin Plots:
    %   J. L. Hintze and R. D. Nelson, "Violin plots: a box
    %   plot-density trace synergism," The American Statistician, vol.
    %   52, no. 2, pp. 181-184, 1998.
    %
    % Violin Properties:
    %    ViolinColor    - Fill color of the violin area and data points.
    %                     Can be either a matrix nx3 or an array of up to two
    %                     cells containing nx3 matrices.
    %                     Defaults to the next default color cycle.
    %    ViolinAlpha    - Transparency of the violin area and data points.
    %                     Can be either a single scalar value or an array of
    %                     up to two cells containing scalar values.
    %                     Defaults to 0.3.
    %    EdgeColor      - Color of the violin area outline.
    %                     Defaults to [0.5 0.5 0.5]
    %    BoxColor       - Color of the box, whiskers, and the outlines of
    %                     the median point and the notch indicators.
    %                     Defaults to [0.5 0.5 0.5]
    %    MedianColor    - Fill color of the median and notch indicators.
    %                     Defaults to [1 1 1]
    %    ShowData       - Whether to show data points.
    %                     Defaults to true
    %    ShowNotches    - Whether to show notch indicators.
    %                     Defaults to false
    %    ShowMean       - Whether to show mean indicator.
    %                     Defaults to false
    %    ShowBox        - Whether to show the box.
    %                     Defaults to true
    %    ShowMedian     - Whether to show the median indicator.
    %                     Defaults to true
    %    ShowWhiskers   - Whether to show the whiskers
    %                     Defaults to true
    %    HalfViolin     - Whether to do a half violin(left, right side) or
    %                     full. Defaults to full.
    %    QuartileStyle - Option on how to display quartiles, with a
    %                     boxplot, shadow or none. Defaults to boxplot.
    %    DataStyle      - Defines the style to show the data points. Opts:
    %                     'scatter', 'histogram' or 'none'. Default is 'scatter'.
    %
    %
    % Violin Children:
    %    ScatterPlot    - scatter plot of the data points
    %    ScatterPlot2   - scatter second plot of the data points
    %    ViolinPlot     - fill plot of the kernel density estimate
    %    ViolinPlot2    - fill second plot of the kernel density estimate
    %    BoxPlot        - fill plot of the box between the quartiles
    %    WhiskerPlot    - line plot between the whisker ends
    %    MedianPlot     - scatter plot of the median (one point)
    %    NotchPlots     - scatter plots for the notch indicators
    %    MeanPlot       - line plot at mean value

    % Copyright (c) 2016, Bastian Bechtold
    % This code is released under the terms of the BSD 3-clause license

    properties (Access=public)
        ScatterPlot     % scatter plot of the data points
        ScatterPlot2    % comparison scatter plot of the data points
        ViolinPlot      % fill plot of the kernel density estimate
        ViolinPlot2     % comparison fill plot of the kernel density estimate
        BoxPlot         % fill plot of the box between the quartiles
        WhiskerPlot     % line plot between the whisker ends
        MedianPlot      % scatter plot of the median (one point)
        NotchPlots      % scatter plots for the notch indicators
        MeanPlot        % line plot of the mean (horizontal line)
        HistogramPlot   % histogram of the data
        ViolinPlotQ     % fill plot of the Quartiles as shadow
    end

    properties (Dependent=true)
        ViolinColor         % fill color of the violin area and data points
        ViolinAlpha         % transparency of the violin area and data points
        MarkerSize          % marker size for the data dots
        MedianMarkerSize    % marker size for the median dot
        LineWidth           % linewidth of the median plot
        EdgeColor           % color of the violin area outline
        BoxColor            % color of box, whiskers, and median/notch edges
        BoxWidth            % width of box between the quartiles in axis space (default 10% of Violin plot width, 0.03)
        MedianColor         % fill color of median and notches
        ShowData            % whether to show data points
        ShowNotches         % whether to show notch indicators
        ShowMean            % whether to show mean indicator
        ShowBox             % whether to show the box
        ShowMedian          % whether to show the median line
        ShowWhiskers        % whether to show the whiskers
        HalfViolin          % whether to do a half violin(left, right side) or full
    end

    methods
        function obj = Violin(data, pos, varargin)
            %Violin plots a violin plot of some data at pos
            %   VIOLIN(DATA, POS) plots a violin at x-position POS for
            %   a vector of DATA points.
            %
            %   VIOLIN(..., 'PARAM1', val1, 'PARAM2', val2, ...)
            %   specifies optional name/value pairs:
            %     'Width'        Width of the violin in axis space.
            %                    Defaults to 0.3
            %     'Bandwidth'    Bandwidth of the kernel density
            %                    estimate. Should be between 10% and
            %                    40% of the data range.
            %     'ViolinColor'  Fill color of the violin area
            %                    and data points.Can be either a matrix
            %                    nx3 or an array of up to two cells
            %                    containing nx3 matrices.
            %     'ViolinAlpha'  Transparency of the violin area and data
            %                    points. Can be either a single scalar
            %                    value or an array of up to two cells
            %                    containing scalar values. Defaults to 0.3.
            %     'MarkerSize'   Size of the data points, if shown.
            %                    Defaults to 24
            % 'MedianMarkerSize' Size of the median indicator, if shown.
            %                    Defaults to 36
            %     'EdgeColor'    Color of the violin area outline.
            %                    Defaults to [0.5 0.5 0.5]
            %     'BoxColor'     Color of the box, whiskers, and the
            %                    outlines of the median point and the
            %                    notch indicators. Defaults to
            %                    [0.5 0.5 0.5]
            %     'MedianColor'  Fill color of the median and notch
            %                    indicators. Defaults to [1 1 1]
            %     'ShowData'     Whether to show data points.
            %                    Defaults to true
            %     'ShowNotches'  Whether to show notch indicators.
            %                    Defaults to false
            %     'ShowMean'     Whether to show mean indicator.
            %                    Defaults to false
            %     'ShowBox'      Whether to show the box
            %                    Defaults to true
            %     'ShowMedian'   Whether to show the median line
            %                    Defaults to true
            %     'ShowWhiskers' Whether to show the whiskers
            %                    Defaults to true
            %     'HalfViolin'   Whether to do a half violin(left, right side) or
            %                    full. Defaults to full.
            %   'QuartileStyle'  Option on how to display quartiles, with a
            %                    boxplot or as a shadow. Defaults to boxplot.
            %     'DataStyle'    Defines the style to show the data points. Opts:
            %                   'scatter', 'histogram' or 'none'. Default is 'Scatter'.

            st = dbstack; % get the calling function for reporting errors
            namefun = st.name;
            args = obj.checkInputs(data, pos, varargin{:});

            if length(data)==1
                data2 = [];
                data = data{1};

            else
                data2 = data{2};
                data = data{1};
            end

            if isempty(args.ViolinColor)
                Release= strsplit(version('-release'), {'a','b'}); %Check release
                if str2num(Release{1})> 2019 || strcmp(version('-release'), '2019b')
                     C = colororder;
                else
                     C = lines;
                end

                if pos > length(C)
                    C = lines;
                end
                args.ViolinColor = {repmat(C,ceil(size(data,2)/length(C)),1)};
            end

            data = data(not(isnan(data)));
            data2 = data2(not(isnan(data2)));
            if numel(data) == 1
                obj.MedianPlot = scatter(pos, data, 'filled');
                obj.MedianColor = args.MedianColor;
                obj.MedianPlot.MarkerEdgeColor = args.EdgeColor;
                return
            end

            hold('on');

            %% Calculate kernel density estimation for the violin
            [density, value, width] = obj.calcKernelDensity(data, args.Bandwidth, args.Width);

            % also calculate the kernel density of the comparison data if
            % provided
            if ~isempty(data2)
                [densityC, valueC, widthC] = obj.calcKernelDensity(data2, args.Bandwidth, args.Width);
            end

            %% Plot the data points within the violin area
            if length(density) > 1
                [~, unique_idx] = unique(value);
                jitterstrength = interp1(value(unique_idx), density(unique_idx)*width, data, 'linear','extrap');
            else % all data is identical:
                jitterstrength = density*width;
            end
            if isempty(data2) % if no comparison data
                jitter = 2*(rand(size(data))-0.5); % both sides
            else
                jitter = rand(size(data)); % only right side
            end
            switch args.HalfViolin % this is more modular
                case 'left'
                    jitter = -1*(rand(size(data))); %left
                case 'right'
                    jitter = 1*(rand(size(data))); %right
                case 'full'
                    jitter = 2*(rand(size(data))-0.5);
            end
            % Make scatter plot
            switch args.DataStyle
                case 'scatter'
                    if ~isempty(data2)
                        jitter = 1*(rand(size(data))); %right
                        obj.ScatterPlot = ...
                            scatter(pos + jitter.*jitterstrength, data, args.MarkerSize, 'filled');
                        % plot the data points within the violin area
                        if length(densityC) > 1
                            jitterstrength = interp1(valueC, densityC*widthC, data2);
                        else % all data is identical:
                            jitterstrength = densityC*widthC;
                        end
                        jitter = -1*rand(size(data2));% left
                        obj.ScatterPlot2 = ...
                            scatter(pos + jitter.*jitterstrength, data2, args.MarkerSize, 'filled');
                    else
                        obj.ScatterPlot = ...
                            scatter(pos + jitter.*jitterstrength, data, args.MarkerSize, 'filled');

                    end
                case 'histogram'
                    [counts,edges] = histcounts(data, size(unique(data),1));
                    switch args.HalfViolin
                        case 'right'
                            obj.HistogramPlot= plot([pos-((counts')/max(counts))*max(jitterstrength)*2, pos*ones(size(counts,2),1)]',...
                                [edges(1:end-1)+max(diff(edges))/2; edges(1:end-1)+max(diff(edges))/2],'-','LineWidth',1, 'Color', 'k');
                        case 'left'
                            obj.HistogramPlot= plot([pos*ones(size(counts,2),1), pos+((counts')/max(counts))*max(jitterstrength)*2]',...
                                [edges(1:end-1)+max(diff(edges))/2; edges(1:end-1)+max(diff(edges))/2],'-','LineWidth',1, 'Color', 'k');
                        otherwise
                            fprintf([namefun, ' No histogram/bar plot option available for full violins, as it would look overcrowded.\n'])
                    end
                case 'none'
            end

            %% Plot the violin
            halfViol= ones(1, size(density,2));
            if isempty(data2) % if no comparison data
                switch args.HalfViolin
                    case 'right'
                        obj.ViolinPlot =  ... % plot color will be overwritten later
                            fill([pos+density*width halfViol*pos], ...
                            [value value(end:-1:1)], [1 1 1]);
                    case 'left'
                        obj.ViolinPlot =  ... % plot color will be overwritten later
                            fill([halfViol*pos pos-density(end:-1:1)*width], ...
                            [value value(end:-1:1)], [1 1 1]);
                    case 'full'
                        obj.ViolinPlot =  ... % plot color will be overwritten later
                            fill([pos+density*width pos-density(end:-1:1)*width], ...
                            [value value(end:-1:1)], [1 1 1]);
                end
            else
                % plot right half of the violin
                obj.ViolinPlot =  ...
                    fill([pos+density*width pos-density(1)*width], ...
                    [value value(1)], [1 1 1]);
                % plot left half of the violin
                obj.ViolinPlot2 =  ...
                    fill([pos-densityC(end)*widthC pos-densityC(end:-1:1)*widthC], ...
                    [valueC(end) valueC(end:-1:1)], [1 1 1]);
            end

            %% Plot the quartiles within the violin
            quartiles = quantile(data, [0.25, 0.5, 0.75]);
            flat= [halfViol*pos halfViol*pos];
            switch args.QuartileStyle
                case 'shadow'
                    switch args.HalfViolin
                        case 'right'
                            w = [pos+density*width halfViol*pos];
                            h= [value value(end:-1:1)];
                        case 'left'
                            w = [halfViol*pos pos-density(end:-1:1)*width];
                            h= [value value(end:-1:1)];
                        case 'full'
                            w = [pos+density*width pos-density(end:-1:1)*width];
                            h= [value value(end:-1:1)];
                    end
                    w(h<quartiles(1))=flat(h<quartiles(1));
                    w(h>quartiles(3))=flat((h>quartiles(3)));
                    obj.ViolinPlotQ =  ... % plot color will be overwritten later
                        fill(w, ...
                        h, [1 1 1]);
                case 'boxplot'
                    obj.BoxPlot = ... % plot color will be overwritten later
                        fill(pos+[-1,1,1,-1]*args.BoxWidth, ...
                        [quartiles(1) quartiles(1) quartiles(3) quartiles(3)], ...
                        [1 1 1]);
                case 'none'
            end

            %% Plot the data mean
            meanValue = mean(data);
            if length(density) > 1
                [~, unique_idx] = unique(value);
                meanDensityWidth = interp1(value(unique_idx), density(unique_idx), meanValue, 'linear','extrap')*width;
            else % all data is identical:
                meanDensityWidth = density*width;
            end
            if meanDensityWidth<args.BoxWidth/2
                meanDensityWidth=args.BoxWidth/2;
            end
            switch args.HalfViolin
                case 'right'
                    obj.MeanPlot = plot(pos+[0,1].*meanDensityWidth, ...
                        [meanValue, meanValue]);
                case 'left'
                    obj.MeanPlot = plot(pos+[-1,0].*meanDensityWidth, ...
                        [meanValue, meanValue]);
                case 'full'
                    obj.MeanPlot = plot(pos+[-1,1].*meanDensityWidth, ...
                        [meanValue, meanValue]);
            end
            obj.MeanPlot.LineWidth = 1;

            %% Plot the median, notch, and whiskers
            IQR = quartiles(3) - quartiles(1);
            lowhisker = quartiles(1) - 1.5*IQR;
            lowhisker = max(lowhisker, min(data(data > lowhisker)));
            hiwhisker = quartiles(3) + 1.5*IQR;
            hiwhisker = min(hiwhisker, max(data(data < hiwhisker)));
            if ~isempty(lowhisker) && ~isempty(hiwhisker)
                obj.WhiskerPlot = plot([pos pos], [lowhisker hiwhisker]);
            end

            % Median
            obj.MedianPlot = scatter(pos, quartiles(2), args.MedianMarkerSize, [1 1 1], 'filled');

            % Notches
            obj.NotchPlots = ...
                scatter(pos, quartiles(2)-1.57*IQR/sqrt(length(data)), ...
                [], [1 1 1], 'filled', '^');
            obj.NotchPlots(2) = ...
                scatter(pos, quartiles(2)+1.57*IQR/sqrt(length(data)), ...
                [], [1 1 1], 'filled', 'v');

            %% Set graphical preferences
            obj.EdgeColor = args.EdgeColor;
            obj.MedianPlot.LineWidth = args.LineWidth;
            obj.BoxColor = args.BoxColor;
            obj.BoxWidth = args.BoxWidth;
            obj.MedianColor = args.MedianColor;
            obj.ShowData = args.ShowData;
            obj.ShowNotches = args.ShowNotches;
            obj.ShowMean = args.ShowMean;
            obj.ShowBox = args.ShowBox;
            obj.ShowMedian = args.ShowMedian;
            obj.ShowWhiskers = args.ShowWhiskers;

            if not(isempty(args.ViolinColor))
                if size(args.ViolinColor{1},1) > 1
                    ViolinColor{1} = args.ViolinColor{1}(pos,:);
                else
                    ViolinColor{1} = args.ViolinColor{1};
                end
                if length(args.ViolinColor)==2
                    if size(args.ViolinColor{2},1) > 1
                        ViolinColor{2} = args.ViolinColor{2}(pos,:);
                    else
                        ViolinColor{2} = args.ViolinColor{2};
                    end
                else
                    ViolinColor{2} = ViolinColor{1};
                end
            else
                % defaults
                if args.scpltBool
                    ViolinColor{1} = obj.ScatterPlot.CData;
                else
                    ViolinColor{1} = [0 0 0];
                end
                ViolinColor{2} = [0 0 0];
            end
            obj.ViolinColor = ViolinColor;

            if not(isempty(args.ViolinAlpha))
                if length(args.ViolinAlpha{1})>1
                    error('Only scalar values are accepted for the alpha color channel');
                else
                    ViolinAlpha{1} = args.ViolinAlpha{1};
                end
                if length(args.ViolinAlpha)==2
                    if length(args.ViolinAlpha{2})>1
                        error('Only scalar values are accepted for the alpha color channel');
                    else
                        ViolinAlpha{2} = args.ViolinAlpha{2};
                    end
                else
                    ViolinAlpha{2} = ViolinAlpha{1}/2;  % default unless specified
                end
            else
                % default
                ViolinAlpha = {1,1};
            end
            obj.ViolinAlpha = ViolinAlpha;

        end

        %% SET METHODS
        function set.EdgeColor(obj, color)
            if ~isempty(obj.ViolinPlot)
                obj.ViolinPlot.EdgeColor = color;
                obj.ViolinPlotQ.EdgeColor = color;
                if ~isempty(obj.ViolinPlot2)
                    obj.ViolinPlot2.EdgeColor = color;
                end
            end
        end

        function color = get.EdgeColor(obj)
            if ~isempty(obj.ViolinPlot)
                color = obj.ViolinPlot.EdgeColor;
            end
        end

        function set.MedianColor(obj, color)
            obj.MedianPlot.MarkerFaceColor = color;
            if ~isempty(obj.NotchPlots)
                obj.NotchPlots(1).MarkerFaceColor = color;
                obj.NotchPlots(2).MarkerFaceColor = color;
            end
        end

        function color = get.MedianColor(obj)
            color = obj.MedianPlot.MarkerFaceColor;
        end

        function set.BoxColor(obj, color)
            if ~isempty(obj.BoxPlot)
                obj.BoxPlot.FaceColor = color;
                obj.BoxPlot.EdgeColor = color;
                obj.WhiskerPlot.Color = color;
                obj.MedianPlot.MarkerEdgeColor = color;
                obj.NotchPlots(1).MarkerFaceColor = color;
                obj.NotchPlots(2).MarkerFaceColor = color;
            elseif  ~isempty(obj.ViolinPlotQ)
                obj.WhiskerPlot.Color = color;
                obj.MedianPlot.MarkerEdgeColor = color;
                obj.NotchPlots(1).MarkerFaceColor = color;
                obj.NotchPlots(2).MarkerFaceColor = color;
            end
        end

        function color = get.BoxColor(obj)
            if ~isempty(obj.BoxPlot)
                color = obj.BoxPlot.FaceColor;
            end
        end

        function set.BoxWidth(obj,width)
            if ~isempty(obj.BoxPlot)
                pos=mean(obj.BoxPlot.XData);
                obj.BoxPlot.XData=pos+[-1,1,1,-1]*width;
            end
        end

        function width = get.BoxWidth(obj)
            width=max(obj.BoxPlot.XData)-min(obj.BoxPlot.XData);
        end

        function set.ViolinColor(obj, color)
            obj.ViolinPlot.FaceColor = color{1};
            obj.ScatterPlot.MarkerFaceColor = color{1};
            obj.MeanPlot.Color = color{1};
            if ~isempty(obj.ViolinPlot2)
                obj.ViolinPlot2.FaceColor = color{2};
                obj.ScatterPlot2.MarkerFaceColor = color{2};
            end
            if  ~isempty(obj.ViolinPlotQ)
                obj.ViolinPlotQ.FaceColor = color{1};
            end
            for idx = 1: size(obj.HistogramPlot,1)
                obj.HistogramPlot(idx).Color = color{1};
            end
        end

        function color = get.ViolinColor(obj)
            color{1} = obj.ViolinPlot.FaceColor;
            if ~isempty(obj.ViolinPlot2)
                color{2} = obj.ViolinPlot2.FaceColor;
            end
        end

        function set.ViolinAlpha(obj, alpha)
            obj.ViolinPlotQ.FaceAlpha = .65;
            obj.ViolinPlot.FaceAlpha = alpha{1};
            obj.ScatterPlot.MarkerFaceAlpha = 1;
            if ~isempty(obj.ViolinPlot2)
                obj.ViolinPlot2.FaceAlpha = alpha{2};
                obj.ScatterPlot2.MarkerFaceAlpha = 1;
            end
        end

        function alpha = get.ViolinAlpha(obj)
            alpha{1} = obj.ViolinPlot.FaceAlpha;
            if ~isempty(obj.ViolinPlot2)
                alpha{2} = obj.ViolinPlot2.FaceAlpha;
            end
        end

        function set.ShowData(obj, yesno)
            if yesno
                obj.ScatterPlot.Visible = 'on';
                for idx = 1: size(obj.HistogramPlot,1)
                    obj.HistogramPlot(idx).Visible = 'on';
                end
            else
                obj.ScatterPlot.Visible = 'off';
                for idx = 1: size(obj.HistogramPlot,1)
                    obj.HistogramPlot(idx).Visible = 'off';
                end
            end
            if ~isempty(obj.ScatterPlot2)
                obj.ScatterPlot2.Visible = obj.ScatterPlot.Visible;
            end

        end

        function yesno = get.ShowData(obj)
            if ~isempty(obj.ScatterPlot)
                yesno = strcmp(obj.ScatterPlot.Visible, 'on');
            end
        end

        function set.ShowNotches(obj, yesno)
            if ~isempty(obj.NotchPlots)
                if yesno
                    obj.NotchPlots(1).Visible = 'on';
                    obj.NotchPlots(2).Visible = 'on';
                else
                    obj.NotchPlots(1).Visible = 'off';
                    obj.NotchPlots(2).Visible = 'off';
                end
            end
        end

        function yesno = get.ShowNotches(obj)
            if ~isempty(obj.NotchPlots)
                yesno = strcmp(obj.NotchPlots(1).Visible, 'on');
            end
        end

        function set.ShowMean(obj, yesno)
            if ~isempty(obj.MeanPlot)
                if yesno
                    obj.MeanPlot.Visible = 'on';
                else
                    obj.MeanPlot.Visible = 'off';
                end
            end
        end

        function yesno = get.ShowMean(obj)
            if ~isempty(obj.BoxPlot)
                yesno = strcmp(obj.BoxPlot.Visible, 'on');
            end
        end

        function set.ShowBox(obj, yesno)
            if ~isempty(obj.BoxPlot)
                if yesno
                    obj.BoxPlot.Visible = 'on';
                else
                    obj.BoxPlot.Visible = 'off';
                end
            end
        end

        function yesno = get.ShowBox(obj)
            if ~isempty(obj.BoxPlot)
                yesno = strcmp(obj.BoxPlot.Visible, 'on');
            end
        end

        function set.ShowMedian(obj, yesno)
            if ~isempty(obj.MedianPlot)
                if yesno
                    obj.MedianPlot.Visible = 'on';
                else
                    obj.MedianPlot.Visible = 'off';
                end
            end
        end

        function yesno = get.ShowMedian(obj)
            if ~isempty(obj.MedianPlot)
                yesno = strcmp(obj.MedianPlot.Visible, 'on');
            end
        end

        function set.ShowWhiskers(obj, yesno)
            if ~isempty(obj.WhiskerPlot)
                if yesno
                    obj.WhiskerPlot.Visible = 'on';
                else
                    obj.WhiskerPlot.Visible = 'off';
                end
            end
        end

        function yesno = get.ShowWhiskers(obj)
            if ~isempty(obj.WhiskerPlot)
                yesno = strcmp(obj.WhiskerPlot.Visible, 'on');
            end
        end

    end

    methods (Access=private)
        function results = checkInputs(~, data, pos, varargin)
            isscalarnumber = @(x) (isnumeric(x) & isscalar(x));
            p = inputParser();
            p.addRequired('Data', @(x)isnumeric(vertcat(x{:})));
            p.addRequired('Pos', isscalarnumber);
            p.addParameter('Width', 0.3, isscalarnumber);
            p.addParameter('Bandwidth', [], isscalarnumber);
            iscolor = @(x) (isnumeric(x) & size(x,2) == 3);
            p.addParameter('ViolinColor', [], @(x)iscolor(vertcat(x{:})));
            p.addParameter('MarkerSize', 24, @isnumeric);
            p.addParameter('MedianMarkerSize', 36, @isnumeric);
            p.addParameter('LineWidth', 0.75, @isnumeric);
            p.addParameter('BoxColor', [0.5 0.5 0.5], iscolor);
            p.addParameter('BoxWidth', 0.01, isscalarnumber);
            p.addParameter('EdgeColor', [0.5 0.5 0.5], iscolor);
            p.addParameter('MedianColor', [1 1 1], iscolor);
            p.addParameter('ViolinAlpha', {0.3,0.15}, @(x)isnumeric(vertcat(x{:})));
            isscalarlogical = @(x) (islogical(x) & isscalar(x));
            p.addParameter('ShowData', true, isscalarlogical);
            p.addParameter('ShowNotches', false, isscalarlogical);
            p.addParameter('ShowMean', false, isscalarlogical);
            p.addParameter('ShowBox', true, isscalarlogical);
            p.addParameter('ShowMedian', true, isscalarlogical);
            p.addParameter('ShowWhiskers', true, isscalarlogical);
            validSides={'full', 'right', 'left'};
            checkSide = @(x) any(validatestring(x, validSides));
            p.addParameter('HalfViolin', 'full', checkSide);
            validQuartileStyles={'boxplot', 'shadow', 'none'};
            checkQuartile = @(x)any(validatestring(x, validQuartileStyles));
            p.addParameter('QuartileStyle', 'boxplot', checkQuartile);
            validDataStyles = {'scatter', 'histogram', 'none'};
            checkStyle = @(x)any(validatestring(x, validDataStyles));
            p.addParameter('DataStyle', 'scatter', checkStyle);

            p.parse(data, pos, varargin{:});
            results = p.Results;
        end
    end

    methods (Static)
        function [density, value, width] = calcKernelDensity(data, bandwidth, width)
            if isempty(data)
                error('Empty input data');
            end
            [density, value] = ksdensity(data, 'bandwidth', bandwidth);
            density = density(value >= min(data) & value <= max(data));
            value = value(value >= min(data) & value <= max(data));
            value(1) = min(data);
            value(end) = max(data);
            value = [value(1)*(1-1E-5), value, value(end)*(1+1E-5)];
            density = [0, density, 0];

            % all data is identical
            if min(data) == max(data)
                density = 1;
                value= mean(value);
            end

            width = width/max(density);
        end
    end
end

author: bastibe

repo: Violinplot-Matlab


text:
function err = procstreamOrderCheckDlg(elem)
    % Check the processing stream order
    [fn_error, missing_args, prereqs] = elem.CheckProcStreamOrder();
    if isa(fn_error, 'FuncCallClass')
        l1 = sprintf('The following function: %s', fn_error.nameUI);
        l2 = sprintf('cannot run because of unavailable input(s) %s.\n', cell2str(missing_args));
        if ~isempty(prereqs)
           l3 = sprintf('Add one of the following prerequisite functions to the processing stream:\n%s', strtrim(prereqs));
        else
           l3 = 'Ensure that a function which outputs the necessary inputs appears before the function in the processing stream.';
        end
        choice = questdlg({l1, l2, l3}, ...
        'Invalid Processing Stream', ...
        'Continue Anyway','Cancel', 'Cancel');
        switch choice
            case 'Continue Anyway'
                err = 0;
            case 'Cancel'
                err = -1;
            otherwise
                err = -1;
        end
    else
       err = 0;
    end
end

author: BUNPC


text:
%book : Signals and Systems Laboratory with MATLAB
%by Alex Palamides & Anastasia Veloni

% FS of periodic signals

%x(t)= 1, 0<t<1   , T=2
%     -1, 1<t<2

%x(t), 0<t<5T
t1=0:.01:1;
t2=1.01:.01:2;
x1=ones(size(t1));
x2=-ones(size(t2));
x=[x1 x2];
xp=repmat(x,1,5);
t=linspace(0,10,length(xp));
plot(t,xp)

%x(t), 0<t<T
figure
syms t
x=heaviside(t)-2*heaviside(t-1);
ezplot(x,[0 2]);

%exponential FS
figure
k=-2:2;
t0 =0;
T=2;
w=2*pi/T;
a=(1/T)*int(x*exp(-j*k*w*t),t,t0,t0+T)
xx=sum(a.*exp(j*k*w*t))
ezplot(xx,[0 10])
title('approximation with 5 terms')
figure
k=-5:5;
a=(1/T)*int(x*exp(-j*k*w*t),t,t0,t0+T);
xx=sum(a.*exp(j*k*w*t));
ezplot(xx,[0 10])
title('approximation with 11 terms')
figure
k=-10:10;
a=(1/T)*int(x*exp(-j*k*w*t),t,t0,t0+T);
xx=sum(a.*exp(j*k*w*t));
ezplot(xx,[0 10])
title('approximation with 21 terms')
figure
k=-30:30;
a=(1/T)*int(x*exp(-j*k*w*t),t,t0,t0+T);
xx=sum(a.*exp(j*k*w*t));
ezplot(xx,[0 10])
title('approximation with 61 terms')

%trigonometric FS
figure
a0=(1/T)*int(x,t0,t0+T);
n=1:2;
b=(2/T)*int(x*cos(n*w*t),t,t0,t0+T);
c=(2/T)*int(x*sin(n*w*t),t,t0,t0+T);
xx=a0+sum(b.*cos(n*w*t))+sum(c.*sin(n*w*t));
ezplot(xx,[0 10])
title('approximation with 3 terms')
figure
n=1:5;
b=(2/T)*int(x*cos(n*w*t),t,t0,t0+T);
c=(2/T)*int(x*sin(n*w*t),t,t0,t0+T);
xx=a0+sum(b.*cos(n*w*t))+sum(c.*sin(n*w*t));
ezplot(xx,[0 10])
title(' approximation with 6 terms')
figure
n=1:40;
b=(2/T)*int(x*cos(n*w*t),t,t0,t0+T);
c=(2/T)*int(x*sin(n*w*t),t,t0,t0+T);
xx=a0+sum(b.*cos(n*w*t))+sum(c.*sin(n*w*t));
ezplot(xx,[0 10])
title('approximation with 41 terms')

%FS in cosine with phase form
figure
for n=1:2
A(n)=sqrt(b(n)^2+c(n)^2);
thita(n)=atan2(-eval(c(n)),eval(b(n)));
end
k=1:2;
xx=a0+sum(A.*cos(k*w*t+thita))
ezplot(xx,[0 10])
title('approximation with 3 terms')
figure
for n=1:5
A(n)=sqrt(b(n)^2+c(n)^2);
thita(n)=atan2(-eval(c(n)),eval(b(n)));
end
k=1:5;
xx=a0+sum(A.*cos(k*w*t+thita));
ezplot(xx,[0 10])
title('approximation with 6 terms')
figure
for n=1:40
A(n)=sqrt(b(n)^2+c(n)^2);
thita(n)=atan2(-eval(c(n)),eval(b(n)));
end
k=1:40;
xx=a0+sum(A.*cos(k*w*t+thita));
ezplot(xx,[0 10])
title('approximation with 41 terms')

author: Sable

repo: mcbench-benchmarks


text:
%% function Xhat = triangulate(x1,x2,P2,homogeneous)
% triangulates the points Xhat, x1,x2 using P1(assumed to be [I|0]),P2
% Assumes the points x1,x2 are normalized so their last dimension is 1
% Xhat is the homogeneous coordinates of the reconstructed X
function Xhat = triangulate(x1,x2,P2,homogeneous)
Xhat = ones(4,size(x1,2));
P1 = [eye(3) zeros(3,1)];
for i = 1 : size(x1,2)
    if homogeneous
        A =[x1(1,i)*P1(3,:) - P1(1,:);
            x1(2,i)*P1(3,:) - P1(2,:);
            x2(1,i)*P2(3,:) - P2(1,:);
            x2(2,i)*P2(3,:) - P2(2,:)];
        [u,s,v] = svd(A);
        Xhat(:,i) = v(:,end);
    else
        A =[x1(1,i)*P1(3,1:3) - P1(1,1:3);
            x1(2,i)*P1(3,1:3) - P1(2,1:3);
            x2(1,i)*P2(3,1:3) - P2(1,1:3);
            x2(2,i)*P2(3,1:3) - P2(2,1:3)];
        b = -[x1(1,i)*P1(3,4) - P1(1,4);
            x1(2,i)*P1(3,4) - P1(2,4);
            x2(1,i)*P2(3,4) - P2(1,4);
            x2(2,i)*P2(3,4) - P2(2,4)];
        [u,s,v] = svd(A);
        bprime = u'*b;
        y = bprime(1:3) ./ (diag(s));
        Xhat(1:3,i) = v*y;
    end
end
Xhat = Xhat./repmat(Xhat(4,:),4,1);
end


text:
function [h_polygons, h_circles] = createStadium(mode)

    %% ----- Initialize the obstacles and the arena

    if( strcmp(mode,'Static') )

        obs = [];
        h_polygons = [];
        h_circles = [];

        %Properties of the arena
        arena_r = 25;                           %The arena radius
        [X,Y] = circle([0, 0], arena_r, 1000);
        arena = patch(X,Y,'w');

        % Properties of the obstacles
        n_points = 20;                          %number of points the obstacles are divided into
        %Offset of the origo (0,0)
        offs = 11;

        % Radius of the circle obstacles
        radius = [2.5, 2.5, 4.5, 2.5];

        % Origo of each circle obstacles
        centre = [offs+2, 8; ...
                -offs+1, -offs+1; ...
                0, 25; -2, 13];

        % Create the circle obstacles
        n_circles = size(centre,1);             %number of obstacles
        for i = 1:n_circles
            [X, Y] = circle(centre(i,:), radius(i), n_points);
            obs(i) = patch(X, Y, 'r');
        end

        %check if there already exists circle obstacles
        if( ~isempty(obs))
            h_circles = [obs, arena];

        % Set h_circles equal the arena no circle obstacles exists
        elseif( isempty(obs) )
            h_circles = arena;

        end

        %square = [10, 0; 15, -4; 15, 4; 10, 2];
        %square1 = [10, -22.9; 15, -20; 8, -11.38; 3, -14.28];
        %square1 = [16, -19.21; 18, -17.35; 13, -11.97; 11, -13.83];
        %square1 = [17, -18.33; 19, -16.25; 12, -9.52; 10, -11.6];
        %triangle1 = [16, -19.21; 18, -17.35; 13, -13.923];
        triangle1 = [10, -22.9; 15, -20; 7, -11.967];
        square2 = [-24.92, 2; -24.92, -2; -15, -2; -15, 2];

        %sq1 = patch(square1(:,1), square1(:,2), 'g');
        tri1 = patch(triangle1(:,1), triangle1(:,2), 'g');
        sq2 = patch(square2(:,1), square2(:,2), 'g');

        h_polygons = [tri1, sq2];

    elseif( strcmp(mode,'Dynamic') )

        obs = [];
        h_polygons = [];
        h_circles = [];

        arena_r = 25;
        [X,Y] = circle([0, 0], arena_r, 1000);
        arena = patch(X,Y,'w');
        n_points = 20;
        offs = 11;

        %radius = [2.5, 2.5, 2.5, 5.5];
        radius = [1.8, 1.8, 1.8, 3.5, 1.8, 1.8];
        centre = [0, norm([offs,offs]); ...
                -offs, -offs;  ...
                offs, -offs;  ...
                0, 25; -8, 0; 8, 0];

        n_circles = size(centre,1);
        for i = 1:n_circles
            [X, Y] = circle(centre(i,:), radius(i), n_points);
            obs(i) = patch(X, Y, 'r');
        end

        if( ~isempty(obs))
            h_circles = [obs, arena];

        elseif( isempty(obs) )
            h_circles = arena;

        end

        %square = [10, 0; 15, -4; 15, 4; 10, 2];
        %square1 = [10, -22.9; 15, -20; 8, -11.38; 3, -14.28];
        %square1 = [16, -19.21; 18, -17.35; 13, -11.97; 11, -13.83];
        %square1 = [17, -18.33; 19, -16.25; 12, -9.52; 10, -11.6];
        triangle1 = [16, -19.21; 18, -17.35; 13, -13.923];
        %triangle1 = [10, -22.9; 15, -20; 7, -11.967];
        %square2 = [-24.92, 2; -24.92, -2; -19, -2; -19, 2];
        square2 = [-24.92, 2; -24.92, -2; -18, -2; -18, 2];

        %sq1 = patch(square1(:,1), square1(:,2), 'g');
        tri1 = patch(triangle1(:,1), triangle1(:,2), 'g');
        sq2 = patch(square2(:,1), square2(:,2), 'g');

        h_polygons = [tri1, sq2];
    end
end

author: kennydl

repo: Reinforcment-Learning-With-Q-Learning


text:
function [depths,parents] = gsp_tree_depths(A,root)

if gsp_check_connectivity(A) == 0
    error('Graph is not connected');
end

N = size(A,1);
assigned = root;
depths = zeros(N,1);
parents = zeros(N,1);

next_to_expand = root;
current_depth = 1;

while ( length(assigned) < N )
    new_entries_whole_round = [];
    for i = 1:length(next_to_expand)
        neighbors = find(A(next_to_expand(i),:)>1e-7);
        new_entries = setdiff(neighbors,assigned);
        parents(new_entries) = next_to_expand(i);
        depths(new_entries)=current_depth;
        assigned=[assigned;new_entries'];
        new_entries_whole_round=[new_entries_whole_round;new_entries'];
    end
    current_depth=current_depth+1;
    next_to_expand=new_entries_whole_round;
end

end

repo: gspbox


text:
% make the structure of an embedded HMM with 2 rows and 3 columns

% 1------------>2
% |\   \        | \  \
% 3->4->5       6->7->8

n = 8;
dag = zeros(n);
dag(1,[2 3 4 5])=1;
dag(2,[6 7 8])=1;
for i=3:4
  dag(i,i+1)=1;
end
for i=6:7
  dag(i,i+1)=1;
end
ns = 2*ones(1,n);
bnet = mk_bnet(dag,ns);
for i=1:n
  bnet.CPD{i}=tabular_CPD(bnet,i);
end
[jtree, root, cliques] =  graph_to_jtree(moralize(bnet.dag), ones(1,n), {}, {});
%[jtree, root, cliques, B, w, elim_order, moral_edges, fill_in_edges] = dag_to_jtree(bnet);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%% Testing OOQP
clc
clear

%% QP1 -2.83333333301227;
H = speye(3);
f = -[2 3 1]';
A = sparse([1 1 1;3 -2 -3; 1 -3 2]);
b = [1;1;1];

opts = [];
opts.display = 2;

[x,fval,e,i,l] = ooqp(H,f,A,-Inf(size(b)),b,[],[],[],[],opts)

%% QP2 -8.22222220552525
clc
H = (sparse([1 -1; -1 2]));
f = -[2 6]';
A = sparse([1 1; -1 2; 2 1]);
b = [2; 2; 3];
lb = [0;0];

opts = [];
opts.display = 2;

[~,fval,e,i] = ooqp(tril(H),f,A,-Inf(size(b)),b,[],[],lb,[],opts)

%% QP3 -6.41379310344827
H = tril(sparse([1 -1; -1 2]));
f = -[2 6]';
A = sparse([1 1; -1 2; 2 1]);
b = [2; 2; 3];
Aeq = sparse([1 1.5]);
beq = 2;
lb = [0;0];
ub = [10;10];

opts = [];
opts.display = 2;
opts.linear_solver = 'pardiso';
[x,f,e,i] = ooqp(H,f,A,-Inf(size(b)),b,Aeq,beq,lb,ub,opts)

%% LP1 -3.75
clc
H = [];
f = -[-1, 2]';
A = sparse([2, 1;-4, 4]);
b = [5, 5]';

opts = [];
opts.display = 2;
[x,f,e,i] = ooqp(H,f,A,-Inf(size(b)),b,[],[],[],[],opts)

%% LP2 3
clc
H = [];
f = [1, 2, 3, 7, 8, 8];
A = -sparse([5, -3, 2, -3, -1, 2; -1, 0, 2, 1, 3, -3;1, 2, -1, 0, 5, -1]);
b = -[-5, -1, 3]';
lb = zeros(6,1);
ub = 10*ones(6,1);

opts = [];
opts.display = 2;
[x,f,e,i,l] = ooqp(H,f,-A,-Inf(size(b)),-b,[],[],lb,ub,opts)

%% LP3 4
clc
n = 40;
t = (0:n-1)';
y = 3.5 -.2*t;
b = y + 0.5*ones(size(y));
m = [ones(n,1),t(:)];
A = sparse([m,-m,eye(n)]);
H = [];spalloc(44,44,0);
f = [sum(m),sum(-m),2*ones(1,n)];
lb = zeros(n+4,1);
ub = [10, 10, 10, 10, 5*ones(1,n)];

opts = [];
opts.display = 2;
[~,f,e,i,l] = ooqp(H,f,-A,-Inf(size(b)),-b,[],[],lb,ub,opts)

%% LARGE LP
clc
prob = coinRead('maros-r7.mps');
opts = optiset('solver','ooqp','display','iter');
Opt = opti(prob,opts);
%Solve
[~,f,e,i] = solve(Opt)

%%
qp1 = qp_prob(1);
opts = optiset('solver','ooqp','display','iter');
Opt = opti(qp1,opts);
%Solve
[~,f,e,i] = solve(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function log_info( input_str )
time_str = datestr(now, 31);
info_str = strcat('[', time_str, '] ', input_str);
disp(info_str);
end

author: XgTu


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

%% -*- texinfo -*-
%% @defmethod  @@symfun private_disp_name (@var{f}, @var{name})
%% A string appropriate for representing the name of this symfun.
%%
%% Private method: this is not the method you are looking for.
%%
%% @end defmethod

function s = private_disp_name(f, input_name)

  if (isempty(input_name))
    s = input_name;
    return
  end

  vars = f.vars;
  if length(vars) == 0
    varstr = '';
  else
    v = vars{1};
    varstr = v.flat;
  end
  for i = 2:length(vars);
    v = vars{i};
    varstr = [varstr ', ' v.flat];
  end
  s = [input_name, '(', varstr, ')'];

end

%!test
%! syms f(x)
%! s = private_disp_name(f, 'f');
%! assert (strcmp (s, 'f(x)'))

%!test
%! syms x y
%! g(y, x) = x + y;
%! s = private_disp_name(g, 'g');
%! assert (strcmp (s, 'g(y, x)'))

%!test
%! syms f(x)
%! assert (isempty (private_disp_name(f, '')))

repo: octsympy


text:
function M = xfm_read(fname)
% M = xfm_read(fname)

%
% xfm_read.m
%
% Original Author: Bruce Fischl
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:13 $
%    $Revision: 1.4 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

fid = fopen(fname) ;
if (fid < 0)
	error(sprintf('could not open file %s', fname));
end

tline = fgetl(fid) ;
while ((length(tline) > 0) & (tline(1) == '%'))
	tline = fgetl(fid) ;
end

tok = strtok(tline);
while (strcmp(tok, 'Linear_Transform') ~= 1)
	tline = fgetl(fid) ;
	tok = strtok(tline);
end

M = zeros(4,4) ; M(4,4) = 1;
for row=1:3
	tline = fgetl(fid) ;  % one row of matrix
	tmp = sscanf(tline, '%f');
	 M(row,:) = tmp';
end

fclose(fid) ;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function output=fcn_LeaderAcceleration(u)
global leaderNum dim samplePos_all nodenum formationType
global rotateFlag tempp2 tempp3

u_all = reshape(u(1:dim*nodenum*3),dim,nodenum*3);
v_all=u_all(:,1:nodenum);
p_all=u_all(:,nodenum+1:2*nodenum);
control_all=u_all(:,2*nodenum+1:3*nodenum);
currentTime=u(end);
%output
a_leader_all=zeros(dim,leaderNum);

% speed up to speed=1
time_start_speedup=0;
time_span_speedup=5;
if currentTime>time_start_speedup && currentTime<time_start_speedup+time_span_speedup
    ax=fcn_sin(time_start_speedup,time_span_speedup,currentTime)/pi;
    a_leader_all(:,1)=[ax,0]';
    a_leader_all(:,2)=[ax,0]';
    a_leader_all(:,3)=[ax,0]';
end

% slow down to speed=0
time_start_slowdown=10;
time_span_slowdown=5;
if currentTime>time_start_slowdown && currentTime<time_start_slowdown+time_span_slowdown
    %ax=fcn_normalDistribution(time_start_slowdown,time_span_slowdown, 2, currentTime);
    ax=-fcn_sin(time_start_slowdown,time_span_slowdown,currentTime)/pi;
    a_leader_all(:,1)=[ax,0]';
    a_leader_all(:,2)=[ax,0]';
    a_leader_all(:,3)=[ax,0]';
end

% rotate
v3=v_all(:,3);
v2=v_all(:,2);
v1=v_all(:,1);
p3=p_all(:,3);
p2=p_all(:,2);
p1=p_all(:,1);
R90=[cos(pi/2),-sin(pi/2);
     sin(pi/2),cos(pi/2)];
time_start_rotate=20;
time_span_rotate=10;
if currentTime>=time_start_rotate && currentTime<=time_start_rotate+time_span_rotate
    if rotateFlag==1
        tempp2=p2;
        tempp3=p3;
        rotateFlag=0;
    end
    % agent 2
    acce_mag=-fcn_SlowDownSpeedUp(time_start_rotate, time_span_rotate, currentTime, 4);
%     acce_vec=(p1+R90'*(tempp2-p1))-tempp2;
%     acce_vec=acce_vec/norm(acce_vec);
    acce_vec=[1 0]';
    a_leader_all(:,2)=acce_mag*acce_vec;
    % agent 3
    acce_mag=-fcn_SlowDownSpeedUp(time_start_rotate, time_span_rotate, currentTime, 4);
%     acce_vec=(p1+R90'*(tempp3-p1))-tempp3;
%     acce_vec=acce_vec/norm(acce_vec);
    acce_vec=[0 1]';
    a_leader_all(:,3)=acce_mag*acce_vec;
end

% move downward
time_start_downward=time_start_rotate;
time_span_downward=time_span_rotate;
if currentTime>=time_start_downward && currentTime<=time_start_downward+time_span_downward
    % agent 1
    ay=fcn_sin(time_start_downward,time_span_downward,currentTime)/pi;
    a_leader_all(:,1)=[0,-ay]';
    % agent 2
    a_leader_all(:,2)=a_leader_all(:,2)+[0,-ay]';
    % agent 3
    a_leader_all(:,3)=a_leader_all(:,3)+[0,-ay]';
end

% slow down to speed=0
time_start_slowdown2=time_start_downward+time_span_downward+5;
time_span_slowdown2=5;
if currentTime>time_start_slowdown2 && currentTime<time_start_slowdown2+time_span_slowdown2
    %ax=fcn_normalDistribution(time_start_slowdown,time_span_slowdown, 2, currentTime);
    ay=-fcn_sin(time_start_slowdown2,time_span_slowdown2,currentTime)/pi;
    a_leader_all(:,1)=[0,-ay]';
    a_leader_all(:,2)=[0,-ay]';
    a_leader_all(:,3)=[0,-ay]';
end

% if currentTime>=time_start_rotate && currentTime<=time_start_rotate+time_span_rotate
%     % agent 3
%     aUnit=(p1-p3)/norm(p1-p3);
%     tht_dot=fcn_sin(time_start_rotate,time_span_rotate,currentTime);
%     speed_dot=fcn_sin(time_start_rotate,time_span_rotate,currentTime);
%     a_leader_all(:,3)=aUnit*tht_dot*norm(v3)+speed_dot*R90*aUnit;
%     % agent 2
%     aUnit=(p1-p2)/norm(p1-p2);
%     %tht_dot=fcn_sin(time_start_rotate,time_span_rotate,currentTime);
%     speed_dot=fcn_sin(time_start_rotate,time_span_rotate,currentTime);
%     a_leader_all(:,2)=aUnit*tht_dot*norm(v2)+speed_dot*R90*aUnit;
% end

% % scaling
% a_leader_all(2,2)=-fcn_TanhManeuver(20, 1, 10, currentTime);
% a_leader_all(2,3)=fcn_TanhManeuver(20, 1, 10, currentTime);

% % rotation
% v3=v_all(:,3);
% v2=v_all(:,2);
% v1=v_all(:,1);
% p3=p_all(:,3);
% p2=p_all(:,2);
% p1=p_all(:,1);
% R90=[cos(pi/2),-sin(pi/2);
%      sin(pi/2),cos(pi/2)];
% time_start_rotate=20;
% % agent 3
% time_span_rotate3=10;
% if currentTime>=time_start_rotate && currentTime<=time_start_rotate+time_span_rotate3
%     if rotateCenterFlag==1
%         roateCenter(:,3)=p3-[0,5]';
%         rotateCenterFlag=0;
%         r3=norm(roateCenter(:,3)-p3);
%         roateCenter(:,2)=roateCenter(:,3);
%         r2=norm(roateCenter(:,2)-p2);
%         roateCenter(:,1)=roateCenter(:,3);
%         r1=norm(roateCenter(:,1)-p1);
%     end
%     tht_dot=pi/time_span_rotate3;
%     speed_dot=0;
%     a_leader_all(:,3)=R90'*v3/norm(v3)*tht_dot*norm(v3)+speed_dot*v3/norm(v3);
% end
% % agent 2
% time_span_rotate2=time_span_rotate3;
% if currentTime>=time_start_rotate && currentTime<=time_start_rotate+time_span_rotate2
%     tht_dot=pi/time_span_rotate2;
%     if currentTime>=time_start_rotate && currentTime<=time_start_rotate+1
%         speed_dot=fcn_sin(time_start_rotate,1,currentTime)/2.6;
%     else
%         speed_dot=0;
%     end
%     a_leader_all(:,2)=R90'*v2/norm(v2)*tht_dot*norm(v2)+speed_dot*v2/norm(v2);
% end
% if currentTime>=time_start_rotate+time_span_rotate2 && currentTime<=time_start_rotate+time_span_rotate2+1
%         speed_dot=-fcn_sin(time_start_rotate,1,currentTime)/2.6;
%         tht_dot=0;
%         a_leader_all(:,2)=R90'*v2/norm(v2)*tht_dot*norm(v2)+speed_dot*v2/norm(v2);
% end
% % agent 1
% time_span_rotate1=time_span_rotate3;
% if currentTime>=time_start_rotate && currentTime<=time_start_rotate+time_span_rotate1
%     tht_dot=pi/time_span_rotate3;%fcn_sin(time_start_rotate,time_span_rotate1,currentTime);
%     if currentTime>=time_start_rotate && currentTime<=time_start_rotate+1
%         speed_dot=fcn_sin(time_start_rotate,1,currentTime)/3;
% %         tht_dot=
%     else
%         speed_dot=0;
%     end
%     aUnit=(roateCenter(:,1)-p1)/norm(roateCenter(:,1)-p1);
%     a_leader_all(:,1)=aUnit*tht_dot*norm(v1)+speed_dot*v1/norm(v1);
% end
% if currentTime>=time_start_rotate+time_span_rotate2 && currentTime<=time_start_rotate+time_span_rotate2+1
%         speed_dot=-fcn_sin(time_start_rotate,1,currentTime)/2.6;
%         tht_dot=0;
%         a_leader_all(:,2)=R90'*v2/norm(v2)*tht_dot*norm(v2)+speed_dot*v2/norm(v2);
% end

% % agent 1: distance-based formation control law
% time_span_rotate1=time_span_rotate3;
% if currentTime>=time_start_rotate && currentTime<=time_start_rotate+time_span_rotate1
%     delta12=norm(p1-p2)-2*sqrt(2);
%     delta13=norm(p1-p3)-2*sqrt(2);
%     a_leader_all(:,1)=5*(delta12*(p2-p1)+delta13*(p3-p1)) + 0.5*(delta12*(v2-v1)+delta13*(v3-v1));
% end

% % agent 2
% time_span_rotate2=26;
% if currentTime>=time_start_rotate && currentTime<=time_start_rotate+time_span_rotate2
%     v2=v_all(:,2);
%     p2=p_all(:,2);
%     R90=[cos(pi/2),-sin(pi/2);
%              sin(pi/2),cos(pi/2)];
%     tht_dot=fcn_sin(time_start_rotate,time_span_rotate2,currentTime);
%     speed_dot=0;
%     a_leader_all(:,2)=R90'*v2/norm(v2)*tht_dot*norm(v2)+speed_dot*v2/norm(v2);
% end
% % agent 1
% time_span_rotate1=18;
% if currentTime>=time_start_rotate && currentTime<=time_start_rotate+time_span_rotate1
%     v1=v_all(:,1);
%     p1=p_all(:,1);
%     R90=[cos(pi/2),-sin(pi/2);
%              sin(pi/2),cos(pi/2)];
%     tht_dot=fcn_sin(time_start_rotate,time_span_rotate1,currentTime);
%     speed_dot=0;
%     a_leader_all(:,1)=R90'*v1/norm(v1)*tht_dot*norm(v1)+speed_dot*v1/norm(v1);
% end
% if currentTime>=time_start_rotate && currentTime<=time_start_rotate+time_span_rotate
%     v3=v_all(:,3);
%     v2=v_all(:,2);
%     v1=v_all(:,1);
%     p3=p_all(:,3);
%     p2=p_all(:,2);
%     p1=p_all(:,1);
%     R90=[cos(pi/2),-sin(pi/2);
%          sin(pi/2),cos(pi/2)];
%     % agent 3
%     tht_dot=-fcn_normalDistribution(time_start_rotate,time_span_rotate, 3.15, currentTime);
%     speed_dot=0;
%     a_leader_all(:,3)=R90'*v3/norm(v3)*tht_dot*norm(v3)+speed_dot*v3/norm(v3);
%     % agent 2
%     tht_dot=-fcn_normalDistribution(time_start_rotate,time_span_rotate*1.5, 3.15, currentTime);
%     speed_dot=0;
%     a_leader_all(:,2)=R90'*v2/norm(v2)*tht_dot*norm(v2)+speed_dot*v2/norm(v2);
%     % agent 1
%     tht_dot=-fcn_normalDistribution(time_start_rotate,time_span_rotate*, 2, currentTime);
%     speed_dot=0;%fcn_normalDistribution(time_start_rotate,time_span_rotate, 1, currentTime);
%     a_leader_all(:,1)=R90'*v1/norm(v1)*tht_dot*norm(v1)+speed_dot*v1/norm(v1);
% end

% if currentTime>=time_start_rotate && currentTime<=time_start_rotate+time_span_rotate
%     offset=4;
%     v3=v_all(:,3);
%     v2=v_all(:,2);
%     v1=v_all(:,1);
%     p3=p_all(:,3);
%     p2=p_all(:,2);
%     p1=p_all(:,1);
%     R90=[cos(pi/2),-sin(pi/2);
%          sin(pi/2),cos(pi/2)];
%     threshhold=0.005;
%     if rotateCenterFlag==1
%         roateCenter(:,3)=p3-[0,offset]';
%         rotateCenterFlag=0;
%         r3=norm(roateCenter(:,3)-p3);
%         rotateTime(3)=pi*r3/norm(v3);
%         roateCenter(:,2)=roateCenter(:,3);
%         r2=norm(roateCenter(:,2)-p2);
%         rotateTime(2)=pi*r2/norm(v2);
% %         roateCenter(:,1)=roateCenter(:,3)+(p1-p3);
%         roateCenter(:,1)=roateCenter(:,2)+[p1(1)-p2(1),0]';
%         r1=norm(roateCenter(:,1)-p1);
%         rotateTime(1)=pi*r1/norm(v1);
%     end
%     % agent 3
%     tht=atan2(v3(2),v3(1));
%     if ~(tht>pi-threshhold || tht<-pi+threshhold)
%         r3=norm(roateCenter(:,3)-p3);
%         aUnit=(roateCenter(:,3)-p3)/norm(roateCenter(:,3)-p3);
%         a_leader_all(:,3)=aUnit*norm(v3)^2/r3;%+speed_dot*v3/norm(v3);
%     end
%     % agent 2
%     tht=atan2(v2(2),v2(1));
%     if ~(tht>pi-threshhold || tht<-pi+threshhold)
%         r2=norm(roateCenter(:,2)-p2);
%         aUnit=(roateCenter(:,2)-p2)/norm(roateCenter(:,2)-p2);
%         speed_dot=-fcn_SlowDownSpeedUp(time_start_rotate, 5, currentTime, 8);
%         a_leader_all(:,2)=aUnit*norm(v2)^2/r2;%+speed_dot*v2/norm(v2);
%     end
%     % agent 1
%     tht=atan2(v1(2),v1(1));
%     if ~(tht>pi-threshhold || tht<-pi+threshhold)
%         r1=norm(roateCenter(:,1)-p1);
%         aUnit=(roateCenter(:,1)-p1)/norm(roateCenter(:,1)-p1);
%         speed_dot=-fcn_SlowDownSpeedUp(time_start_rotate, 5, currentTime, 8);
%         a_leader_all(:,1)=aUnit*norm(v1)^2/r1;%+speed_dot*v1/norm(v1);
%     end
% end

% % slow down and speed up to align the agents
% time_start_align=20;
% time_span=5;
% if currentTime>time_start_align && currentTime<time_start_align+time_span
% %     % agent 3: back 4 meters
% %     distance_change=4;
% %     a_leader_all(2,3)=-fcn_SlowDownSpeedUp(time_start_align, time_span, currentTime, distance_change);
% %     % agent 2: foward 2 meters
% %     distance_change=2.3;
% %     a_leader_all(2,2)=fcn_SlowDownSpeedUp(time_start_align, time_span, currentTime, distance_change);
%     % agent 1 slow down
%     distance_change=2;
%     a_leader_all(1,1)=fcn_SlowDownSpeedUp(time_start_align, time_span, currentTime, distance_change);
% end

% output: replace the leaders' velocity with the one prescribed above
control_all(:,1:leaderNum) = a_leader_all;
output=reshape(control_all, dim*nodenum, 1);

repo: paper-simulation


text:
function varargout = ft_warning(varargin)

% FT_WARNING prints a warning message on screen, depending on the verbosity
% settings of the calling high-level FieldTrip function. This function works
% similar to the standard WARNING function, but also features the "once" mode.
%
% Use as
%   ft_warning(...)
% with arguments similar to fprintf, or
%   ft_warning(msgId, ...)
% with arguments similar to warning.
%
% You can switch of all warning messages using
%   ft_warning off
% or for specific ones using
%   ft_warning off msgId
%
% To switch them back on, you would use
%   ft_warning on
% or for specific ones using
%   ft_warning on msgId
%
% Warning messages are only printed once per timeout period using
%   ft_warning timeout 60
%   ft_warning once
% or for specific ones using
%   ft_warning once msgId
%
% You can see the most recent messages and identifier using
%   ft_warning last
%
% You can query the current on/off/once state for all messages using
%   ft_warning query
%
% See also FT_ERROR, FT_WARNING, FT_NOTICE, FT_INFO, FT_DEBUG, ERROR, WARNING

% Copyright (C) 2012-2017, Robert Oostenveld, J?rn M. Horschig
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargout
  varargout{:} = ft_notification(varargin{:});
elseif isequal(varargin, {'last'})
  % return an answer anyway
  varargout{1} = ft_notification(varargin{:});
else
  ft_notification(varargin{:});
end

author: spm


text:
function Process_ConvertOptitrack2Pos(fbasename,varargin)

% USAGE
%
%   Process_ConvertOptitrack2Pos(fbasename,varargin)
%
%
%
%    fbasename   basename of the video file (should be filebasenmae.avi)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'syncDatFile' name of binary file where sync signal is stored
%                   (default = filebasename_digitalin.dat)
%     'syncSampFq'  sampling freqeuncy of the sync signal (default= 20kHz)
%     'syncChan'    sync channel to read (default = 1)
%     'syncNbCh'    number of channels in the sync file (default = 1)
%     'posSampFq'   sampling frequency of the final pos file (after
%                   interpolation)
%    =========================================================================
%
%
% This function assumes that a Motive (Optitrack) behavioral session
% has been exported to a .CSV file with default settings.
% (http://wiki.optitrack.com/index.php?title=Data_Export:_CSV)
%
%
%

if strcmp(fbasename,'')
    d = dir('*sessionInfo.mat');
    load(d.name);
    fbasename = sessionInfo.FileName;
end

warning('this functions is now deprecated and has been replaced by processConvertOptitrack2Behav.m')

syncDatFile = ['digitalin.dat']; % defaults if args aren't given... assuming a single digitalin channel on Intan
syncSampFq  = 20000;
syncChan    = 1;
syncNbCh    = 1;
posSampFq   = 120; %in Hz

% Parse options
for i = 1:2:length(varargin)
  if ~isa(varargin{i},'char')
    error(['Parameter ' num2str(i+3) ' is not a property .']);
  end
  switch(lower(varargin{i}))
    case 'syncdatfile'
      syncDatFile = varargin{i+1};
      if ~isa(duration,'char')
        error('Incorrect value for property ''syncDatFile'' .');
      end
    case 'syncsampfq'
      syncSampFq = varargin{i+1};
      if ~isa(frequency,'numeric')
        error('Incorrect value for property ''syncsampfq'' .');
      end
    case 'syncchan'
      syncChan = varargin{i+1};
      if ~isa(start,'numeric')pw
        error('Incorrect value for property ''syncChan'' .');
      end
        if start < 0, start = 0; end
    case 'syncnbch'
      syncNbCh = varargin{i+1};
      if ~isa(syncNbCh,'numeric')
        error('Incorrect value for property ''syncNbCh'' .');
      end
    case 'possampfq'
      posSampFq = varargin{i+1};
      if ~isa(posSampFq,'numeric')
        error('Incorrect value for property ''posSampFq'' .');
      end
    otherwise
      error(['Unknown property ''' num2str(varargin{i}) ''' .']);
  end
end

if exist([fbasename '.csv'],'file') % assumes csv naming matches basename
    dat = importdata([fbasename '.csv']);
    dat = scrubTracking(dat);
elseif ~isempty(dir('*.csv')) % looks for any csv file in the current recording folder
    csv = dir('*.csv');
    if length(csv) > 1
        dat=[];
       for i=1:length(csv)
           dat = [dat; importdata(csv(i).name)];
       end
    else
        dat = importdata(csv.name);
    end

    if isstruct(dat)
        d=[];
        for i=1:length(csv)
             d = [d;dat(i).data];
        end
        clear dat; dat.data = d;
    end
    dat = scrubTracking(dat);
elseif ~isempty(dir('Session*')) % Looks for a /Session*/ folder that Motive/Optitrack has created
    d = dir('Session*');
    cd(d.name)
    csv = dir('*.csv');
    if length(csv) > 1
        dat=[];
       for i=1:length(csv)
           dat = [dat; importdata(csv(i).name)];
       end
    else
        dat = importdata(csv.name);
    end
    if isstruct(dat)
        d=[];
        for i=1:length(csv)
             d = [d;dat(i).data];
        end
        clear dat; dat.data = d;
    else
        d.data = dat;
        dat = d;
    end
    dat = scrubTracking(dat);
    cd ..
end

pos = dat.data;

fid = fopen(syncDatFile);
dig = fread(fid,[syncNbCh inf],'int16=>int16');  % default type for Intan digitalin
dig = dig(syncChan,:);

t = (0:length(dig)-1)'/syncSampFq;

dPos = find(diff(dig)==1);
dNeg = find(diff(dig)==-1);

if length(dPos) == length(dNeg)+1
    dPos = dPos(1:end-1);
elseif length(dPos) > length(dNeg)+1 || length(dPos) < length(dNeg)
     keyboard
end

% Frame timing is the middle of shuter opening
frameT = (t(dPos)+t(dNeg))/2;

% The system sometimes (rarely) keeps on recording a few frames after software stopped
% recording. So we skip the last frames of the TTL

if length(frameT)<size(pos,1)
    warning('Too many video frames!')
    keyboard
elseif length(frameT)>size(pos,1)
    frameT = frameT(1:size(pos,1));
end

% We now interpolate the data at 120 Hz (or sampling fqcy specified in
% arguments)
% recDuration = length(dig)/syncSampFq;
%
% timestamps = (0:1/posSampFq:recDuration-1/posSampFq)';
% pos(pos==-1) = NaN;
% newPos = interp1(frameT,pos,timestamps,'linear');

% timestamps(isnan(newPos(:,2))) = [];
% newPos(isnan(newPos(:,2)),:)=[];

dlmwrite([fbasename '.pos'],[frameT pos],'delimiter','\t', 'precision', 32);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function incidence_matrix = gr_incidence_matrix ( node_num, ...
  node_coordinates, edge_num, edge_nodes )

%*****************************************************************************80
%
%% GR_INCIDENCE_MATRIX computes the incidence matrix.
%
%  Discussion:
%
%    The incidence matrix is of order EDGE_NUM by NODE_NUM.
%
%    A(I,J) = 1 if edge I uses node J.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_COORDINATES(2,NODE_NUM), the coordinates of the nodes.
%
%    Input, integer EDGE_NUM, the number of edges.
%
%    Input, integer EDGE_NODES(2,EDGE_NUM), the indices of the two nodes
%    that form each edge.
%
%    Output, integer INCIDENCE_MATRIX(EDGE_NUM,NODE_NUM), the incidence
%    matrix.
%
  incidence_matrix = zeros ( edge_num, node_num );

  for e = 1 : edge_num
    for i = 1 : 2
      n = edge_nodes(i,e);
      incidence_matrix(e,n) = incidence_matrix(e,n) + 1;
    end
  end

  return
end


text:
function[]=noxlabels(str)
%NOXLABELS   Remove some or all x-axis tick mark labels.
%
%   NOXLABELS by itself removes all tick mark labels on the current x-axis.
%
%   NOXLABELS MAX removes only the label of the maximum x-tick mark.
%
%   NOXLABELS MIN removes only the label of the minimum x-tick mark.
%
%   NOXLABELS ALL is the same as the default.
%
%   See also NOYLABELS.
%
%   Usage:  noxlabels
%           noxlabels xmin
%           noxlabels ymin
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004--2008 J.M. Lilly --- type 'help jlab_license' for details

if nargin==0
    str='all';
end

if strcmpi(str(1:3),'min')
    xlabels=get(gca,'xticklabel');
    xlabels(1,:)=' ';
    set(gca,'xticklabel',xlabels)
elseif strcmpi(str(1:3),'max')
    xlabels=get(gca,'xticklabel');
    xlabels(end,:)=' ';
    set(gca,'xticklabel',xlabels)
elseif strcmpi(str(1:3),'all')
    set(gca,'xticklabel',[])
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function h=m_gshhs_h(varargin)
% M_GSHHS_H Add a coastline to a given map using the 'high' resolution of
%           the Global Self-consistant Hierarchical High-resolution
%           Shorelines.
%
%         M_GSHHS_H( (standard line option,...,...) ) draws the coastline
%         as a simple line.
%         M_GSHHS_H('patch' ( ,standard patch options,...,...) ) draws the
%         coastline as a number of patches.
%
%         M_GSHHS_H('save',FILENAME) saves the extracted coastline data
%         for the current projection in a file FILENAME. This allows
%         speedier replotting using M_USERCOAST(FILENAME).
%
%         See also M_PROJ, M_GRID, M_COAST, M_GSHHS_L, M_GSHHS_I, M_GSHHS_C
%         M_USERCOAST

% Rich Pawlowicz (rich@ocgy.ubc.ca) 15/June/98
%
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.

% Changes
%  Nov/2017 - changed this into a stub calling gshhs.m (kept
%             for backwards compatability)

h=m_gshhs('hc',varargin{:});

author: goGPS-Project


text:
function [ loglik ] = ml_LogLikelihood_gmm(X,Priors,Mu,Sigma)
%LOGLIKELIHOOD_GMM
%
%   input------------------------------------------------------------------
%
%       o X : (D x N), set of N data points of dimension  D
%       o Priors:   1 x K array representing the prior probabilities of the
%               K GMM components.
%       o Mu:       D x K array representing the centers of the K GMM components.
%       o Sigma:    D x D x K array representing the covariance matrices of the
%                   K GMM components.
%
%   output ----------------------------------------------------------------
%
%       o loglik : (N x 1) , loglikelihood
%

N = size(X,2);
K = size(Priors,2);

Pxi = zeros(N,K);
for i=1:K
    Pxi(:,i) = ml_gaussPDF(X, Mu(:,i), Sigma(:,:,i));
end

%Compute the log likelihood
F            = Pxi*Priors';
F(F<realmin) = realmin;
F            = Pxi*Priors';
loglik       = sum(log(F));

end

author: epfl-lasa


text:
clear all; close all; clc

L=30; n=512; x2=linspace(-L/2,L/2,n+1); x=x2(1:n); % spatial discretization
k=(2*pi/L)*[0:n/2-1 -n/2:-1].';    % wavenumbers for FFT
V=x.^2.';             % potential
t=0:0.2:20;           % time domain collection points

u=exp(-0.2*(x-1).^2); % initial conditions
ut=fft(u);            % FFT initial data
[t,utsol]=ode45('harm_rhs',t,ut,[],k,V); % integrate PDE
for j=1:length(t)
  usol(j,:)=ifft(utsol(j,:));      % transforming back
end

u2=exp(-0.2*(x-0).^2);
ut=fft(u2);
[t,utsol]=ode45('harm_rhs',t,ut,[],k,V);
for j=1:length(t)
  usol2(j,:)=ifft(utsol(j,:));
end

figure(1)
subplot(2,2,2)
surfl(x,t,abs(usol)+2), shading interp, colormap(gray) %colormap([0 0 0])
hold on
pcolor(x,t,abs(usol)), shading interp, %colormap(gray)
view(20,25)
tv=0*x+20;
Vx=x.^2;
plot3(x(12:end-12),tv(12:end-12),Vx(12:end-12)/100+2,'k','Linewidth',[2])
set(gca,'Xlim',[-15 15],'Ylim',[0 20],'Zlim',[0 4],'Xtick',[-15 -10 -5 0 5 10 15],'Ytick',[0 10 20],'Ztick',[2 3 4], ...
    'Zticklabel',{'0','1','2'},'Fontsize',[15]);

figure(1)
subplot(2,2,1)
surfl(x,t,abs(usol2)+2), shading interp, colormap(gray) %colormap([0 0 0])
hold on
pcolor(x,t,abs(usol2)), shading interp, %colormap(gray)
view(20,25)
tv=0*x+20;
Vx=x.^2;
plot3(x(12:end-12),tv(12:end-12),Vx(12:end-12)/100+2,'k','Linewidth',[2])
set(gca,'Xlim',[-15 15],'Ylim',[0 20],'Zlim',[0 4],'Xtick',[-15 -10 -5 0 5 10 15],'Ytick',[0 10 20],'Ztick',[2 3 4], ...
    'Zticklabel',{'0','1','2'},'Fontsize',[15]);

%surfl(x,t,abs(usol)), colormap(gray)

for j=1:length(t)
    usol3(j,:)=usol(j,n:-1:1);
end

usym=[usol; usol3];

[U,S,V]=svd(usol.');
[U2,S2,V2]=svd(usol2.');
[U3,S3,V3]=svd(usym.');

figure(1)
subplot(4,2,6)
plot(100*diag(S)/sum(diag(S)),'ko','Linewidth',[2]), grid on
set(gca,'Xlim',[0 20],'Ylim',[0 80],'Ytick',[0 40 80],'Fontsize',[15])
xlabel('')

subplot(4,2,5)
plot(100*diag(S2)/sum(diag(S2)),'ko','Linewidth',[2]), grid on
set(gca,'Xlim',[0 20],'Ylim',[0 80],'Ytick',[0 40 80],'Fontsize',[15])
xlabel('')

figure(2)
subplot(3,1,3)
sn=[-1 1 1 1 -1];
for j=1:5
  nrm=sqrt(trapz(x,U(:,j).^2));
  Up(:,j)=real(U(:,j))*sn(j)/nrm;
end
plot(x,real(Up(:,1:5)),'Linewidth',[2])
set(gca,'Xlim',[-5 5],'Xtick',[-5 -2.5 0 2.5 5],'Ylim',[-1 1],'Ytick',[-1 -0.5 0 0.5 1],'Fontsize',[15])
legend('mode 1','mode 2','mode 3','mode 4','mode 5')

subplot(3,1,2)
sn=[-1 -1 1 -1 -1];
for j=1:5
  nrm=sqrt(trapz(x,U2(:,j).^2));
  Up2(:,j)=real(U2(:,j))*sn(j)/nrm;
end
plot(x,real(Up2(:,1:5)),'Linewidth',[2])
set(gca,'Xlim',[-5 5],'Xtick',[-5 -2.5 0 2.5 5],'Ylim',[-1 1],'Ytick',[-1 -0.5 0 0.5 1],'Fontsize',[15])
% subplot(4,1,4)
% sn=[-1 1 1 1 -1];
% for j=1:5
%   Up3(:,j)=real(U3(:,j))*sn(j)/norm(U3(:,j));
% end
% plot(x,real(Up3(:,1:5)),'Linewidth',[2])
% set(gca,'Xlim',[-5 5],'Xtick',[-5 -2.5 0 2.5 5],'Ylim',[-0.2 0.2],'Ytick',[-0.2 -0.1 0 0.1 0.2],'Fontsize',[15])

h=[1+0*x; 2*x; 4*x.^2-2; 8*x.^3-12*x; 16*x.^4-48*x.^2+12];
for j=0:4
   phi(:,j+1)=(1/(sqrt(factorial(j)*(2^j)*sqrt(pi))).*exp(-x.^2/2).*h(j+1,:)).';
   nrm=sqrt(trapz(x,phi(:,j+1).^2));
   phi2(:,j+1)=phi(:,j+1)/nrm;
end
subplot(3,1,1), plot(x,phi2,'Linewidth',[2]),
set(gca,'Xlim',[-5 5],'Xtick',[-5 -2.5 0 2.5 5],'Ylim',[-1 1],'Ytick',[-1 -0.5 0 0.5 1],'Fontsize',[15])

author: dynamicslab


text:
function epsil = precise_rand(n, delta, logm, shutup)
% This function computes (almost) precise random coding bound for BSC. The only imprecision is due to
% union bound.
% epsil <= sum_k  \delta^k (1-delta)^(n-k) P[one of M-1 other codewords is inside sphere <= k]

Ks = 0:n;
bino_coeffs = (gammaln(n+1) - gammaln(Ks+1) - gammaln(n-Ks+1))/log(2);
terms_spheres = bino_coeffs + logm - n;
[tmp log_spheres] = sumlog2(terms_spheres);

% cut off stupid >1 values (this what gallager rho-trick is supposed to do also)
log_spheres = min(log_spheres, 0);

terms = bino_coeffs + log2(delta)*Ks + log2(1-delta)*(n-Ks) + log_spheres;

epsil = 2^sumlog2(terms);

if(nargin < 4) || (isempty(shutup))
	disp(sprintf('-- precise_rand(n = %d, delta = %g, logm = %.1f): epsil = %.3g', n, delta, logm, epsil));
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function p = predict(Theta1, Theta2, X)
%PREDICT Predict the label of an input given a trained neural network
%   p = PREDICT(Theta1, Theta2, X) outputs the predicted label of X given the
%   trained weights of a neural network (Theta1, Theta2)

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned neural network. You should set p to a
%               vector containing labels between 1 to num_labels.
%
% Hint: The max function might come in useful. In particular, the max
%       function can also return the index of the max element, for more
%       information see 'help max'. If your examples are in rows, then, you
%       can use max(A, [], 2) to obtain the max for each row.
%

a1 = [ones(m,1) X];
z2 = a1 *Theta1';
a2 = [ones(size(z2),1) sigmoid(z2)];
z3 = a2*Theta2';
a3 = sigmoid(z3);

[predict_max, index_max] = max(a3, [], 2);

p = index_max;

% =========================================================================

end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function options = lipsolset(varargin)
%LIPSOLSET  Create or alter the options for Optimization with LIPSOL
%
% options = lipsolset('param1',value1,'param2',value2,...) creates an LIPSOL
% options structure with the parameters 'param' set to their corresponding
% values in 'value'. Parameters not specified will be set to the LIPSOL
% default.
%
% options = lipsolset(oldopts,'param1',value1,...) creates a copy of the old
% options 'oldopts' and then fills in (or writes over) the parameters
% specified by 'param' and 'value'.
%
% options = lipsolset() creates an options structure with all fields set to
% lipsolset defaults.
%
% lipsolset() prints a list of all possible fields and their function.
%
% See supplied LIPSOL Documentation for further details of these options.

%   Copyright (C) 2013 Jonathan Currie (I2C2)

% Print out possible values of properties.
if((nargin == 0) && (nargout == 0))
    printfields();
    return
end

%Names and Defaults
Names = {'phi0';'tau0';'monitor';'verb';'big';'cache_size';'loop_level'};
Defaults = {1e-5,0.9995,0,1,1e32,16,8}';

%Enter and check user args
try
    options = opticheckset(Names,Defaults,@checkfield,varargin{:});
catch ME
    throw(ME);
end

function checkfield(field,value)
%Check a field contains correct data type
switch lower(field)
    %Scalar double
    case 'big'
        err = opticheckval.checkScalarDbl(value,field);
    %Scalar non negative double
    case 'phi0'
        err = opticheckval.checkScalarNonNeg(value,field);
    %Scalar double 0 < val <= 1
    case 'tau0'
        err = opticheckval.checkScalarBoundLLE(value,field,0,1);
    %Scalar 0/1
    case 'monitor'
        err = opticheckval.checkScalar01(value,field);
    %Scalar Int 0-2
    case 'verb'
        err = opticheckval.checkScalarIntBoundLELE(value,field,0,2);
    %Scalar Int 1 - 100
    case 'cache_size'
        err = opticheckval.checkScalarIntBoundLELE(value,field,1,100);
    %Scalar Int 1-8
    case 'loop_level'
        err = opticheckval.checkScalarIntBoundLELE(value,field,1,8);
    otherwise
        err = MException('OPTI:SetFieldError','Unrecognized parameter name ''%s''.', field);
end
if(~isempty(err)), throw(err); end

function printfields()
%Print out fields with defaults
fprintf('       phi0: [ Relative Termination Tolerance: {1e-5} ] \n');
fprintf('       tau0: [ Initial Line-Search Step-Size Multiplier: {0.9995} ] \n');
fprintf('        big: [ Bound Infinite Value: {1e32} ] \n');
fprintf('    monitor: [ Graphical Iteration Monitor: On (1), Off {0} ] \n');
fprintf('       verb: [ Verbosity: On (1), Off {0} ] \n');
fprintf(' cache_size: [ SYMFCT Cache Size: {16} ] \n');
fprintf(' loop_level: [ BLKFCT Loop Unroll Level: 1,2,4,{8} ] \n');

fprintf('\n');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


function [l,c]=count(a);
%
% function [l,c]=count(a);
%
% INPUT:
%
% a : sorted list of doubles
%
% OUTPUT:
%
% l : sorted list of doubles with only unique elements
% c : count how many times each element appears
%
% EXAMPLE:
%
% [l,c]=count([1 1 2 2 3 3 3 4 5 5]);
%
% l=1 2 3 4 5
% c=2 2 3 1 2
%
%
% copyright 2005 by Kilian Q. Weinberger
% University of Pennsylvania
% kilianw@seas.upenn.edu
% ********************************************

error('not yet implemented!\n\n\n');
[l,index]=unique(a(1,:));
l=a(:,index);
c=[index(1) diff(index)];


text:
function [x_guess Iter] = VMPCM(ode,tau,x_guess,omega1,omega2,errTol,varargin)
%Purpose:
%Generic Function wrapper for the Vectorized Picard Chebyshev Method
%-------------------------------------------------------------------------%
%                                                                         %
% Inputs:                                                                 %
%--------
%ode                    object                          function name to
%                                                       evaluate i.e.
%                                                       @jatForces
%
%tau                    [N x 1]                         transformed time
%                                                       domain vector
%
%x_guess                [N x M]                         Initial Guess of
%                                                       solution values for
%                                                       the Picard
%                                                       Chebyshev Method
%
%omega1                 double                          First Omega Term
%
%omega2                 double                          Second Omega Term
%
%errTol                 double                          Error Tolerance of
%                                                       solution
%
%varargin                                                Additional inputs
%                                                       which are needed to
%                                                       be passed to the
%                                                       evaluation function
%                                                       ode
%
% Outputs:
%---------                                                                %
%x_guess                  [N x M]                        Refined solution
%                                                        meeting the error
%                                                        tolerances defined
%                                                        by errTol
%
%--------------------------------------------------------------------------
% Programmed by Darin Koblick 03-04-2012                                  %
%--------------------------------------------------------------------------
tau = tau(:)';
N = numel(tau)-1;
err1 = Inf;
err2 = Inf;
Iter = 0;
MaxIter = 300;
%disp('Running The Picard Chebyshev Algorithm')
tic;
%Initialize Constant Matricies up front so we don't have to recompute them
%which would slow down the iterative process
%--------------------------------------------------------------------------
T = ChebyshevPolynomial((0:N+1)',tau);
V = ones(1,length(tau))./N;
V(2:end-1) = V(2:end-1).*2;
TV1 = bsxfun(@times,T(1:N,:),V);
TV2 = bsxfun(@times,T(3:N+2,:),V);
TV = bsxfun(@rdivide,(TV1-TV2),(2.*(1:N))');
TV(end,:) = TV1(end,:)./(2*N)';
S = 2.*((-1).^((1:N)+1));
Cx = T(1:N+1,1:N+1)';
Cx(:,1) = Cx(:,1)./2;
%--------------------------------------------------------------------------
while (any(errTol < err1 ) || any(errTol < err2)) && Iter < MaxIter
    Iter = Iter + 1;
    %disp(['-------- Iteration # ',num2str(Iter),' -----------']);
    input = {omega2.*tau+omega1,x_guess,varargin{:}};
    F = ode(input{:}).*omega2;
    if size(F,2) == 1 || size(F,1) == 1
       F = F';
    end
    Beta_r = NaN(size(TV,1),size(F,2));
    Beta_k = NaN(size(TV,1) + 1, size(F,2));
    x_new = Beta_k;
    %Matrix Multiply
    for i=1:size(F,2)
        Beta_r(:,i) = TV*F(:,i);
        Beta_k(:,i) = [S*Beta_r(:,i) + 2.*x_guess(1,i); Beta_r(:,i)];
        x_new(:,i) = (Cx*Beta_k(:,i));
    end
    if any(size(x_new) ~= size(F))
       x_new = x_new';
    end
    err2 = err1;
    err1 = max(abs(x_new - x_guess),[],2);
    disp(['Max Error Is found to be: ', num2str(max(err1))]);
    x_guess = x_new;
end
end

function Tk = ChebyshevPolynomial(k,tau)
%The Chebyshev polynomial,T, corresponding to degree k
Tk = cos(bsxfun(@times,k,acos(tau)));
end


text:
function [ res ] = bellman_residual(Qt,Qtmp)
%BELLMAN_RESIDUAL Computes the Bellman residual
%
%   input -------------------------------------------------------
%
%       o Q: Q-value function at time t
%
%       o Qtmp: Q-value function at time t-1
%

vt   = compute_value_function(Qt);
vtmp =  compute_value_function(Qtmp);

res = (vt - vtmp).^2;
res = sum(res(:)) / size(vt(:),1);

end

function v = compute_value_function(nfq_Q)

u           = BuildActionList();
num_actions = 2;
nbSamples   = 100;
xs          = linspace(-1,1,nbSamples);
[Xs,Ys]     = meshgrid(xs,xs);
test        = [Xs(:),Ys(:)];
vs          = zeros(size(test,1),num_actions);

M           = size(test,1);

if iscell(nfq_Q)
    for i=1:num_actions
        vs(:,i) = nfq_Q{i}.f(test);
    end
else
    for i=1:num_actions
        test_ui = [test,repmat(u(i),M,1)];
        vs(:,i) = nfq_Q.f(test_ui);
    end
end

v          = min(vs,[],2);

end


text:
classdef Metric < handle
    properties (Access = public)
        flops = 0;
        memory = 0;
        outputSize = 0;
    end

    methods (Access = public)
        function obj = Metric(flops,memory,outputSize)
            if nargin > 0
                obj.flops = flops;
                obj.memory = memory;
                obj.outputSize = outputSize;
            end
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function H = boxper(sequence,isplot,boxnumber)
%
% 'boxper' estimate the hurst parameter of a given sequence with modified
%     periodogram method.
%
% Inputs:
%     sequence: the input sequence for estimate
%     isplot: whether display the plot. without a plot if isplot equal to 0
% Outputs:
%     H: the estimated hurst coeffeient of the input sequence

%  Author: Chu Chen
%  Version 1.0,  03/10/2008
%  chen-chu@163.com
%

if nargin == 1
    boxnumber = 50;
    isplot = 0;
end

if nargin == 2
    boxnumber = 50;
end

if boxnumber < 30 || boxnumber > 100
     error('The input argument boxnumber must be a integer between [30,100]');
end

n = length(sequence);
Xk = fft(sequence);
P_origin = abs(Xk).^2/(2*pi*n);
P = P_origin(1:floor(n/2)+1);

cut_min = ceil(0.001*n/2);
M = floor(logspace(log10(cut_min),log10(0.1*n-cut_min),boxnumber+1));
M = unique(M);
N = length(M)-1;

x = zeros(1,N);
y = zeros(1,N);
for i = 1:N
    m1 = M(i) + cut_min;
    m2 = M(i+1) + cut_min;
    x(i) = log10((pi * (m2 - m1))/(n));
    y(i) = log10(sum(P(m1:m2))/(m2-m1+1));
end

X = x;
Y = y;
p1 = polyfit(X,Y,1);
Yfit = polyval(p1,X);
H = (1-(Yfit(end)-Yfit(1))/(X(end)-X(1)))/2;

if isplot ~= 0
    figure,clf,hold on;
    plot(x,y,'b.');
    plot(X,Yfit,'r-','LineWidth',2);
    xlabel('Log10(Frequency)'),ylabel('Log10(Periodogram)'),title('Boxed Periodogram Method');
end


text:
function [ y, symm ] = cvx_s_banded( m, n, symm, lower, upper )

%CVX_S_BANDED (U,L)-banded matrices.

if nargin < 4,
    error( 'Bandwidth arguments missing.' );
end
if ~isnumeric( lower ) || length( lower ) ~= 1 || lower < 0 || lower ~= floor( lower ),
    error( 'Bandwidth arguments must be nonnegative integers.' );
elseif nargin < 5,
    upper = lower;
elseif ~isnumeric( upper ) || length( upper ) ~= 1 || upper < 0 || upper ~= floor( upper ),
    error( 'Bandwidth arguments must be nonnegative integers.' );
end

stflag = length(symm) == 2;
if stflag,
    toep = symm(2);
    symm = symm(1);
else
    toep = false;
end

if symm,
    lower = min( lower, upper );
    upper = 0;
end

c    = 0 : n - 1;
c    = c( ones( 1, m ), : );
r    = ( 0 : m - 1 )';
r    = r( :, ones( 1, n ) );
temp = r - c;
temp = temp <= lower & temp >= -upper;
r    = r( temp );
c    = c( temp );
nu   = length( r );

if toep,
    v = r - c;
    v = abs( v ) + max( v ) * ( v < 0 ) + 1;
    nu = max( v );
    toep = false;
else
    v = ( 1 : nu )';
end

if symm,
    tt = r ~= c;
    r = [ r ; c(tt) ];
    c = [ c ; r(tt) ];
    v = [ v ; v(tt) ];
    symm = false;
end

if stflag,
    symm = [ symm, toep ];
end

y = sparse( v, r + m * c + 1, 1, nu, m * n );

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% Compute the link density of a graph, defined as the number of edges divided by
% number_of_nodes(number_of_nodes-1)/2 where the latter is the maximum possible number of edges.
%
% Inputs: adjacency matrix, nxn
% Outputs: link density, a float between 0 and 1
%
% Note 1: The graph has to be non-trivial (more than 1 node).
% Note 2: Routine works for both directed and undirected graphs.
%
% Other routines used: numNodes.m, numEdges.m, isDirected.m
% GB: last update Sep 19, 2012

function d=linkDensity(adj)

n = numNodes(adj);

coeff = 2;
if isDirected(adj); coeff = 1; end

d = coeff*numEdges(adj)/(n*(n-1));

author: aeolianine

repo: octave-networks-toolbox


text:
function A_sym = symetrizeMatrix(A)
[i,j,k]=find(A);
W=sparse([i; j], [j; i], ones(size(k,1)*2,1));
A_help=sparse([i; j], [j; i], [k; k]);
[i,j,k]=find(A_help);
[i,j,kk]=find(W);
A_sym=sparse(i,j,(kk.^(-1)).*k); %Now Kantennr_sym is a symetric form of Kantennr

author: Sable


text:
function t=getInvHaarTransform(im,N)
h=haarmtx(N);
t=h'*im*h;