text:
function varargout = process_corr1n_time( varargin )
% PROCESS_COHERE1N_TIME: Compute the time-resolved coherence between all the pairs of signals, in one file.
%
% USAGE:   OutputFiles = process_corr1n_time('Run', sProcess, sInputA)
%                        process_corr1n_time('Test')

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Time-resolved correlation NxN';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 653;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    sProcess.isSeparator = 1;

    % === CONNECT INPUT
    sProcess = process_corr1n('DefineConnectOptions', sProcess, 1);
    % Time window
    sProcess.options.win.Comment = 'Estimation window length:';
    sProcess.options.win.Type    = 'value';
    sProcess.options.win.Value   = {.350, 'ms', []};
    % overlap
    sProcess.options.overlap.Comment = 'Sliding window overlap:';
    sProcess.options.overlap.Type    = 'value';
    sProcess.options.overlap.Value   = {50, '%', []};
    % === TITLE
    sProcess.options.label2.Comment = '
Estimator options:';
    sProcess.options.label2.Type    = 'label';
%     % === P-VALUE THRESHOLD
%     sProcess.options.pthresh.Comment = 'Metric significativity:     p<';
%     sProcess.options.pthresh.Type    = 'value';
%     sProcess.options.pthresh.Value   = {0.05,'',4};
    % === SCALAR PRODUCT
    sProcess.options.scalarprod.Comment    = 'Compute scalar product instead of correlation
(do not remove average of the signal)';
    sProcess.options.scalarprod.Type       = 'checkbox';
    sProcess.options.scalarprod.Value      = 0;
    % === OUTPUT MODE
    sProcess.options.label3.Comment = '
Output configuration:';
    sProcess.options.label3.Type    = 'label';
%     sProcess.options.outputmode.Comment = {'Save individual results (one file per input file)', 'Save average connectivity matrix (one file)'};
%     sProcess.options.outputmode.Type    = 'radio';
%     sProcess.options.outputmode.Value   = 1;
    % === OUTPUT COMMENT TAG
    sProcess.options.commenttag.Comment = 'Comment tag: ';
    sProcess.options.commenttag.Type    = 'text';
    sProcess.options.commenttag.Value   = '';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA) %#ok<DEFNU>
    % Initialize returned values
    OutputFiles = {};
    % Forcing the concatenation of the inputs
    sProcess.options.outputmode.Comment = {'Save individual results (one file per input file)', 'Save average connectivity matrix (one file)'};
    sProcess.options.outputmode.Type    = 'radio';
    sProcess.options.outputmode.Value   = 2;
    % Input options
    OPTIONS = process_corr1n('GetConnectOptions', sProcess, sInputA);
    if isempty(OPTIONS)
        return
    end
    CommentTag = sProcess.options.commenttag.Value;
    % Metric options
    OPTIONS.Method     = 'corr';
    OPTIONS.pThresh    = 0.05;
    OPTIONS.RemoveMean = ~sProcess.options.scalarprod.Value;
    OPTIONS.isSave     = 0;
    % Sliding time windows options
    WinLength  = sProcess.options.win.Value{1};
    WinOverlap = sProcess.options.overlap.Value{1};

    % Read time information
    TimeVector   = in_bst(sInputA(1).FileName, 'Time');
    sfreq        = round(1/(TimeVector(2) - TimeVector(1)));
    % Select input time window
    TimeVector = TimeVector((TimeVector >= OPTIONS.TimeWindow(1)) & (TimeVector <= OPTIONS.TimeWindow(2)));
    nTime      = length(TimeVector);

    % Compute sliding windows length
    Lwin  = round(WinLength * sfreq);
    Loverlap = round(Lwin * WinOverlap / 100);
    Nwin = floor((nTime - Loverlap) ./ (Lwin - Loverlap));
    % If window is bigger than the data
    if (Lwin > nTime)
        bst_report('Error', sProcess, sInputA, 'Sliding window for the coherence estimation is too long compared with the epochs in input.');
        return;
    end

    % Loop over all the time windows
    for iWin = 1:Nwin
        % Select time window
        iTimes = (1:Lwin) + (iWin-1)*(Lwin - Loverlap);
        OPTIONS.TimeWindow = TimeVector(iTimes([1,end]));
        % Compute metric
        ConnectMat = bst_connectivity({sInputA.FileName}, [], OPTIONS);
        % Processing errors
        if isempty(ConnectMat) || ~iscell(ConnectMat) || ~isstruct(ConnectMat{1}) || isempty(ConnectMat{1}.TF)
            bst_report('Error', sProcess, sInputA, 'Correction for the selected time segment could not be calculated.');
            return;
        end
        % Start a new brainstorm structure
        if (iWin == 1)
            NewMat = ConnectMat{1};
            NewMat.Time = OPTIONS.TimeWindow(1);
            NewMat.TimeBands = [];
        % Add next time point
        else
            NewMat.TF(:,iWin,:) = ConnectMat{1}.TF;
            NewMat.Time(end+1) = OPTIONS.TimeWindow(1);
        end
    end

    % Fix time vector
    if (length(NewMat.Time) == 1)
        bst_report('Warning', sProcess, sInputA, 'Only one sliding time window could be estimated.');
        NewMat.Time = [TimeVector(1), TimeVector(end)];
    end
    % Add comment tag
    if ~isempty(CommentTag)
        NewMat.Comment = [NewMat.Comment ' | ' CommentTag];
    end
    % Output filename
    sOutputStudy = bst_get('Study', OPTIONS.iOutputStudy);
    OutputFiles{1} = bst_process('GetNewFilename', bst_fileparts(sOutputStudy.FileName), 'timefreq_connectn_corr_time');
    % Save file
    bst_save(OutputFiles{1}, NewMat, 'v6');
    % Add file to database structure
    db_add_data(OPTIONS.iOutputStudy, OutputFiles{1}, NewMat);
end

author: brainstorm-tools


text:
function [ xtab_new, ytab_new, weight_new ] = ...
  triangle_rule_adjust ( xval, yval, norder, xtab, ytab, weight )

%*****************************************************************************80
%
%% TRIANGLE_RULE_ADJUST adjusts a unit quadrature rule to an arbitrary triangle.
%
%  Integration region:
%
%    Points (X,Y) such that:
%
%      (X,Y) = ALPHA * (X1,Y1) + BETA * (X2,Y2) + ( 1 - ALPHA - BETA ) * (X3,Y3)
%      0 <= ALPHA <= 1 - BETA
%      0 <= BETA <= 1 - ALPHA
%
%  Discussion:
%
%    This routine accepts as input abscissas and weights appropriate for
%    quadrature in a unit triangle, and returns abscissas and weights
%    appropriate for quadrature in a given triangle.
%
%    Once this routine has been called, an integral over the given triangle
%    can be approximated as:
%
%      QUAD = Sum ( 1 <= I <= NORDER ) WTAB2(I) * FUNC ( XTAB2(I), YTAB2(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XVAL(3), YVAL(3), the coordinates of the nodes.
%
%    Input, integer NORDER, the order of the rule.
%
%    Input, real XTAB(NORDER), YTAB(NORDER), the abscissas for
%    a unit triangle.
%
%    Input, real WEIGHT(NORDER), the weights for a unit triangle.
%
%    Output, real XTAB_NEW(NORDER), YTAB_NEW(NORDER), the adjusted
%    abscissas.
%
%    Output, real WEIGHT_NEW(NORDER), the adjusted weights.
%
  volume = triangle_volume ( xval, yval );

  xtab_new(1:norder) = ...
                  xtab(1:norder)                    * xval(1) ...
                                 + ytab(1:norder)   * xval(2) ...
    + ( 1.0E+00 - xtab(1:norder) - ytab(1:norder) ) * xval(3);

  ytab_new(1:norder) = ...
                  xtab(1:norder)                    * yval(1) ...
                                 + ytab(1:norder)   * yval(2) ...
    + ( 1.0E+00 - xtab(1:norder) - ytab(1:norder) ) * yval(3);

  weight_new(1:norder) = weight(1:norder) * 2.0E+00 * volume;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [regression_label] = fast_rcnn_bbox_transform(ex_boxes, gt_boxes)
% [regression_label] = fast_rcnn_bbox_transform(ex_boxes, gt_boxes)
% --------------------------------------------------------
% Fast R-CNN
% Reimplementation based on Python Fast R-CNN (https://github.com/rbgirshick/fast-rcnn)
% Copyright (c) 2015, Shaoqing Ren
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

    ex_widths = ex_boxes(:, 3) - ex_boxes(:, 1) + 1;
    ex_heights = ex_boxes(:, 4) - ex_boxes(:, 2) + 1;
    ex_ctr_x = ex_boxes(:, 1) + 0.5 * (ex_widths - 1);
    ex_ctr_y = ex_boxes(:, 2) + 0.5 * (ex_heights - 1);

    gt_widths = gt_boxes(:, 3) - gt_boxes(:, 1) + 1;
    gt_heights = gt_boxes(:, 4) - gt_boxes(:, 2) + 1;
    gt_ctr_x = gt_boxes(:, 1) + 0.5 * (gt_widths - 1);
    gt_ctr_y = gt_boxes(:, 2) + 0.5 * (gt_heights - 1);

    targets_dx = (gt_ctr_x - ex_ctr_x) ./ (ex_widths+eps);
    targets_dy = (gt_ctr_y - ex_ctr_y) ./ (ex_heights+eps);
    targets_dw = log(gt_widths ./ ex_widths);
    targets_dh = log(gt_heights ./ ex_heights);

    regression_label = [targets_dx, targets_dy, targets_dw, targets_dh];
end

author: ShaoqingRen


text:
function r=transpose(r)
% transpose array of rotations

r = transpose@quaternion(r);
r.i = r.i.';

author: mtex-toolbox

repo: mtex

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


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
% ##2
%==============================================================================
% function [Tc,dT] = linearInter(T,omega,x,varargin)
%
% linear interpolator for the data T given on a cell-centered grid
% to be evaluated at locations x
%
% Input
%  T        coefficients for image model img : \R^{dim*n} \to\R^{n)},
%           for j=1:n, img(j)  = sum T_i b_i(x_j); end;
%  omega    specification of domain
%           Omega = (omega(1),omega(2)) x ... x (omega(2*dim-1,omega(2*dim)),
%           spatial dimension dim = length(omega/2)
%  x        location evaluation points
%           x = reshape(x,n,dim), x(j,:) coordinates of j-th point
%  varargin additional parameters such as doDerivative,
%##4 cleanup
%==============================================================================

function [Tc,dT] = linearInter(T,omega,x,varargin)

Tc = mfilename('fullpath'); dT = [];

if nargin == 0,
  help(mfilename);
  testOneImgModel(mfilename);
  return;
elseif nargin == 1 && isempty(T),
  return;
end;

% flag for computing the derivative
doDerivative = (nargout>1);
matrixFree   = 0;
for k=1:2:length(varargin), % overwrite default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

%------------------------------------------------------------------------------

% get data size m, cell size h, dimension d, and number n of interpolation points
dim = length(omega)/2;
m   = size(T);         if dim == 1, m = numel(T); end;
h   = (omega(2:2:end)-omega(1:2:end))./m;
n   = length(x)/dim;
x   = reshape(x,n,dim);

% map x from [h/2,omega-h/2] -> [1,m],
for i=1:dim, x(:,i) = (x(:,i)-omega(2*i-1))/h(i) + 0.5; end;

Tc = zeros(n,1); dT = [];                     % initialize output
if doDerivative, dT = zeros(n,dim);  end;     % allocate memory in column format
Valid = @(j) (0<x(:,j) & x(:,j)<m(j)+1);      % determine indices of valid points

switch dim,
  case 1, valid = find( Valid(1) );
  case 2, valid = find( Valid(1) & Valid(2) );
  case 3, valid = find( Valid(1) & Valid(2) & Valid(3) );
end;

if isempty(valid),
  if doDerivative, dT = sparse(n,dim*n); end; % allocate memory in column format
  return;
end;

pad = 1; TP = zeros(m+2*pad);                 % pad data to reduce cases

P = floor(x); x = x-P;                        % split x into integer/remainder
p = @(j) P(valid,j); xi = @(j) x(valid,j);

% increments for linearized ordering
i1 = 1; i2 = size(T,1)+2*pad; i3 = (size(T,1)+2*pad)*(size(T,2)+2*pad);

switch dim,
  case 1,
    TP(pad+(1:m)) = reshape(T,m,1);
    clear T;
    p = pad + p(1);
    Tc(valid) = TP(p).* (1-xi(1)) + TP(p+1).*xi(1);   % compute weighted sum

    if ~doDerivative, return; end;
    % compute and format the derivative
    dT(valid) = TP(p+1)-TP(p);
  case 2,
    TP(pad+(1:m(1)),pad+(1:m(2))) = T;
    clear T;
    p  = (pad + p(1)) + i2*(pad + p(2) - 1);
    % compute Tc as weighted sum
    Tc(valid) = (TP(p)   .* (1-xi(1)) + TP(p+i1)    .*xi(1)) .* (1-xi(2)) ...
      + (TP(p+i2) .* (1-xi(1)) + TP(p+i1+i2) .*xi(1)) .* (xi(2));

    if ~doDerivative, return; end;
    dT(valid,1) = (TP(p+i1)-TP(p)).*(1-xi(2)) + (TP(p+i1+i2)-TP(p+i2)).*xi(2);
    dT(valid,2) = (TP(p+i2)-TP(p)).*(1-xi(1)) + (TP(p+i1+i2)-TP(p+i1)).*xi(1);
  case 3,
    TP(pad+(1:m(1)),pad+(1:m(2)),pad+(1:m(3))) = T;
    clear T;
    p  = (pad + p(1)) + i2*(pad + p(2) - 1) + i3*(pad + p(3) -1);
    % compute Tc as weighted sum
    Tc(valid) = ((TP(p).*(1-xi(1))+TP(p+i1).*xi(1)).*(1-xi(2))...
      +(TP(p+i2).*(1-xi(1))+TP(p+i1+i2).*xi(1)).*(xi(2))).*(1-xi(3)) ...
      +((TP(p+i3).*(1-xi(1))+TP(p+i1+i3).*xi(1)).*(1-xi(2)) ...
      +(TP(p+i2+i3).*(1-xi(1))+TP(p+i1+i2+i3).*xi(1)).*(xi(2))).*(xi(3));

    if ~doDerivative, return; end;
    dT(valid,1) = ((TP(p+i1)-TP(p)).*(1-xi(2))+(TP(p+i1+i2)-TP(p+i2)).*xi(2)).*(1-xi(3)) ...
      +((TP(p+i1+i3)-TP(p+i3)).*(1-xi(2))+(TP(p+i1+i2+i3)-TP(p+i2+i3)).*xi(2)).*(xi(3));
    dT(valid,2) = ((TP(p+i2)-TP(p)).*(1-xi(1))+(TP(p+i1+i2)-TP(p+i1)).*xi(1)).*(1-xi(3)) ...
      +((TP(p+i2+i3)-TP(p+i3)).*(1-xi(1))+(TP(p+i1+i2+i3)-TP(p+i1+i3)).*xi(1)).*(xi(3));
    dT(valid,3) = ((TP(p+i3).*(1-xi(1))+TP(p+i1+i3).*xi(1)).*(1-xi(2)) ...
      +(TP(p+i2+i3).*(1-xi(1))+TP(p+i1+i2+i3).*xi(1)).*(xi(2))) ....
      -((TP(p).*(1-xi(1))+TP(p+i1).*xi(1)).*(1-xi(2)) ...
      +(TP(p+i2).*(1-xi(1))+TP(p+i1+i2).*xi(1)).*(xi(2)));
end;
if doDerivative
    for i=1:dim, dT(:,i) = dT(:,i)/h(i); end
    if not(matrixFree)
        dT = spdiags(dT,n*(0:(dim-1)),n,dim*n);
    end
end
%==============================================================================


text:
function x0 = secondary_source_tapering(x0,conf)
%SECONDARY_SOURCE_TAPERING applies a tapering window to the secondary sources
%
%   Usage: x0 = secondary_source_tapering(x0,conf)
%
%   Input parameters:
%       x0          - secondary sources / m
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       x0          - secondary sources / m, containing the applied tapering
%                     window in its weights in x0(:,7)
%
%   SECONDARY_SOURCE_TAPERING(x0,conf) applies a tapering window (depending on
%   the conf.usetapwin and conf.tapwinlen settings) to the secondary sources.
%   It is applied to the weights stored in x0(:,7).
%
%   See also: secondary_source_positions, tapering_window

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 2;
nargmax = 2;
narginchk(nargmin,nargmax);

%% ===== Calculation ====================================================
% Apply tapering window to secondary sources
x0(:,7) = x0(:,7) .* tapering_window(x0,conf);

author: sfstoolbox

repo: sfs-matlab


text:
function test_suite=test_montecarlo_cluster_stat_distribution
% probability uniformity tests for cosmo_montecarlo_cluster_stat
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_mccs_uniformity_slow()
    show_progress=true;

    if show_progress
        cleaner=onCleanup(@()progress_helper('-clear'));
        progress_helper(sprintf('Slow test in %s ',...
                                    mfilename()));
    end

    % test for uniformity of p-values of monte_carlo_cluster_stat
    %
    % this test aims to verify that there is not an excessive number of
    % false positives under the null hypothesis. It does so by running
    % monte_carlo_cluster_stat several times on random data, storing the p
    % values for each iteration, and comparing those to a uniform
    % distribution.
    %
    % because running this test is slow, it can perform the same test
    % multiple times, each time increasing the number of iterations. This
    % is repeated until enough evidence is gathered that p values are
    % uniform or not.

    max_attempts=8;
    grow_niter=1.2;

    % benchmark values obtained by runnning
    % helper_mccs_get_correlation_with_uniform multiple times with
    % 'correct' monte carlo custer stat function.
    uniform_c_mu=.985;
    uniform_c_sd=.05;

    % the null hypothesis is that p values are uniformly distributed;
    % earlier versions of monte_carlo_cluster_stat would fail this test,
    % with a lower correlation value as a result. The correlation value
    % is converted to a z-score
    %
    pass_min_z=-1;
    fail_max_z=-5;

    min_pass_or_fail_count=2;
    count_pass=0;
    count_fail=0;

    ps_cell=cell(1,max_attempts);

    niter=20;
    for attempt=1:max_attempts
        ps_cell{attempt}=helper_mccs_get_pvalues(niter,show_progress);
        ps=sort(cat(1,ps_cell{:}));

        % compute correlation with expected p-values
        n_ps=numel(ps);
        ps_uniform=(.5:n_ps)'/n_ps;
        c=cosmo_corr(ps,ps_uniform);

        z=sqrt(niter)*(c-uniform_c_mu)/uniform_c_sd;

        if z>pass_min_z
            count_pass=count_pass+1;
            count_fail=0;

        elseif z<fail_max_z
            count_fail=count_fail+1;
            count_pass=0;
        else
            count_fail=0;
            count_pass=0;
        end

        if count_pass>=min_pass_or_fail_count
            finalize_test_helper(show_progress);

            % test passes
            return;

        elseif count_fail>=min_pass_or_fail_count
            % enough evidence that p values are non-uniform, fail
            finalize_test_helper(show_progress);

            error(['Found z=%d, indicating that probability values '...
                        'are probably not uniform'],z);
        end

        % not enough evidence for either uniform or non-uniform, redo the
        % test with more iterations
        niter=ceil(niter*grow_niter);
        progress_helper('#');

    end

    finalize_test_helper(show_progress);
    error('Maximum number of attempts reached');

function finalize_test_helper(show_progress)
    if show_progress
        fprintf('\n');
    end

function ps=helper_mccs_get_pvalues(niter,show_progress)
    % output: correlation between expected uniform distribution of p values
    % and those obtained from monte_carl_cluster_stat

    niter_tfce=25;
    nsubj=10;

    ps=zeros(niter,1);

    % dataset with single features
    ds=struct();
    ds.samples=randn(nsubj,1);
    ds.sa.targets=ones(nsubj,1);
    ds.sa.chunks=(1:nsubj)';
    ds.fa=struct();
    ds.a=struct();

    % trivial (singleton) neighborhood
    nh=struct();
    nh.neighbors={1};
    nh.fa=struct();
    nh.fa.sizes=1;
    nh.a=struct();
    nh.origin.fa=ds.fa;
    nh.origin.a=ds.a;

    for iter=1:niter
        opt=struct();
        opt.niter=niter_tfce;
        opt.progress=false;
        opt.h0_mean=0;
        opt.dh=.1;

        ds.samples=randn(nsubj,1);

        z=cosmo_montecarlo_cluster_stat(ds,nh,opt);

        ps(iter)=normcdf(z.samples);

        if show_progress
            progress_helper(':');
        end
    end

function progress_helper(what)
    % helper to show progress during the test, which then flushes at the
    % end
    persistent delete_count;

    if isempty(delete_count)
        delete_count=0;
    end

    if strcmp(what,'-clear')
        to_print=repmat(sprintf('\b'),1,delete_count);
        delete_count=0;
    else
        to_print=what;
        delete_count=delete_count+numel(what);
    end

    fprintf(to_print);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
%trotz SE(3) rotation about Z axis
%
% T = trotz(THETA) is a homogeneous transformation (4x4) representing a rotation
% of THETA radians about the z-axis.
%
% T = trotz(THETA, 'deg') as above but THETA is in degrees.
%
% Notes::
% - Translational component is zero.
%
% See also rotz, trotx, troty, trot2, SE3.Rz.

%## 3d homogeneous rotation

% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function T = trotz(t, varargin)
	T =    [rotz(t, varargin{:}) [0 0 0]'; 0 0 0 1];

author: petercorke

repo: spatialmath-matlab


text:
% ----------------------------------------------------------------------------
% function hfssCreateRelativeCS(fid, Name, Origin, Units)
%
% Description :
% -------------
% Create a relative Coordinate System at "origin" point relative to the
% Global Coordinate System.
%
% Parameters :
% ------------
% fid     - file identifier of the HFSS script file.
% Name    -
% Origin  - point relative to the Global Coordinate System where the
%           Relative Coordinate System will be created.
% Units   - units of the points (specify using either 'in', 'mm', 'meter'
%           or anything else defined in HFSS).
%
% Note :
% ------
% Use this function to facilitate drawing stuff. To return to use
% another CS (i.e. the Global CS), you can use the function hfssSetWCS
%
% Example :
% ---------
% fid = fopen('myantenna.vbs', 'wt');
% ...
% hfssInsertDesign(fid, 'Dipole_SingleElement');
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% Written by Daniel R. Prado
% danysan@gmail.com / drprado@tsc.uniovi.es
% 20 September 2012
% ----------------------------------------------------------------------------

function hfssCreateRelativeCS(fid, Name, Origin, Units)

% arguments processor.
if (nargin < 4)
	error('Insufficient number of arguments !');
end

% Preamble.
fprintf(fid, '\n');
fprintf(fid, 'oEditor.CreateRelativeCS _\n');

% CS Parameters
fprintf(fid, 'Array("NAME:RelativeCSParameters", _\n');
fprintf(fid, '"OriginX:=", "%.4f%s", _\n', Origin(1), Units);
fprintf(fid, '"OriginY:=", "%.4f%s", _\n', Origin(2), Units);
fprintf(fid, '"OriginZ:=", "%.4f%s", _\n', Origin(3), Units);
fprintf(fid, '"XAxisXvec:=", "1%s", _\n', Units);
fprintf(fid, '"XAxisYvec:=", "0%s", _\n', Units);
fprintf(fid, '"XAxisZvec:=", "0%s", _\n', Units);
fprintf(fid, '"YAxisXvec:=", "0%s", _\n', Units);
fprintf(fid, '"YAxisYvec:=", "1%s", _\n', Units);
fprintf(fid, '"YAxisZvec:=", "0%s"), _\n', Units);

% CS Attributes
fprintf(fid, 'Array("NAME:Attributes", _\n');
fprintf(fid, '"Name:=", "%s")\n', Name);

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function gbasis=groebner(polyset,ord,varnames,tol)
% GROEBNER - calculate the reduced Groebner basis of a set of polynomials
%  usage: gbasis = groebner(polyset,ord)
%         gbasis = groebner(polyset,ord,varnames)
%         gbasis = groebner(polyset,ord,varnames,tol)
%
% INPUTS: polyset, ord, varnames, tol
%  polyset is a cell array of multivariate polynomial coefficients
%    the cell elements can be all strings or all arrays.
%    If the polynomials are specified with arrays, each row of the array
%    represents a single term.  The first element is the coefficient and
%    the others are the degrees of each unknown in the term.
%    e.g. [-1,2,1;1,0,1] represents -x1^2*x2+x2
%    If the polynomials are specified as strings, they must use variable
%    names 'x1','x2',... (unless varnames is provided) and be valid inputs
%    for str2poly.m
%  ord is a string specifying the monomial ordering:
%    'lex': lexicographical, order by highest power of most significant
%      indeterminate, e.g. x1>x2^2, x1*x2>x1.
%    'grlex': graded lex, order by total degree then lex,
%      e.g. x1^2*x2 > x1*x2^2 > x1^2
%    'grevlex': graded reverse lex, order by total degree then by lowest
%      power of least significant indeterminate, e.g. x1*x2^2 > x1*x2*x3
%    (for all cases have x1>x2>...>xn)
%    note: revlex is not a well-ordering because 1>x1>x1^2>...
%  varnames (optional) is a cell array of variable names if polyset is a
%    cell array of strings and the variable names are not 'x1', 'x2', etc.
%    If specified, there must always be at least 2 variables.
%  tol (optional) is the zero tolerance for coefficients, otherwise
%    catastrophic cancellation may occur.  Default value is 0
%
% OUTPUTS: gbasis
%  gbasis is the reduced Groebner basis of the set of polynomials, in the
%    same format (strings or matrices) as polyset{1}.
%
% ALGORITHM:
%  A modified Buchberger's algorithm is used.  In Buchberger's algorithm,
%  at each step calculate Sij = (Lij/ai)*Pi - (Lij/aj)*Pj, for all pairs of
%  polynomials Pi, Pj, where ai and aj are the leading terms of Pi and Pj,
%  and Lij is the least common multiple of ai and aj; then reduce Sij with
%  respect to the polynomials and add it to the set (note that Sij always
%  reduces to zero if ai and aj have no variables in common).  The
%  algorithm terminates when no new polynomials can be added to the set.
%  This function reduces the polynomials at each step in order to find the
%  reduced Groebner basis.
%
% EXAMPLE: simplify equations x^2+2xy^2=0, xy+2y^3=1
%  A1=[1,2,0;2,1,2]; % 1*x^2*y^0 + 2*x^1*y^2
%  A2=[-1,0,0;1,1,1;2,0,3]; % -1*x^0*y^0 + 1*x^1*y^1 + 2*x^0*y^3
%  y=groebner({A1,A2},'lex');
%  poly2str(y) % returns {'x2^3-0.5', 'x1'}
%
% EXAMPLE: same, but with equations specified as strings:
%  groebner({'x^2+2*x*y^2','x*y+2*y^3-1'},'lex',{'x','y'})
%  % returns {'y^3-0.5','x'}
%
% EXAMPLE: show that a set of equations is inconsistent
%  groebner({'x + y', 'x^2 - 1', 'y^2 - 2*x'}, 'lex', {'x', 'y'})
%  % returns {'1'}
%
%
% KNOWN ISSUES
%
%  0. GENERAL: while a significant effort has been made to ensure that the
%    program works as described, unidentified bugs may remain and the results
%    should be regarded with caution.  In particular, the suite of tests
%    performed was by no means comprehensive.
%
%  1. ACCURACY: the algorithm by default works in floating-point which
%    ultimately cannot perform exact arithmetic, so unexpected silent errors may
%    occur even when the coefficients are specified as integers.  A workaround
%    is to use other data types but this has its own issues (see below).
%
%  1a. EXAMPLE (roundoff - linearly independent polynomials)
%    eps % returns 2.2204e-16 (32-bit PC windows machine epsilon)
%    groebner({'x1+x2','x1+1.000000000000001*x2'},'lex') % ok, returns {'x2','x1'}
%    groebner({'x1+x2','x1+1.0000000000000001*x2'},'lex') % not ok, returns {'x2+x1'}
%
%  1b. EXAMPLE (roundoff - linearly dependent polynomials)
%    2/sqrt(2)-sqrt(2) % expect 0, returns -2.2204e-016
%    poly2str(groebner({[1,1,0;sqrt(2),0,1],[1/sqrt(2),1,0;1,0,1]},'lex')) % ok, returns {'1.41421*x2+x1'}
%    poly2str(groebner({[1,1,0;sqrt(2),0,1],[sqrt(2),1,0;2,0,1]},'lex')) % not ok, returns {'x2','x1'}
%    poly2str(groebner({[1,1,0;sqrt(2),0,1],[sqrt(2),1,0;2,0,1]},'lex',{},1e-14)) % ok, returns {'1.41421*x2+x1'}
%
%  1c. EXAMPLE (courtesy of Christophe Lauwerys - catastrophic cancellation & tolerance bug)
%    groebner({'t^3+x+y','t^2+0.5*x^2-x-z^2','t^2+y-z^2'},'lex',{'t','x','y','z'}) % ans{1} contains Inf
%
%  1d. EXAMPLE (from Mathematica documentation - catastrophic cancellation because of large intermediate coefficients)
%    % expect 3 polynomials returns, one of degree 8 in z; instead can get '1'
%    groebner({'x^2 + y^2 + z^2 - 1', 'x*y - z + 2', 'z^2 - 2*x + 3*y'}, 'lex', {'x', 'y','z'})
%
%  2. NON-DOUBLE INPUTS: use of non-double coefficients is not supported.  While
%    it may be tempting to introduce symbolic coefficients, it is not clear how
%    they should be consistently treated - for example 'A*x' could reduce to '0'
%    or 'x' or both, depending upon the assumptions on the coefficient A.
%
%  3. EFFICIENCY: The reduction algorithms aren't terribly efficient.  One of
%    the major bottlenecks at present is the calculation of the leading term,
%    which currently involves a call to sortrows each time.  It could be more
%    efficient to maintain an ordered list or even a priority queue (e.g.
%    implemented as a heap).  Also lattice reduction algorithms such as LLL
%    might speed things up.
%
%
% SEE ALSO:
%  poly2str, str2poly, polynsolve

% Author: Ben Petschel 19/6/2009
%
% Change history:
%  19/6/2009 - first release (ord='grlex' and ord='grevlex' not
%    implemented yet)
%  22/6/2009 - implemented ord='grlex' and ord='grevlex'
%    reduction algorithm now reduces lower-order terms as well
%  23/6/2009 - removed ord='revlex' because it is not a well-ordering
%  17/7/2009 - fixed bug in handling polynomials with >=3 variables
%  20/3/2010 - changed polynomial representation from multidimensional
%    arrays to 2d arrays (more memory efficient with large number of vars)
%  11/10/2010 - bugfix for handling error tolerance

if (numel(polyset)>0) && ischar(polyset{1}),
  charout = true;
  if nargin<3,
    polyset = str2poly(polyset);
  else
    polyset = str2poly(polyset,varnames);
  end;
else
  charout = false;
end;

if nargin<4,
  tol = 0;
end;

% determine the number of degrees
nd=max(cellfun(@(x)size(x,2)-1,polyset)); % error if = 0
if nd == 0
  error('number of variables must be >= 1');
end;

gbasis = fullreduce(polyset,ord,tol,nd);
oldgbasis = {};

while ~isequal(oldgbasis,gbasis),
  oldgbasis = gbasis;
  Qset = SPoly(gbasis,ord,tol,nd);
  gbasis = fullreduce([gbasis,Qset],ord,tol,nd);
end;

if charout,
  if nargin<3,
    gbasis = poly2str(gbasis);
  else
    gbasis = poly2str(gbasis,varnames);
  end;
end;

end % main function groebner(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Qset = SPoly(Pset,ord,tol,nd)
% calculates all reduced S-polynomials of Pset (Pset must be reduced)
%   Sij = (Lij/ai)*Pi - (Lij/aj)*Pj,
% for all pairs of polynomials Pi, Pj, where ai and aj are the leading
% terms of Pi and Pj, and Lij is the least common multiple of ai and aj
% then reduces Sij wrt Pset

Qset = {};
for i=1:numel(Pset)-1,
  for j=i+1:numel(Pset),
    [c1,d1] = leadterm(Pset{i},ord,tol,nd);
    [c2,d2] = leadterm(Pset{j},ord,tol,nd);
    if any(and(d1>0,d2>0)),
      % leading terms have a variable in common, so S-poly is nontrivial
      L = max(d1,d2); % LCM of leading terms
      S = addpolys(multiplyterm(Pset{i},L-d1),multconst(multiplyterm(Pset{j},L-d2),-1));
      S = reduceset(S,Pset,ord,tol,nd);
      %if any(S(:)~=0),
      if abs(leadterm(S,ord,tol,nd))>tol
        Qset{end+1}=S; % add S to the list
      end;
    end;
  end;
end;

end % helper function SPoly(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Qset = fullreduce(Pset,ord,tol,nd)
% reduces all polynomials wrt each other
Qset = Pset;
oldQset = {};

while ~isequal(oldQset,Qset),
  oldQset = Qset;
  for i=1:numel(Qset)-1,
    for j=i+1:numel(Qset),
      Qset{i}=reduce(Qset{i},Qset{j},ord,tol,nd);
      Qset{j}=reduce(Qset{j},Qset{i},ord,tol,nd);
    end;
  end;
end;

% keep only the nonzero results (to within zero tolerance tol)
Qset = Qset(cellfun(@(x)abs(leadterm(x,ord,tol,nd))>tol,Qset));

if numel(Qset)==1,
  % special case, a single equation can slip through without being reduced
  c = leadterm(Qset{1},ord,tol,nd);
  % coefficients occupy first column of the array
  Qset{1} = multconst(Qset{1},1/c); % know c~=0 because have already removed zero eqns
end;

end % helper function fullreduce(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q = reduceset(P,Pset,ord,tol,nd)
% reduces P wrt all polynomials in Pset
Q = P;
oldQ = [];

while ~isequal(oldQ,Q),
  oldQ = Q;
  for i=1:numel(Pset),
    Q=reduce(Q,Pset{i},ord,tol,nd);
  end;
end;

end % helper function reduceset(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q1=reduce(P1,P2,ord,tol,nd)
% reduces a polynomial P1 by subtracting multiples of P2 so that leading
% term coefficient is 1 and leading term of P2 does not divide that of P1

Q2=P2;
[c2,d2] = leadterm(Q2,ord,tol,nd);
if abs(c2)<=tol,
  % ignore terms less than tol
  keepgoing = false;
  Q1 = P1;
else
  keepgoing = true;
  Q2 = multconst(Q2,1/c2);
  Q1 = 0; % will successively add terms to Q
  % remove the leading terms, in order to reduce wrt lower terms
  remain = P1;
end;

while keepgoing,
  [c1,d1] = leadterm(remain,ord,tol,nd);
  if abs(c1)<=tol,
    % ignore terms less than tol
    keepgoing = false;
  else
    %reduce remain by Q2, if possible, or remove leading term
    if d1>=d2,
      % can reduce Q1 by Q2, so subtract a multiple of Q2 from Q1
      remain = addpolys(remain,multconst(multiplyterm(Q2,d1-d2),-c1));
    else
      % cannot reduce any further wrt leading term, but try lower terms
      lead = multiplyterm(c1,d1);
      remain = addpolys(remain,multconst(lead,-1));
      Q1 = addpolys(Q1,lead); % add irreducible terms to Q1
    end;
  end;
end;

% reduce leading coefficient, if possible
c1 = leadterm(Q1,ord,tol,nd);
if abs(c1)>tol,
  Q1 = multconst(Q1,1/c1);
end;

end % helper function reduce(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=multconst(P,c)
% returns c*P where c is a constant

Q = P;
Q(:,1)=Q(:,1)*c;

end % helper function multconst(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=addpolys(P1,P2)
% adds two polynomials

s1 = size(P1,2);
s2 = size(P2,2);

if s1>s2
  Q = [P1;P2,zeros(size(P2,1),s1-s2)];
else
  Q = [P1,zeros(size(P1,1),s2-s1);P2];
end;

Q = sortrows(Q,2:size(Q,2));
i=1;
while i<size(Q,1)
  % merge rows representing terms with the same exponents
  if Q(i,1)==0,
    Q(i,:)=[];
  elseif all(Q(i,2:end)==Q(i+1,2:end))
    Q(i,1)=Q(i,1)+Q(i+1,1);
    Q(i+1,:)=[];
    if Q(i,1)==0
      % remove terms that add to zero
      Q(i,:)=[];
    end;
  else
    i=i+1;
  end;
end;

end % helper function addpolys

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=multiplyterm(P,d)
% multiplies a polynomial by a term (no coefficient)

% padding a matrix at the end is same as multiplying by x
d=d(:).';
sP = size(P,2);
sd = size(d,2)+1; % d doesn't include coefficient
if sP<sd,
  Q=[P,zeros(size(P,1),sd-sP)];
else
  Q=P;
  d=[(d(:)).',zeros(1,sP-sd)];
end;
Q=Q+repmat([0,d],size(Q,1),1);

end % helper function multiplyterm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [coeff,deg]=leadterm(P,ord,tol,nd)
% returns the leading term of P
% deg is a row vector of length size(P,2)-1
% coeff is the leading coefficient, deg(i) is the degree of xi
% ignores terms with coefficients <= tol
% nd is the minimum number of degrees (in order to handle cases where deg=0)

d=size(P,2);

if all(abs(P(:,1))<=tol),
  % P=0
  coeff=0;
  deg=zeros(1,d-1);
else
  switch ord
    case 'lex'
      P=sortrows(P,2:d);

    case 'revlex'
      error('revlex is not a well-ordering');

    case 'grlex'
      P = [P,sum(P(:,2:d),2)];
      P = sortrows(P,[d+1,2:d]);
      P = P(:,1:d);

    case 'grevlex'
      P = [P,sum(P(:,2:d),2)];
      P = sortrows(P,[d+1,-(d:-1:2)]);
      P = P(:,1:d);

    otherwise
      error('ordering %s not recognized',ord);

  end; % switch ord

  ind = find(abs(P(:,1))>tol,1,'last');
  if isempty(ind),
    coeff = 0;
    deg = zeros(1,d-1);
  else
    coeff = P(ind,1);
    deg = P(ind,2:end);
  end;

end; % if all(P(:)==0)

if nd>d-1,
  deg = [deg,zeros(1,nd-d+1)];
end;

end % helper function leadterm(...)

author: Sable

repo: mcbench-benchmarks


text:
function []= fill_legend(obj,title_text,legend_text,legend_type,legend_color,legend_marker,legend_linestyle,legend_line_size,legend_point_size)

legend_y_step=1;
legend_y_additional_step=0.5;

obj.legend_text_handles=[obj.legend_text_handles...
    text(1,obj.legend_y,title_text,...
    'FontWeight','bold',...
    'Interpreter',obj.text_options.interpreter,...
    'FontName',obj.text_options.font,...
    'FontSize',obj.text_options.base_size*obj.text_options.legend_title_scaling,...
    'Parent',obj.legend_axe_handle)];

obj.legend_y=obj.legend_y-legend_y_step;

for ind_legend=1:length(legend_text)

    if size(legend_color,1)>1
        temp_color = legend_color(ind_legend,:);
    else
        temp_color = legend_color;
    end

    if length(legend_marker)>1
        temp_marker = legend_marker{ind_legend};
    else
        temp_marker = legend_marker;
    end

    if length(legend_linestyle)>1
        temp_linestyle = legend_linestyle{ind_legend};
    else
        temp_linestyle = legend_linestyle;
    end

   if length(legend_line_size)>1
        temp_line_size = legend_line_size(ind_legend);
    else
        temp_line_size = legend_line_size;
   end

   if length(legend_point_size)>1
        temp_point_size = legend_point_size(ind_legend);
    else
        temp_point_size = legend_point_size;
    end

    if strcmp(legend_type,'point')
        plot(1.5,obj.legend_y,temp_marker,...
            'MarkerEdgeColor','none',...
            'MarkerFaceColor', temp_color,...
            'Parent',obj.legend_axe_handle,...
            'MarkerSize',temp_point_size);
    end
    if strcmp(legend_type,'line') || strcmp(legend_type,'both')
        plot([1 2],[obj.legend_y obj.legend_y],temp_linestyle,...
            'Color',temp_color,...
            'lineWidth',temp_line_size,...
            'Parent',obj.legend_axe_handle)
    end
    if strcmp(legend_type,'both')
        plot(1,obj.legend_y,temp_marker,...
            'MarkerEdgeColor','none',...
            'MarkerFaceColor', temp_color,...
            'Parent',obj.legend_axe_handle,...
            'MarkerSize',temp_point_size);
    end

    obj.legend_text_handles=[obj.legend_text_handles...
        text(2.5,obj.legend_y,legend_text{ind_legend},...
        'Interpreter',obj.text_options.interpreter,...
        'FontName',obj.text_options.font,...
        'FontSize',obj.text_options.base_size*obj.text_options.legend_scaling,...
        'Parent',obj.legend_axe_handle)];

    obj.legend_y=obj.legend_y-legend_y_step;
end

obj.legend_y=obj.legend_y-legend_y_additional_step;

end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
% This script will pass sample parameters to gradient descent function and get Theta.

% I have used data consisting of populations as variables and profits as target values.

% Navigate to directory where you keep these three files on octave cli and then execute.

data = load('data_.txt');

Y = data(:, 2); % Y is a column vector with all entries in second column of data matrix.

n = length(Y);

x = [ones(n, 1), data(:,1)]; % x is the matrix with first column as ones and second data.

Theta = zeros(2, 1); % Initial Theta vector with both theta values set to zero (column).

Alpha = 0.01;

noi = 1500;

Theta = gradientdescent(x, Y, Theta, Alpha, noi);

fprintf('New theta vector is below:\n');

fprintf('%f\n', Theta);

fprintf('Expected Theta vector (approx)\n');

fprintf(' -3.6303\n 1.1664\n\n');

author: TheAlgorithms

repo: MATLAB-Octave


text:
% Analysis of the Gradient of Convolution
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     01/03/2020
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0; %<! Continue from Question 1
figureCounterSpec   = '%04d';

generateFigures = ON;

CONVOLUTION_SHAPE_FULL         = 1;
CONVOLUTION_SHAPE_SAME         = 2;
CONVOLUTION_SHAPE_VALID        = 3;

DIFF_MODE_FORWARD   = 1;
DIFF_MODE_BACKWARD  = 2;
DIFF_MODE_CENTRAL   = 3;
DIFF_MODE_COMPLEX   = 4;

%% Simulation Parameters

numCoefficients = 1;
numSamples      = 14;

diffMode    = DIFF_MODE_CENTRAL;
epsVal      = 1e-6;

%% Generate Data

vH = rand(numCoefficients, 1);
vX = randn(numSamples, 1);

vY = conv2(vX, vH, 'same');

hObjFun = @(vX) 0.5 * sum((conv2(vX, vH, 'same') - vY) .^ 2);

%% Numerical Derivative
vX0 = randn(numSamples, 1);

vDRef = CalcFunGrad(vX0, hObjFun, diffMode, epsVal);
vDD = conv2(vH(end:-1:1), (conv2(vX0, vH, 'same') - vY), 'full');

firstIdx = floor((numSamples + numCoefficients - 1 - numSamples) / 2) + 1;
lastIdx = firstIdx + numSamples - 1;
vD = vDD(firstIdx:lastIdx);

max(abs(vDRef - vD))

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function p05_header ( )

%*****************************************************************************80
%
%% P05_HEADER prints some information about problem 05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    None.
%
  m = 2;

  center1 = [  0.0, 0.0 ];
  center2 = [ -0.4, 0.0 ];
  r1 = 1.0;
  r2 = 0.55;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'P05:\n' );
  fprintf ( 1, '  Strang and Persson example #5\n' );
  fprintf ( 1, '  The horn.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Circle C1 has center = (0,0)\n' );
  fprintf ( 1, '  Radius R1 = %f\n', r1 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Circle C2 has center = (-0.4,0)\n' );
  fprintf ( 1, '  Radius R2 = %f\n', r2 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Points in the region are:\n' );
  fprintf ( 1, '    in C1 and not in C2 and have 0 <= Y.\n' );
  fprintf ( 1, '  A uniform mesh density is requested.\n' );
  fprintf ( 1, '  Element sizes tried were 0.4, 0.2, 0.1.\n' );

  fprintf ( 1, '\n' );
  boundary_segment_num = p05_boundary_segment_num ( );
  fprintf ( 1, '  Number of boundary segments = %d\n', boundary_segment_num );
  fixed_num = p05_fixed_num ( );
  fprintf ( 1, '  Number of fixed points = %d\n', fixed_num );
  hole_num = p05_hole_num ( );
  fprintf ( 1, '  Number of holes = %d\n', hole_num );

  return
end


text:
function test_suite = test_sphericalAngle
%TESTSPHERICALANGLE  One-line description here, please.
%
%   output = testSphericalAngle(input)
%
%   Example
%   testSphericalAngle
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-06-21,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testP2At100(testCase) %#ok<*DEFNU>

p2 = [1 0 0];
p1 = [0 1 0];
p3 = [0 0 1];

alpha = sphericalAngle(p1, p2, p3);

testCase.assertEqual(pi/2, alpha);

% try in the other direction
alpha = sphericalAngle(p3, p2, p1);

testCase.assertEqual(3*pi/2, alpha, 'AbsTol', 1e-4);

function testP2At010(testCase)

p1 = [1 0 0];
p2 = [0 1 0];
p3 = [0 0 1];

alpha = sphericalAngle(p1, p2, p3);

testCase.assertEqual(3*pi/2, alpha);

% try in the other direction
alpha = sphericalAngle(p3, p2, p1);

testCase.assertEqual(pi/2, alpha, 'AbsTol', 1e-4);

function testSpherical(testCase)

sph1 = [.1 0];
sph2 = [0 0];
sph3 = [0 .1];

alpha = sphericalAngle(sph1, sph2, sph3);
testCase.assertEqual(pi/2, alpha, 'AbsTol', 1e-4);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ x, qraux, ipvt ] = cqrdc ( x, ldx, n, p, ipvt, job )

%*****************************************************************************80
%
%% CQRDC computes the QR factorization of an N by P complex matrix.
%
%  Discussion:
%
%    CQRDC uses Householder transformations to compute the QR factorization
%    of an N by P matrix X.  Column pivoting based on the 2-norms of the
%    reduced columns may be performed at the user's option.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex X(LDX,P), the matrix whose decomposition
%    is to be computed.
%
%    Input, integer LDX, the leading dimension of X.  N <= LDX.
%
%    Input, integer N, the number of rows of the matrix.
%
%    Input, integer P, the number of columns in the matrix X.
%
%    Input, integer IPVT(P), integers that control the
%    selection of the pivot columns.  The K-th column X(K) of X is placed
%    in one of three classes according to the value of IPVT(K):
%      IPVT(K) > 0, then X(K) is an initial column.
%      IPVT(K) == 0, then X(K) is a free column.
%      IPVT(K) < 0, then X(K) is a final column.
%    Before the decomposition is computed, initial columns are moved to the
%    beginning of the array X and final columns to the end.  Both initial
%    and final columns are frozen in place during the computation and only
%    free columns are moved.  At the K-th stage of the reduction, if X(K)
%    is occupied by a free column it is interchanged with the free column
%    of largest reduced norm.
%
%    Input, integer JOB, initiates column pivoting.
%    0, no pivoting is done.
%    nonzero, pivoting is done.
%
%    Output, complex X(LDX,P); the upper triangle contains the upper
%    triangular matrix R of the QR factorization.  Below its diagonal, X
%    contains information from which the unitary part of the decomposition
%    can be recovered.  If pivoting has been requested, the decomposition is
%    not that of the original matrix X, but that of X with its columns
%    permuted as described by IPVT.
%
%    Output, complex QRAUX(P), further information required to recover
%    the unitary part of the decomposition.
%
%    Output, integer IPVT(P); the index of the column of the
%    original matrix that has been interchanged into
%    the K-th column, if pivoting was requested.
%    IPVT is not referenced if JOB == 0.
%
  pl = 1;
  pu = 0;

  if ( job ~= 0 )
%
%  Pivoting has been requested.  Rearrange the columns according to IPVT.
%
    for j = 1 : p

      swapj = ( 0 < ipvt(j) );
      negj = ( ipvt(j) < 0 );

      if ( negj )
        ipvt(j) = -j;
      else
        ipvt(j) = j;
      end

      if ( swapj )

        if ( j ~= pl )
          temp      = x(1:n,pl);
          x(1:n,pl) = x(1:n,j);
          x(1:n,j)  = temp;
        end

        ipvt(j) = ipvt(pl);
        ipvt(pl) = j;
        pl = pl + 1;

      end

    end

    pu = p;

    for jj = 1 : p

      j = p - jj + 1;

      if ( ipvt(j) < 0 )

        ipvt(j) = -ipvt(j);

        if ( j ~= pu )

          temp      = x(1:n,pu);
          x(1:n,pu) = x(1:n,j);
          x(1:n,j)  = temp;

          i        = ipvt(pu);
          ipvt(pu) = ipvt(j);
          ipvt(j)  = i;

        end

        pu = pu - 1;

      end

    end

  end
%
%  Compute the norms of the free columns.
%
  for j = pl : pu
    qraux(j) = scnrm2 ( n, x(1:n,j), 1 );
    work(j) = qraux(j);
  end
%
%  Perform the Householder reduction of X.
%
  lup = min ( n, p );

  for l = 1 : lup
%
%  Locate the column of largest norm and bring it
%  into the pivot position.
%
    if ( pl <= l & l < pu )

      maxnrm = 0.0;
      maxj = l;

      for j = l : pu
        if ( maxnrm < real ( qraux(j) ) )
          maxnrm = real ( qraux(j) );
          maxj = j;
        end
      end

      if ( maxj ~= l )

        temp        = x(1:n,l);
        x(1:n,l)    = x(1:n,maxj);
        x(1:n,maxj) = temp;

        qraux(maxj) = qraux(l);
        work(maxj)  = work(l);

        i          = ipvt(maxj);
        ipvt(maxj) = ipvt(l);
        ipvt(l)    = i;

      end

    end

    qraux(l) = 0.0;

    if ( l ~= n )
%
%  Compute the Householder transformation for column L.
%
      nrmxl = scnrm2 ( n-l+1, x(l:n,l), 1 );

      if ( cabs1 ( nrmxl ) ~= 0.0 )

        if ( cabs1 ( x(l,l) ) ~= 0.0 )
          nrmxl = csign2 ( nrmxl, x(l,l) );
        end

        x(l:n,l) = x(l:n,l) / nrmxl;
        x(l,l) = 1.0 + x(l,l);
%
%  Apply the transformation to the remaining columns,
%  updating the norms.
%
        for j = l+1 : p

          t = - ( conj ( transpose ( x(l:n,l) ) ) * x(l:n,j) ) / x(l,l);
          x(l:n,j) = x(l:n,j) + t * x(l:n,l);

          if ( j < pl | pu < j )
            continue
          end

          if ( cabs1 ( qraux(j) ) == 0.0 )
            continue
          end

          tt = 1.0 - ( abs ( x(l,j) ) / real ( qraux(j) ) )^2;
          tt = max ( tt, 0.0 );
          t = tt;
          tt = 1.0 + 0.05 * tt * ( real ( qraux(j) ) / real ( work(j) ) )^2;

          if ( tt ~= 1.0 )
            qraux(j) = qraux(j) * sqrt ( t );
          else
            qraux(j) = scnrm2 ( n-l, x(l+1:n,j), 1 );
            work(j) = qraux(j);
          end

        end
%
%  Save the transformation.
%
        qraux(l) = x(l,l);
        x(l,l) = -nrmxl;

      end

    end

  end

  return
end


text:
function fg = dtiFiberCombineFiberGroups(fgList,outName,type)
%
% fg = dtiFiberCombineGroups([fgList=uiselectfile],...
%                            [outName='combinedFiberGroup'],[type='pdb'])
%
% Combine N fiber groups (pdb or mat) into a single fiber group and save.
% This routine is a wrapper for dtiMergeFiberGroups, making it simple to
% comine an arbitrary number of groups.
%
% INPUTS:
%   fgList  - CELL ARRAY of fiber names.
%   outName - Name of fiber group to save.
%   type    - 'mat' or 'pdb'
%
%
% (C) Stanford University, VISTA Lab
%

%% Check INPUTS

if notDefined('fgList')
    [names p] = uigetfile({'*.pdb';'*.mat'}, 'Select fiber groups','MultiSelect','on');
    for ff=1:numel(names)
        fgList{ff} = fullfile(p, names{ff});
    end
end

if ~iscell(fgList)
    error('fgList must be a cell array of fiber names');
end

if notDefined('type')
    type = 'pdb';
end

if notDefined('outName')
        outName = fullfile(pwd,['combinedFiberGroup.' type]);
        name = 'combinedFiberGroup';
else
    [~, name ~] = fileparts(outName);
end

%% Combine the fiber groups

disp('Combining fiber groups...');

for ii=1:numel(fgList)
    if ii==1
        fg = fgRead(fgList{ii});
    else
        fgM = fgRead(fgList{ii});
        fg = dtiMergeFiberGroups(fg,fgM);
    end

end

fg.name = name;

%% Save

fprintf('\n Saved: %s\n', outName);
fgWrite(fg,outName,type);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [stc] = mne_read_stc_file(filename)
%
% [stc] = mne_read_stc_file(filename)
%
% Reads an stc file. The returned structure has the following fields
%
%     tmin           The first time point of the data in seconds
%     tstep          Time between frames in seconds
%     vertices       vertex indices (0 based)
%     data           The data matrix (nvert * ntime)
%
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.7  2006/04/23 15:29:41  msh
%   Added MGH to the copyright
%
%   Revision 1.6  2006/04/10 23:26:54  msh
%   Added fiff reading routines
%
%   Revision 1.5  2005/12/05 20:23:21  msh
%   Added fiff_save_evoked. Improved error handling.
%
%   Revision 1.4  2005/11/21 05:40:24  msh
%   Required number of arguments was wrong.
%
%   Revision 1.3  2005/11/21 03:19:12  msh
%   Improved error handling
%
%   Revision 1.2  2005/11/21 02:15:51  msh
%   Added more routines
%
%   Revision 1.1  2005/11/21 01:41:57  msh
%   Introduced structures and start all function names with mne_
%
%
me='MNE:mne_read_stc_file';
if(nargin ~= 1)
   error(me,'usage: [stc] = mne_read_stc_file(filename)');
end

[fid,message] = fopen(filename,'r','ieee-be');
if fid == -1
    error(me,message)
end

STATUS = fseek(fid, 0, 'eof');
filelen= ftell(fid);
STATUS = fseek(fid, 0, 'bof');

% read tmin in ms
[stc.tmin count]=fread(fid,1,'float32');
stc.tmin=stc.tmin/1000.0;
% read sampling rate in ms
[stc.tstep count]=fread(fid,1,'float32');
stc.tstep=stc.tstep/1000.0;
% read number of vertices/sources
[vertices_n count]=fread(fid,1,'uint32');
% read the source vector
[stc.vertices count]=fread(fid,vertices_n,'uint32');
stc.vertices=uint32(stc.vertices);
% read the number of timepts
[data_n count]=fread(fid,1,'uint32');

if (mod((filelen/4-4-vertices_n),data_n*vertices_n) ~= 0)
   error(me,'incorrect stc file size')
end

% read the data matrix
[stc.data count]=fread(fid,vertices_n*data_n,'float32');
stc.data=squeeze(reshape(stc.data,vertices_n,data_n));
% close the file
fclose(fid);
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [] = test()

x0 = [3; 2];    % Inital value
bl = [-2; 0];   % lower bounds
bu = [inf; 3];  % upper bounds

param = [];
param.maxIter = 10;     % max number of iterations
param.maxFnCall = 100;  % max number of calling the function
param.relCha = 1e-5;      % tolerance of constraint satisfaction
param.tolPG = 1e-5;   % final objective function accuracy parameter
param.m = 10;

function [f, g] = toy_func(x)
    cen = [-1; 1];
    f = 0.5 * norm(x - cen)^2;
    g = x - cen;
end

% At present, only box constraints are implemented in the mex.
[x, f, iter, numCall, flag] = lbfgsb(x0, bl, bu, @toy_func, [], @genericcallback, param) %

end

author: HKUST-KnowComp

repo: FMG

save_path: github-repos/MATLAB/HKUST-KnowComp-FMG


text:
% 此函数非本人编写
% 参考自：https://www.cnblogs.com/GarfieldEr007/p/5598861.html
function [ dest ] = drawRect( src, pt, wSize,  lineSize, color )
%简介：
% %将图像画上有颜色的框图，如果输入是灰度图，先转换为彩色图像，再画框图
% 图像矩阵
% 行向量方向  是  y
% 列向量方向  是  x
%----------------------------------------------------------------------
%输入：
% src：        原始图像，可以为灰度图，可为彩色图
% pt：         左上角坐标   [x1, y1]
% wSize：   框的大小      [wx, wy]
% lineSize： 线的宽度
% color：     线的颜色      [r,  g,  b]
%----------------------------------------------------------------------
%输出：
% dest：           画好了的图像
%----------------------------------------------------------------------

%flag=1: 有缺口的框
%flag=2: 无缺口的框
flag = 1;

%判断输入参数个数
if nargin < 5
    color = [255 255 0];
end

if nargin < 4
    lineSize = 1;
end

if nargin < 3
    disp('输入参数不够 !!!');
    return;
end

%判断框的边界问题
[yA, xA, z] = size(src);
x1 = pt(1);
y1 = pt(2);
wx = wSize(1);
wy = wSize(2);
if  x1>xA || ...
        y1>yA||...
        (x1+wx)>xA||...
        (y1+wy)>yA

    disp('画的框将超过图像 !!!');
    return;
end

%如果是单通道的灰度图，转成3通道的图像
if 1==z
    dest(:, : ,1) = src;
    dest(:, : ,2) = src;
    dest(:, : ,3) = src;
else
    dest = src;
end

%开始画框图
for c = 1 : 3                 %3个通道，r，g，b分别画
    for dl = 1 : lineSize   %线的宽度，线条是向外面扩展的
        d = dl - 1;
        if  1==flag %有缺口的框
            try
                dest(  y1-d ,            x1:(x1+wx) ,  c  ) =  color(c); %上方线条
                dest(  y1+wy+d ,     x1:(x1+wx) , c  ) =  color(c); %下方线条
                dest(  y1:(y1+wy) ,   x1-d ,           c  ) =  color(c); %左方线条
                dest(  y1:(y1+wy) ,   x1+wx+d ,    c  ) =  color(c); %左方线条
            catch
                continue
            end
        elseif 2==flag %无缺口的框
            dest(  y1-d ,            (x1-d):(x1+wx+d) ,  c  ) =  color(c); %上方线条
            dest(  y1+wy+d ,    (x1-d):(x1+wx+d) ,  c  ) =  color(c); %下方线条
            dest(  (y1-d):(y1+wy+d) ,   x1-d ,           c  ) =  color(c); %左方线条
            dest(  (y1-d):(y1+wy+d) ,   x1+wx+d ,    c  ) =  color(c); %左方线条
        end
    end
end %主循环尾

end %函数尾

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [foundFiles, iDepth] = file_find( baseDir, filePattern, maxDepth, isSingle, iDepth )
% FILE_FIND: Find a file recursively.
%
% USAGE:  [foundFiles, iFoundDepth] = file_find( baseDir, filePattern, maxDepth=Inf, isSingle=1 )
%
% INPUT:
%    - baseDir     : Full path to the directory to search
%    - filePattern : Name of the target file (wild chars allowed)
%    - maxDepth    : Maximum folder depth from the baseDir
%    - isSingle    : If 1, only search for one file, exit when it's found
%    - iDepth      : Current recursion level
%
% OUTPUT:
%    - foundFiles : Full path to files, or [] if no file were found

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

% Parse inputs
if (nargin < 5) || isempty(iDepth)
    iDepth = 1;
end
if (nargin < 4) || isempty(isSingle)
    isSingle = 1;
end
if (nargin < 3) || isempty(maxDepth)
    maxDepth = Inf;
end
if (nargin < 2)
    error('foundFiles = file_find( baseDir, filePattern, maxDepth, isSingle );');
end
% Default return value
if isSingle
    foundFiles = [];
else
    foundFiles = {};
end

% Base dir not valid
if isempty(baseDir) || ~isdir(baseDir) || (baseDir(1) == '.')
    return;
% Pattern is a subfolder of the base dir
elseif isSingle && ~any(filePattern == '*') && isdir(fullfile(baseDir, filePattern))
    foundFiles = fullfile(baseDir, filePattern);
    return;
else
    % Try to find required file
    listDir = dir(fullfile(baseDir, filePattern));
    if ~isempty(listDir)
        if isSingle
            foundFiles = fullfile(baseDir, listDir(1).name);
        else
            for i = 1:length(listDir)
                foundFiles = cat(2, foundFiles, {fullfile(baseDir, listDir(i).name)});
            end
        end
        return
    end
    % If reached the recursion limit
    if (maxDepth <= 1)
        return
    end
    % Get subdirectories
    listDir = dir(fullfile(baseDir, '*'));
    listDir([listDir.isdir] == 0) = [];
    % Process each subdirectory
    for i = 1:length(listDir)
        if (listDir(i).name(1) ~= '.')
            newDir = fullfile(baseDir, listDir(i).name);
            [foundRec, iFoundDepth] = file_find(newDir, filePattern, maxDepth - 1, isSingle, iDepth + 1);
            if ~isempty(foundRec)
                iDepth = iFoundDepth;
                % Looking for one file only: return
                if isSingle
                    foundFiles = foundRec;
                    return;
                % Otherwise: Keep searching
                else
                    foundFiles = cat(2, foundFiles, foundRec);
                end
            end
        end
    end
end
iDepth = [];

author: brainstorm-tools


text:
function [names, ids, dat, descrips] = concatenate(D, varargin)
% Concatenates Subject-level and Event-level data across all subjects
%
% :Usage:
% ::
%
%    [names ids dat] = concatenate(D, [optional inputs])
%
% :Inputs:
%
%   **D:**
%        canlab_dataset object
%
% :Optional Inputs:
%   **a logical array**
%       a vector of 1/0 values to use as wh_keep
%
%
% :Outputs:
%
%   **names:**
%        cell array of variable names
%
%   **ids:**
%        subject IDs matching data rows in dat
%
%   **dat:**
%        subjects*events x variables matrix
%           - subject number, event number are included
%           - all subject-level and event-level data are included
%           - this format appropriate for, e.g., SAS/HLM
%
%   **descrips:**
%        cell array of variable descriptions
%
% :Examples:
% ::
%
%    [names, ids, flatdat] = concatenate(D);
%    id_numbers = flatdat(:, 1);
%
%    wh_subjs = true(size(D.Subj_Level.id));
%    wh_subjs([13 18 19]) = false;
%    [names, ids, dat] = concatenate(D, wh_subjs);
%

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

wh_ids = true(length(D.Subj_Level.id), 1); % select subjects
if length(varargin) > 0, wh_ids = logical(varargin{1}); end
whfind = find(wh_ids);

names = {'subj_num' D.Subj_Level.names{:} 'Event_number' D.Event_Level.names{:}};
descrips = {'subj_num' D.Subj_Level.descrip{:} 'Event_number' D.Event_Level.descrip{:}};

n = length(D.Subj_Level.id(wh_ids));

x = D.Event_Level.data(wh_ids);

for i = 1:n

    e = size(D.Event_Level.data{whfind(i)}, 1);

    ids{i} = repmat(D.Subj_Level.id{whfind(i)}, e, 1);

    %           subj ID                     subj level data              trialID  event data
    x{i} = [repmat(i, e, 1) repmat(D.Subj_Level.data(whfind(i), :), e, 1) (1:e)' x{i}];

end

ids = strvcat(ids{:});

dat = cat(1, x{:});

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [train_eval,test_eval]  = ml_kcv_regression_eval( X,y,g,train,test,k,train_eval,test_eval)
%ML_KCV_REGRESSION_EVAL
%
%   Evaluation of the Goodness of Fit of your regression model. Have look
%   at http://uk.mathworks.com/help/curvefit/evaluating-goodness-of-fit.html
%   for a reference.
%
%  input ------------------------------------------------------------------
%
%       o   X        : (N x D),  input variables
%
%       o   y        : (N x 1),  output variables
%
%       o   g        : function handle, y = g(x)
%
%       o   train    : (M x 1),  set of indicies of X to use as train data.
%
%       o   test     : (P x 1),  set of indicies of X to use for testing.
%
%       o   k        : (1 x 1),  current iteration of k-fold
%

%   output ----------------------------------------------------------------
%
%       o   train_eval : struct
%
%       o   test_eval  : struct

%% Evaluate regression function on training data

hy = g(X(train(:),:));

gf = gfit2(y(train(:)),hy,'all');

train_eval.mse(k)   = gf(1);    %  '1'  - mean squared error (mse)
train_eval.nmse(k)  = gf(2);    %  '2'  - normalised mean squared error (nmse)
train_eval.rmse(k)  = gf(3);    %  '3'  - root mean squared error (rmse)
train_eval.nrmse(k) = gf(4);    %  '4'  - normalised root mean squared error (nrmse)
train_eval.mae(k)   = gf(5);    %  '5'  - mean absolute error (mae)
train_eval.mare(k)  = gf(6);    %  '6'  - mean  absolute relative error  (mare)
train_eval.r(k)     = gf(7);    %  '7'  - coefficient of correlation (r)
train_eval.d(k)     = gf(8);    %  '8'  - coefficient of determination (d)
train_eval.e(k)     = gf(9);    %  '9'  - coefficient of efficiency (e)
train_eval.me(k)    = gf(10);   %  '10' - maximum absolute error
train_eval.mre(k)   = gf(11);   %  '11' - maximum absolute relative error

%% Evaluate regression function on test data

hy = g(X(test(:),:));

gf = gfit2(y(test(:)),hy,'all');

test_eval.mse(k)   = gf(1);    %  '1'  - mean squared error (mse)
test_eval.nmse(k)  = gf(2);    %  '2'  - normalised mean squared error (nmse)
test_eval.rmse(k)  = gf(3);    %  '3'  - root mean squared error (rmse)
test_eval.nrmse(k) = gf(4);    %  '4'  - normalised root mean squared error (nrmse)
test_eval.mae(k)   = gf(5);    %  '5'  - mean absolute error (mae)
test_eval.mare(k)  = gf(6);    %  '6'  - mean  absolute relative error  (mare)
test_eval.r(k)     = gf(7);    %  '7'  - coefficient of correlation (r)
test_eval.d(k)     = gf(8);    %  '8'  - coefficient of determination (d)
test_eval.e(k)     = gf(9);    %  '9'  - coefficient of efficiency (e)
test_eval.me(k)    = gf(10);   %  '10' - maximum absolute error
test_eval.mre(k)   = gf(11);   %  '11' - maximum absolute relative error

end

author: epfl-lasa


text:
function [newTree, derTree] = splitTree(tree, maxOrder)
%SPLITTREE   Split syntax trees into derivative part and non-derivative part
%   Calling sequence:
%      [NEWTREE, DERTREE] = SPLITTREE(TREE, MAXORDER)
%   where the inputs are:
%      TREE:       The syntax tree to be split.
%      MAXORDER:   A vector that contains the the maximum differential order of
%                  each variable that appears in the problem.
%   and the outputs are
%      NEWTREE:    A syntax tree which describes the factor in which the highest
%                  order derivative appears. E.g. if we split the expression
%                  5*diff(u) + sin(u), NEWTREE is the syntax tree corresponding
%                  to 5*diff(u).
%      DERTREE:    A syntax tree which describes the remaining factors not
%                  included in NEWTREE. In the example above, DERTREE is the
%                  syntax tree corresponding to sin(u).

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Find what variables we have actually computed derivatives of.
diffVar = maxOrder > 0;

% If our input is not a syntax tree, we can be certain we don't need to split it
% (this can e.g. happen if TREE is a CHEBFUN or a scalar).
if ( ~isstruct(tree) )
    newTree = tree;
    derTree = [];
    return
end

% Find the diffOrders in the input tree.
treeDiffOrder = tree.diffOrder;

% If our tree doesn't have any derivatives of highest order that appear in the
% problem, we can return the subtree. But only need to check for those variables
% that we actually differentiate with respect to, e.g. if maxOrder = 0, we can
% return the subtree as well.
if ( all(treeDiffOrder(diffVar) < maxOrder(diffVar)) )
    newTree = tree;
    derTree = [];
    return
end

% Split the tree recursively, depending on how many inputs the operator takes.
switch tree.numArgs
    case 1
        % We're dealing with a unary operator.
        if ( strcmp(tree.method, 'uminus') )
            % For a unary minus, convert it to a "-1*" tree and split, which
            % will take care of negating the resulting trees:
            tempTree = struct('method', 'times', 'numArgs', 2, ...
                'left', -1, 'right', tree.center, ...
                'diffOrder', tree.diffOrder, 'height', tree.height, ...
                'ID', tree.ID, 'hasTerms', tree.hasTerms);
            [newTree, derTree] = treeVar.splitTree(tempTree, maxOrder);
        else
            % Otherwise, we can simply split the center node:
            [newTree, derTree] = treeVar.splitTree(tree.center, maxOrder);
        end
    case 2
        if ( any(strcmp(tree.method, {'diff', 'times', 'rdivide'})) )
            % We have already been through the expandTree() method, which
            % guarantees that there are no terms left in the tree which include
            % highest order derivatives along with other expressions in which
            % the unknown variable appears, e.g. 5*(diff(u) + u). Since we don't
            % reach this point in the code unless the input TREE contains the
            % maximum diffOrder of the problem, we can safely assume that the
            % input TREE is the derivative tree, and return it.
            newTree = [];
            derTree = tree;
        else
            % We're at + or -, split the left and right children trees
            % recursively.
            [newTreeLeft, derTreeLeft] = ...
                treeVar.splitTree(tree.left, maxOrder);
            [newTreeRight, derTreeRight] = ...
                treeVar.splitTree(tree.right, maxOrder);

            if ( isempty(newTreeLeft) )
                % We only had a derivative part on the left. Thus, the
                % non-derivative parts will only consist of the non-derivative
                % part of the right tree. However, if the method of our current
                % tree is MINUS(), we must add a UMINUS() in front of
                % NEWTREERIGHT before we can return it. That will be done in the
                % ONETREEFROMRIGHT() method. Simply put, we're converting a
                % binary minus between an empty tree on left and a non-empty
                % tree on the right to a UMINUS on the non-empty tree.
                newTree = oneTreeFromRight(newTreeRight, tree.method);
            elseif (isempty(newTreeRight) )
                % We only had a derivative part on the right, so the
                % non-derivative part only consists of the left part.
                newTree = newTreeLeft;
            elseif ( ~isstruct(newTreeLeft) && ~isstruct(newTreeRight) )
                % Both left and right trees were actually a CHEBFUN or scalar,
                % combine them in a single CHEBFUN/scalar and return as the
                % NEWTREE.
                newTree = eval([tree.method, '(newTreeLeft, newTreeRight)']);
            else
                % Had non derivative parts on both left and right, potentially a
                % combination of CHEBFUN/scalars and syntax trees. Need to
                % combine them.
                if ( ~isstruct(newTreeLeft) )
                    % Left tree was actually a CHEBFUN or scalar, so the new
                    % diffOrders and height only depend on the right tree.
                    newDiffOrder = newTreeRight.diffOrder;
                    newHeight = newTreeRight.height;
                elseif ( ~isstruct(newTreeRight) )
                    % Right tree was actually a CHEBFUN or scalar, so the new
                    % diffOrders and height only depend on the right tree.
                    newDiffOrder = newTreeLeft.diffOrder;
                    newHeight = newTreeLeft.height;
                else
                    % Both left and right tree were actually syntax trees, so
                    % find the maximum heights and diffOrders.
                    newDiffOrder = max(newTreeLeft.diffOrder, ...
                        newTreeRight.diffOrder);
                    newHeight = max(newTreeLeft.height, newTreeRight.height);
                end

                % Construct a new syntax tree from the new left and right child
                % trees (the non-derivative parts).
                newTree = struct('method', tree.method, ...
                    'numArgs', tree.numArgs, ...
                    'left', newTreeLeft, 'right', newTreeRight, ...
                    'diffOrder', newDiffOrder, ...
                    'height', newHeight);
            end

            if ( isempty(derTreeLeft) )
                % Left child tree only consisted of non-derivative part. Like
                % above, return the right tree, but add a UMINUS() on top of it
                % if needed.
                derTree = oneTreeFromRight(derTreeRight, tree.method);
            elseif (isempty(derTreeRight) )
                % Right child tree only consisted of non-derivative part.
                derTree = derTreeLeft;
            else
                % Combine the left and right derivative subtrees.
                derTree = struct('method', tree.method, ...
                    'numArgs', tree.numArgs, ...
                    'left', derTreeLeft, 'right', derTreeRight, ...
                    'diffOrder', max(derTreeLeft.diffOrder, ...
                    derTreeRight.diffOrder), ...
                    'height', max(derTreeLeft.height, derTreeRight.height));
            end

        end
end
end

function ot = oneTreeFromRight(tree, operator)
%ONETREEFROMRIGHT   Convert binary minus to uminus when left tree is empty.
%   Add a UMINUS in front of a syntax tree in the case where the left child tree
%   of MINUS() only consisted of the derivative part. Simply put, we're
%   converting a binary minus between an empty tree on left and a non-empty tree
%   on the right to a UMINUS on the non-empty tree.
if ( strcmp(operator, 'minus') )
    % Only need to worry if the original operator was a -.
    if ( isstruct(tree) )
        % If the input is a syntax tree, add a UMINUS on the top of it.
        ot = struct('method', 'uminus', 'numArgs', 1, 'center', tree, ...
            'diffOrder', tree.diffOrder, 'ID', tree.ID, ...
            'height', tree.height + 1);
    else
        % Input was a CHEBFUN or a scalar, can simply negate it.
        ot = -tree;
    end
else
    % Didn't have a MINUS(), so can simply return the input.
    ot = tree;
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function this = hs_optical_flow(varargin)
%
%HS_OPTICAL_FLOW   Optical flow computation with Horn & Schunck method
%       B.K.P. Horn and B.G. Schunck. Determining optical flow.
%       Artificial Intelligence, 16:185-203, Aug. 1981.
%       http://people.csail.mit.edu/bkph/papers/Optical_Flow_OPT.pdf
%
%   HS_OPTICAL_FLOW([IMGS]) constructs a HS optical flow object
%   with the optional image sequence IMGS ([n x m x 2] array).
%   HS_OPTICAL_FLOW(O) constructs HS optical flow object by copying O.
%
%   This is a member function of the class 'hs_optical_flow'.

%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-10-30 $
%   $Revision: $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

error(nargchk(0, 1, length(varargin)));

  switch (length(varargin))
    case 0

      this.images               = [];
      this.lambda               = 80;   %5e-2;
      this.lambda_q             = 80;   %not used, for consistency with other program
      this.gnc_iters            = 1;    %not used, for consistency with other program
      this.pyramid_levels       = 4; %

      this.pyramid_spacing      = 2; %         % default to the prolongation/restriction filter

      this.max_warping_iters    = 10;           % # of warping/linearization per pyramid level
      this.median_filter_size   = [];
      this.texture              = false;        % apply to the image pyramid;

      this.solver               = 'backslash';  % pcg, sor
%       this.solver               =  'sor';      % MAYANK's.
      this.sor_max_iters        = 1e1;          % 100 seems sufficient
      this.interpolation_method = 'cubic';      % 'bi-cubic', 'cubic', 'bi-linear'
      this.deriv_filter         = [1 -8 0 8 -1]/12;
      this.display              = false;
      this.limit_update         = true;
      this.sigmaD2              = 1;                 % data term
      this.sigmaS2              = 1;                 % spatial term

      this.sigmaP               = 1;                 % sigma of Gaussian for performing "texture" decomposition
      this.weightP              = 1;                 % images -weightPxGaussian*images

      this.mf_iter              = 1;

      method = 'quadratic';
      this.spatial_filters = {[1 -1], [1; -1]};
      for i = 1:length(this.spatial_filters);
          this.rho_spatial_u{i}   = robust_function(method, 1); % 0.1
          this.rho_spatial_v{i}   = robust_function(method, 1);
      end;
      this.rho_data        = robust_function(method, 1); % 6.3

      this.color_images     = [];

      this = class(this, 'hs_optical_flow');

    case 1
      if isa(varargin{1}, 'hs_optical_flow')
        this = other;

      else
          this = hs_optical_flow;
      end

    otherwise

      error('Incompatible arguments!');

  end

author: kristinbranson

repo: JAABA


text:
function a = hankel ( n, x )

%*****************************************************************************80
%
%% HANKEL returns a HANKEL matrix.
%
%  Formula:
%
%    A(I,J) = X(I+J-1)
%
%  Example:
%
%    N = 5, X = ( 1, 2, 3, 4, 5, 6, 7, 8, 9 )
%
%    1  2  3  4  5
%    2  3  4  5  6
%    3  4  5  6  7
%    4  5  6  7  8
%    5  6  7  8  9
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is a Hankel matrix: constant along anti-diagonals.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X(2*N-1), the vector that defines A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for j = 1 : n
    a(1:n,j) = x(j:j+n-1);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ResultsAllCellLines,OverViewResults] = setQuantConstraints(model, samples, tol, minGrowth, obj, no_secretion, no_uptake, medium, addExtraExch, addExtraExch_value, path)
% This function takes a model and quantitative extracellular metabolomic data and returns a model in which the data is integrated as constraints described
% in "MetaboTools: A Comprehensive Toolbox for Analysis of Genome-Scale Metabolic Models", 2016, Front. Physiol., and the supplemental data "Tutorial II: Workflow for the integration
% of quantitative extracellular metabolomic data into the network context." (Data Sheet 3.pdf)
%
% USAGE:
%
%    [ResultsAllCellLines, OverViewResults] = setQuantConstraints(model, samples, tol, minGrowth, obj, no_secretion, no_uptake, medium, addExtraExch, addExtraExch_value, path, epsilon)
%
% INPUTS:
%       model:                Global metabolic model (Recon)
%       samples:              Vector specifying the samples used (there must be an output file of function .... for each sample)
%       tol:                  Cutoff value for small numbers (e.g., -1e-8). All number smaller than tol will be treated as zero
%       minGrowth:            Will be the lower bound of the objective function (e.g., 0.008). Forces the output model(s) to be able to produce a minimal objective value
%       obj:                  Objective function, e.g. `biomass_reaction2`
%       no_secretion:         Define metabolites that should not be secreted (e.g., {`EX_o2(e)`})
%       no_uptake:            Define metabolites that should not be consumed (e.g., {`EX_o2s(e)`, `EX_h2o2(e)`})
%       medium:               Define if certain exchanges should be excluded from minimization of exchanges (e.g., {}, if no medium except the exometabolomic data has been defined)
%       addExtraExch:         After adding secretions, models are still not growing, this variable allows one to recover exchanges with a defined small value
%       addExtraExch_value:   e.g. 1 as arbitrary small flux value / the resulting ub = 1, lb = -1.
%       path:                 Location of the .mat files for samples.
%
%
% OUTPUTS:
%       ResultsAllCellLines:  Structure that contains pruned and unpruned model, Vector of the Exchange_reactions, the exchange reactions added by `minExCard`, `minFLux` and `maxFlux` of the added reactions, maximal objective value, and the results of the gene deletion
%       OverViewResults:      Overview of model statistics, e.g., number of reactions, metabolites, genes, number of essential genes, min and max objective values for easy comparison between sets of models
%
%
% Depends on `changeRxnBounds`, `fluxVariability`, `optimizeCbModel`, `generateCompactExchModel` as well its dependent functions `pruneModel`, `findMinCardModel`, `findOptExchRxns`
%
% .. Author: - Maike K. Aurich 18/02/15

cntO=1;
% Set overview variable
OverViewResults{1,cntO} ='cell line';cntO = cntO+1;
OverViewResults{1,cntO} = 'num Added Rxns';cntO = cntO+1;
OverViewResults{1,cntO} = 'num rxns pruned model';cntO = cntO+1;
OverViewResults{1,cntO} = 'num mets pruned model';cntO = cntO+1;
OverViewResults{1,cntO} = 'num genes pruned model';cntO = cntO+1;
OverViewResults{1,cntO} = 'max growth rate';cntO = cntO+1;
OverViewResults{1,cntO} = 'num exch rxns';cntO = cntO+1;
OverViewResults{1,cntO} = 'O2 requirement';cntO = cntO+1;
OverViewResults{1,cntO} = 'extra Exchanges Added';cntO = cntO+1;
OverViewResults{1,cntO} = 'num recovered exchanges';cntO = cntO+1;

for j = 1:length(samples)

    j
    %   load([path filesep 'model.mat']); %% uncomment for Recon1
    ExtraExchAdded = 0;

    FILENAME = char(samples(j,1));
    load([path filesep FILENAME '.mat'])%, 'uptake_value', 'secr_value' , 'uptake' ,'secretion', 'cell_line' );

    % Map exo-metabolomic data
    model2=model;
    model2.lb(find(ismember(model2.rxns,obj)))=minGrowth;% based on slowlest cell line in data

    % Map uptake
    for k=1:length(uptake)
        Uptake_rxns=uptake(k);
        lb = uptake_value(k,2);
        ub = uptake_value(k,3);
        model2 = changeRxnBounds(model2,Uptake_rxns,ub,'u'); %enforce uptake of metabolites taken up by cells in the experiment
        model2 = changeRxnBounds(model2,Uptake_rxns,lb,'l'); %enforce uptake of metabolites taken up by cells in the experiment

    end
    clear A sol
    model2ori=model2;
    SecretionRxnsRecovered =[];
    cnt2=1;

    % Map secretion
    for k=1:length(secretion)
        secretion_rxns= secretion(k);
        lb = secr_value(k,3);
        ub = secr_value(k,2);
        model2 = changeRxnBounds(model2, secretion_rxns,lb,'l');% enforce secretion of metabolites that are secreted by the cells in the experiment
        model2 = changeRxnBounds(model2, secretion_rxns,ub,'u');% enforce secretion of metabolites that are secreted by the cells in the experiment

        % check if model works after adding constraints
        SOL = optimizeCbModel(model2);
        if SOL.f<=abs(tol)
            sol(k,1)=SOL.f;
            secretion_rxns
            model2=model2ori; % model2ori is overwritten after each successful iteration so its not the original model (without any secretions)
            model2 = changeRxnBounds(model2, secretion_rxns,ub,'u');% set secretion of metabolite to measured upper bound nevertheless
            SecretionRxnsRecovered{cnt2,1}=secretion_rxns{1,:}; cnt2 = cnt2 +1;
        else
            sol(k,1)=SOL.f;
            model2ori=model2;
        end
    end

    % sets bounds to the lowest measured uptake/secretion rate
    for k=1:length(No_upt_secr)
        No_upt_secr_rxns= No_upt_secr(k);
        bound = min(abs([max(uptake_value(:,3));min(secr_value(:,3))]));
        model2 = changeRxnBounds(model2, No_upt_secr_rxns,-bound,'l');
        model2 = changeRxnBounds(model2, No_upt_secr_rxns,bound,'u');% sets bounds to the lowest measured uptake/secretion rate
    end

    % do not allow secretion or uptake of certain metabolites
    model2 = changeRxnBounds(model2, no_secretion,0,'u');% enforce secretion of metabolites that are secreted by the cells in the experiment
    model2 = changeRxnBounds(model2, no_uptake,0,'l');% enforce secretion of metabolites that are secreted by the cells in the experiment

    %% test if model can grow AT ALL
    SOL = optimizeCbModel(model2);
    if abs(SOL.f) < abs(tol)

        model2 = changeRxnBounds(model2, addExtraExch,-addExtraExch_value,'l');% sets bounds to the lowest measured uptake/secretion rate
        model2 = changeRxnBounds(model2, addExtraExch,addExtraExch_value,'u');% sets bounds to the lowest measured uptake/secretion rate
        ExtraExchAdded = 1;
    end

    SOL = optimizeCbModel(model2);
    if abs(SOL.f) > abs(tol)
        %Generate submodel
        [modelMin, modelPruned, Ex_Rxns] = generateCompactExchModel(model2,minGrowth);

        %added 23/07/2015
        modelPruned.c = zeros(length(modelPruned.rxns),1);
        modelPruned.c(find(ismember(modelPruned.rxns,obj)),1)=1;
        %added 23/07/2015

        sol= optimizeCbModel(modelPruned);
        Ex_RxnsAdded = Ex_Rxns;

        % remove reactions that are in uptake and secretion
        US = unique([secretion;uptake]);
        Ex_RxnsAdded(ismember(Ex_RxnsAdded,US))=[];
        [a(:,1),a(:,2)]=fluxVariability(modelPruned,1,[],Ex_RxnsAdded);

        % print Results

        ResultsAllCellLines.(samples{j}).modelPruned = modelPruned;
        ResultsAllCellLines.(samples{j}).modelMin = modelMin;
        ResultsAllCellLines.(samples{j}).Ex_Rxns = Ex_Rxns;
        ResultsAllCellLines.(samples{j}).Ex_RxnsAdded = Ex_RxnsAdded;
        ResultsAllCellLines.(samples{j}).MinMaxAddedRxns = a;
        ResultsAllCellLines.(samples{j}).maxBiomass = sol;

        ResultsAllCellLines.(samples{j}).SecretionRxnsRecovered = SecretionRxnsRecovered;

        cntO=1;
        OverViewResults{j+1,cntO} = samples{j};cntO = cntO+1;% cell line
        OverViewResults{j+1,cntO} = num2str(length(ResultsAllCellLines.(samples{j}).Ex_RxnsAdded));cntO = cntO+1;% num Added Rxns
        OverViewResults{j+1,cntO} = num2str(length(ResultsAllCellLines.(samples{j}).modelPruned.rxns));cntO = cntO+1;% num rxns
        OverViewResults{j+1,cntO} = num2str(length(ResultsAllCellLines.(samples{j}).modelPruned.mets));cntO = cntO+1;% num mets
        OverViewResults{j+1,cntO} = num2str(length(ResultsAllCellLines.(samples{j}).modelPruned.genes));cntO = cntO+1;% num genes
        OverViewResults{j+1,cntO} = num2str(ResultsAllCellLines.(samples{j}).maxBiomass.f);cntO = cntO+1;% max growth rate
        OverViewResults{j+1,cntO} = num2str(length(ResultsAllCellLines.(samples{j}).Ex_Rxns));cntO = cntO+1;% num Exchange rxns
        if ~isempty(strmatch('EX_o2(e)',ResultsAllCellLines.(samples{j}).Ex_RxnsAdded))
            OverViewResults{j+1,cntO} = num2str(1);cntO = cntO+1;% O2 requirement
        else
            OverViewResults{j+1,cntO} = num2str(0);cntO = cntO+1;% O2 requirement
        end

        OverViewResults{j+1,cntO} = num2str(ExtraExchAdded);cntO = cntO+1;% num Exchange rxns
        OverViewResults{j+1,cntO} = num2str(length(SecretionRxnsRecovered));cntO = cntO+1;% num Exchange rxns

    else
        ResultsAllCellLines.(samples{j}).model = [];
        ResultsAllCellLines.(samples{j}).AddedExchange = [];
        ResultsAllCellLines.(samples{j}).maxBiomass = [];
        cntO=1;
        OverViewResults{j+1,cntO} = samples{j};cntO = cntO+1;% cell line
    end
    clear sol Extra* addExtraExch* model2 modelMin* modelU* modelP* model2ori No* upt* secr* gr* Rxn* cnt* Blocked* h* FBA* j k m bound ans lb Close* Ex_Rxns1* Ex_Rxns2*  Sol U* Ori* i* ma* rev* t un* us* w a SOL a1 ub ToDe* ReP* Ex_* Secr* AddedExchang*
    save([path filesep 'setQuantConstraints.mat'], '-v7.3');
end

end

author: opencobra

repo: cobratoolbox


text:
function y = qinvjmul(labx,frmx,b,K)
% y = qinvjmul(labx,frmx,b,K)
%
% QINVJMUL  Inverse of Jordan multiply for Lorentz blocks
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

lorN = length(K.q);
if lorN == 0
    y = zeros(0,1);
    return
end
if length(labx) > 2*lorN
    labx = labx(K.l+1:K.l+2*lorN);
end
detx = labx(1:lorN) .* labx(lorN+1:end);
x = qframeit(labx,frmx,K);
ix = K.mainblks;
if length(b) == ix(3)-ix(1);   % lorentz only ?
    ix = (1-ix(1)) + ix;
end
% ------------------------------------------------------------
% Let y1(k) = xk'Jbk/(sqrt2*detxk)
% ------------------------------------------------------------
y1 = x(1:lorN).*b(ix(1):ix(2)-1) - ddot(x(lorN+1:end),b,K.qblkstart);
y1 = y1./(sqrt(2)*detx);
% ------------------------------------------------------------
% Let y2[k] = (sqrt2/x1)*b2[k] - (y1/x1) * x2[k]
% ------------------------------------------------------------
y = [y1; qblkmul(sqrt(2)./x(1:lorN),b,K.qblkstart)...
    - qblkmul(y1./x(1:lorN),x(lorN+1:end),K.qblkstart)];

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   R = euler2rot(abc, convention)
%   Returns the Rotation angle for three euler angles abc=[alpha beta gamma].
%
%   The convention specifies the order and axes of the rotations. For
%   example, for the 'XYZ' convention, we have
%   R = Rot(alpha,'x')*Rot(beta,'y')*Rot(gamma,'z')
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche. email:
%   arturo.gil@umh.es date:   11/11/2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function [R] = euler2rot(abc, convention)

if convention =='XYZ'
    ax = abc(1);
    ay = abc(2);
    az = abc(3);

    Rx = Rot(ax, 'x');
    Ry = Rot(ay, 'y');
    Rz = Rot(az, 'z');

    R = Rx*Ry*Rz;
else
    'Unknown convention'
end

repo: ARTE