text:
function polycz() % autogenerated function wrapper
    % crosssel.m                      Alexander Allmann
    % function to select earthquakes in a cross-section and make them the
    % current catalog in the main map windo
    %
    %
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    global newa

    report_this_filefun(mfilename('fullpath'));

    figNumber=findobj('Type','Figure','-and','Name','Z-Value-Cross-section');

    figure(figNumber);

    ax=findobj(gcf,'Tag','mainmap_ax');
    [x,y, mouse_points_overlay] = select_polygon(ax);

    hold on
    plot(x,y,'b-');
    YI = -newa.Depth;          % this substitution just to make equation below simple
    XI = newa(:,end);
    lb = polygon_filter(x,y, XI, YI, 'inside');

    %plot the selected eqs and mag freq curve
    ZG.newt2 = newa.subset(lb);
    newcat = newa.subset(lb);
    pl = plot(ZG.newt2(:,length(ZG.newt2.subset(1))),-ZG.newt2.Depth,'xk'); %FIXME huh? what are they trying to do? looks like indexing are mixed up
    set(pl,'MarkerSize',8, 'LineWidth',1)
    timeplot()
end

author: CelsoReyes


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%
classdef MultiStageProblem < handle

  properties
    solver
    stageList
  end

  methods

    function self = MultiStageProblem(varargin)
      % ocl.Solver(stages, transitions, casadi_options)

      ocl.utils.checkStartup()

      % ocl.Solver(stages, transitions, opt)
      p = ocl.utils.ArgumentParser;

      p.addKeyword('stages', {}, @(el) iscell(el) || isa(el, 'ocl.Stage'));
      p.addKeyword('transitions', {}, @(el) iscell(el) || ishandle(el) );

      p.addParameter('nlp_casadi_mx', false, @islogical);
      p.addParameter('controls_regularization', true, @islogical);
      p.addParameter('controls_regularization_value', 1e-6, @isnumeric);

      p.addParameter('casadi_options', ocl.casadi.CasadiOptions(), @(el) isstruct(el));
      p.addParameter('verbose', true, @islogical);
      p.addParameter('userdata', [], @(in)true);
      p.addParameter('transition_type', 1, @isnumeric);

      r = p.parse(varargin{:});

      stageList = r.stages;
      transitionList = r.transitions;

      nlp_casadi_mx = r.nlp_casadi_mx;
      controls_regularization = r.controls_regularization;
      controls_regularization_value = r.controls_regularization_value;

      casadi_options = r.casadi_options;
      verbose = r.verbose;
      userdata = r.userdata;
      transition_type = r.transition_type;

      solver = ocl.casadi.CasadiSolver(stageList, transitionList, ...
                                       nlp_casadi_mx, controls_regularization, ...
                                       controls_regularization_value, casadi_options, ...
                                       verbose, userdata, transition_type);

      % set instance variables
      self.stageList = stageList;
      self.solver = solver;
    end

    function [sol,times,info] = solve(self, ig)
      % [sol, times] = solve()
      % [sol, times] = solve(ig)

      s = self.solver;

      if nargin==1
        % ig InitialGuess
        ig = self.solver.getInitialGuessWithUserData();
      end

      [sol,times,solver_info] = s.solve(ig);

      if nargout >=3
        info = solver_info;
      end

      ocl.utils.warningNotice()
    end

    function r = timeMeasures(self)
      r = self.solver.timeMeasures;
    end

    function ig = ig(self)
      ig = self.getInitialGuess();
    end

    function igList = getInitialGuess(self)
      igList = self.solver.getInitialGuess();
    end
  end
end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function [x] = parse(inStr)
% parse is a function which takes in a string vector of blank separated text
% and parses out the individual string items into a n item matrix, one row
% for each string.
%
% function [x] = parse(inStr)
% x     - the return matrix of strings
% inStr - the blank separated string vector

% Binary and Real-Valued Simulation Evolution for Matlab
% Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% 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 1, 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. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

sz=size(inStr);
strLen=sz(2);
x=blanks(strLen);
wordCount=1;
last=0;
for i=1:strLen,
  if inStr(i) == ' '
    wordCount = wordCount + 1;
    x(wordCount,:)=blanks(strLen);
    last=i;
  else
    x(wordCount,i-last)=inStr(i);
  end
end

author: HuangCongQing


text:
%kvnthin 'N-1 compositions of Thinnings'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vnthin.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'input image'
% InputFile: str1 'Struct Element 1', required: 'Structuring Element'
% InputFile: str2 'Struct Element 2', required: 'Structuring Element'
% Integer: n 'Number of Thinnings', default: 5: 'number of thinnings'
% MultiChoice: r 'Struct El. Rotation Step', default: 1: 'Rotation step'
%    Choices are:
%   1: '1 '
%   2: '2 '
%   3: '4 '
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvnthin({i, str1, str2}, {'i','';'str1','';'str2','';'n',5;'r',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vnthin - N-1 compositions of Thinnings
%
%  DESCRIPTION
% .I vnthin
% will do N successive thinnings of an image.
%
% The output image data type will be the same as that of the highest data type
% of the input image.
%
% Some useful pairs of structuring elements:
%
%                         0 0 0    1 1 1
%   Homotopic Skeleton    0 1 0    0 0 0
%                         1 1 1    0 0 0
%
%                         0 0 0    1 0 0
%   Homotopic Marking     0 1 1    1 0 0
%                         0 0 0    1 0 0
%
%                         0 0 0    0 0 0
%   Skeleton Pruning      0 1 0    1 0 1
%                         0 0 0    1 1 1
%
%                         1 1 1    0 0 0
%   Boundary              1 1 1    0 0 0
%                         1 1 1    0 0 0
%
%
%
%  EXAMPLES
% vnthin -i feath.xv -str1 b1.str -str2 b2.str -n 5 -r 1 -o outimage.xv
%
% Will do 5 successive thinnings (vnthin) of the input image "feath.xv" by the structuring element b1.str and b2.str, rotating the struct. els. 1 step after each thinning. The resulting image will be written to  "outimage.xv".
%
%  "SEE ALSO"
% vcast1(1)
%
%  RESTRICTIONS
% .I vnthin
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [no,el]=fillsurf(node,face)
%
% [no,el]=fillsurf(node,face)
%
% calculate the enclosed volume for a closed surface
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%    node:  node coordinates
%    face:  surface triangle list
%
% output:
%    vol:   total volume of the enclosed space
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

[no,el]=surf2mesh(node,face,[],[],1,1,[],[],0,'tetgen','-YY');

author: fangq


text:
function tests = test_checkcode
% Test for possible problems in all of MATLAB code files
%__________________________________________________________________________
% Copyright (C) 2017 Wellcome Trust Centre for Neuroimaging

% $Id: test_checkcode.m 7111 2017-06-16 09:01:09Z guillaume $

tests = functiontests(localfunctions);

function test_checkcode_all(testCase)
testCase.assumeTrue(strcmp(spm_check_version,'matlab'));

f = cellstr(spm_select('FPlistRec',spm('Dir'),'\.m$')); % discard external?
s = checkcode(f,'-id');
checked = false(size(s));
for i=1:numel(s)
    for j=1:numel(s{i})
        if isequal(s{i}(j).id,'SYNER') % Syntax error
            checked(i) = true;
        end
    end
end
testCase.verifyFalse(any(checked));
if any(checked)
    fprintf('Problematic files are:\n');
    disp(f(checked));
end

function t_e_s_t_analyzeCodeCompatibility(testCase) % disabled for now
testCase.assumeFalse(~strcmp(spm_check_version,'matlab') | ...
    spm_check_version('matlab','9.3') < 0);

import matlab.unittest.constraints.*
r = analyzeCodeCompatibility(spm('Dir'),'IncludeSubfolders',true); % discard external?
idx = find(r.Recommendations.Severity == 'Error');
testCase.verifyThat(numel(idx), IsEqualTo(0));
for i=1:numel(idx)
    fprintf('%s:\n',r.Recommendations.File(idx(i)));
    fprintf('  %s\n',r.Recommendations.Description(idx(i)));
end

author: spm


text:
function tests = test_ehsan_yale_faces
  tests = functiontests(localfunctions);
end

function setupOnce(testCase)
    close all;
end

function teardownOnce(testCase)
    close all;
end

function test_basic(testCase)
    yf = spx.data.image.EhsanYaleFaces();
    verifyEqual(testCase, yf.image_height, 48);
    verifyEqual(testCase, yf.image_width, 42);
    verifyEqual(testCase, yf.num_subjects, 38);
    verifyEqual(testCase, yf.num_faces_per_subject, 64);
    verifyEqual(testCase, yf.num_total_faces, 64*38);
end

function test1(testCase)
    yf = spx.data.image.EhsanYaleFaces();
    mf = spx.graphics.Figures();
    mf.new_figure('Example faces');
    canvas = yf.create_random_canvas(5, 5);
    imshow(canvas);
    colormap(gray);
    axis image;
    axis off;
    pause(.2);
end

author: indigits

repo: sparse-plex


text:
function C = cs_symperm (A,p)                                               %#ok
%CS_SYMPERM symmetric permutation of a symmetric matrix.
%   C = cs_symperm(A,p) computes C = A(p,p), but accesses only the
%   upper triangular part of A, and returns C upper triangular (A and C are
%   symmetric with just their upper triangular parts stored).  A must be square.
%
%   Example:
%       Prob = UFget ('HB/bcsstk01') ; A = Prob.A ;
%       p = cs_amd (A) ;
%       C = cs_symperm (A, p) ;
%       cspy (A (p,p)) ;
%       cspy (C) ;
%       C - triu (A (p,p))
%
%   See also CS_PERMUTE, SUBSREF, TRIU.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_symperm mexFunction not found') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [] = advance_random()
%   Creates the next batch of 55 random numbers.

global oldrand ;

for j1 = 0:23
    new_random = oldrand(j1 + 1) - oldrand(j1 + 31 + 1);
    if (new_random < 0.0)
        new_random = new_random + 1.0 ;
    end
    oldrand(j1 + 1) = new_random ;
end

for j1 = 24:54
    new_random = oldrand(j1 + 1) - oldrand((j1 - 24) + 1);
    if (new_random < 0.0)
        new_random = new_random + 1.0 ;
    end
    oldrand(j1 + 1) = new_random ;
end
end

author: chudur-budur


text:
function [caNodeIndices, vResolution_] = ex_CreateIndexCatalog(mCatalog, mPolygon, bMap, nGriddingMode, nNumberEvents, fRadius, fSizeRectHorizontal, fSizeRectDepth)
% function [caNodeIndices] = ex_CreateIndexCatalog(mCatalog, mPolygon, bMap, nGriddingMode,
%                                                  nNumberEvents, fRadius, fSizeRectHorizontal, fSizeRectDepth)
% -------------------------------------------------------------------------------------------------------------
% Creates a cell-array with subcatalogs for every grid node defined by mPolygon. These subcatalogs
%   contain only indices to the earthquake "rows" in mCatalog.
%
% Input parameters:
%   mCatalog              Earthquake catalog
%   mPolygon              Polygon (defined by ex_selectgrid)
%   bMap                  Calculate cell-array for a map (=1) or a cross-section (=0)
%   nGriddingMode         Mode of creating grid node subcatalogs
%                         0: Constant number of events
%                         1: Constant radius
%                         2: Rectangular grid node samples
%   nNumberEvents         Number of events per grid node (nGriddingMode == 0)
%   fRadius               Radius of grid node sample (nGriddingMode == 1)
%   fSizeRectHorizontal   Latitude/horizontal size of rectangle (nGriddingMode == 2)
%   fSizeRectDepth        Longitude/depth size of rectangle (nGriddingMode == 2)
%
% Output parameters:
%   caNodeIndices         Cell-array with index-catalogs per grid node of mPolygon
%
% Danijel Schorlemmer
% June 17, 2002

global bDebug;
if bDebug
  report_this_filefun(mfilename('fullpath'));
end

% Create the catalogs for each node with pointers to the overall catalog
nNumberNodes_ = length(mPolygon(:,1));
caNodeIndices = cell(nNumberNodes_, 1);

% If cross-section calculate the length along cross-section
if ~bMap
  [nRow_, nColumn_] = size(mCatalog);
  vXSecX_ = mCatalog(:,nColumn_);  % length along x-section
  vXSecY_ = (-1) * mCatalog(:,7);  % depth of hypocenters
end

% vResolution give the radius (for nGriddingMode = 0) and no. of events
% (for nGriddingMode = 1)
vResolution_(:,1)=ones(nNumberNodes_,1)*nan;

% Loop over all points of the polygon
for nNode_ = 1:nNumberNodes_
  % Get the grid node coordinates
  fX_ = mPolygon(nNode_, 1);
  fY_ = mPolygon(nNode_, 2);

  if (nGriddingMode == 0) | (nGriddingMode == 1)  % Fixed radius or fixed number
    % Calculate distance from center point
    if bMap
      vDistances_ = sqrt(((mCatalog(:,1)-fX_)*cos(pi/180*fY_)*111).^2 + ((mCatalog(:,2)-fY_)*111).^2);
    else
      vDistances_ = sqrt(((vXSecX_ - fX_)).^2 + ((vXSecY_ - fY_)).^2);
    end
    if nGriddingMode == 0 % Fixed number
      if length(mCatalog(:,1)) == 0
        caNodeIndices{nNode_} = [];
        % NaN for no events
        vResolution_(nNode_) = nan;
      elseif nNumberEvents > length(mCatalog(:,1))
        caNodeIndices{nNode_} = vIndices(1:length(mCatalog(:,1)));
        % take the maximal distance for all eq. in the catalog
        vResolution_(nNode_) = max(vdistances_);
      else
        % Use first nNumberEvents events
        [vTmp, vIndices] = sort(vDistances_);
        caNodeIndices{nNode_} = vIndices(1:nNumberEvents);
        % radius of the nNumberEvents-th event in the sorted vDistances_
        vResolution_(nNode_) = vTmp(nNumberEvents);
      end
    else % Fixed radius
      % Use all events within fRadius
      caNodeIndices{nNode_} = find(vDistances_ <= fRadius);
      vResolution_(nNode_) = length(find(vDistances_ <= fRadius));
    end
  else % Rectangular gridding (nGriddingMode == 2)
    if bMap
      vSel_ = ((mCatalog(:,1) >= (fX_ - fSizeRectHorizontal/2)) & (mCatalog(:,1) < (fX_ + fSizeRectHorizontal/2)) & ...
        (mCatalog(:,2) >= (fY_ - fSizeRectDepth/2)) & (mCatalog(:,2) < (fY_ + fSizeRectDepth/2)));
      vResolution_(nNode_) = length(find(vSel_ > 0))
    else
      vSel_ = ((vXSecX_ >= (fX_ - fSizeRectHorizontal/2)) & (vXSecX_ < (fX_ + fSizeRectHorizontal/2)) & ...
        (vXSecY_ >= (fY_ - fSizRectDepth/2)) & (vXSecY_ < (fY_ + fSizeRectDepth/2)));
      vResolution_(nNode_) = length(find(vSel_ > 0))
    end
    caNodeIndices{nNode_} = find(vSel_ == 1);
  end
end; % of for nNode_


text:
function quad_error = simplex_unit_monomial_quadrature ( dim_num, expon, ...
  point_num, x, w )

%*****************************************************************************80
%
%% SIMPLEX_UNIT_MONOMIAL_QUADRATURE: quadrature of monomials in a unit simplex.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer EXPON(DIM_NUM), the exponents.
%
%    Input, integer POINT_NUM, the number of points in the rule.
%
%    Input, real X(DIM_NUM,POINT_NUM), the quadrature points.
%
%    Input, real W(POINT_NUM), the quadrature weights.
%
%    Output, real QUAD_ERROR, the quadrature error.
%

%
%  Get the exact value of the integral of the unscaled monomial.
%
  scale = simplex_unit_monomial_int ( dim_num, expon );
%
%  Evaluate the monomial at the quadrature points.
%
  value = monomial_value ( dim_num, point_num, x, expon );
%
%  Compute the weighted sum and divide by the exact value.
%
  volume = simplex_unit_volume ( dim_num );
  quad = volume * ( w * value' ) / scale;
%
%  Error:
%
  exact = 1.0;
  quad_error = abs ( quad - exact );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
 function ir_stem_spectrum(x, y, varargin)
%function ir_stem_spectrum(x, y, varargin)
%|
%| default spectrum style for engn 100

arg.thresh = 0;
arg.color = [0 0.5 0];
arg = vararg_pair(arg, varargin);

k = y > max(y(:)) * arg.thresh;

stem(x(k), y(k), ...
	'color', arg.color, ...
	'linewidth', 3.0, ...
	'marker', 'none')
ylabelf('amplitude')

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


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

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * compression, tension, compressive, tensile
% * displacement control, displacement boundary condition
% * hexahedral elements, hex8
% * cube, box, rectangular
% * static, solid
% * multigeneration, multi-generation
% * multi-step analysis
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

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

%% Control parameters

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

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

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

%Define applied displacement
appliedStrain=0.3; %Linear strain (Only used to compute applied stretch)
loadingOption='tension'; % or 'tension'
switch loadingOption
    case 'compression'
        stretchLoad=1-appliedStrain; %The applied stretch for uniaxial loading
    case 'tension'
        stretchLoad=1+appliedStrain; %The applied stretch for uniaxial loading
end
displacementMagnitude=(stretchLoad*sampleHeight)-sampleHeight; %The displacement magnitude

%Material parameters
k_factor=50;
c1=2;
m1=2;
k=c1*k_factor;
c1_g=[c1/1000 c1*2];
k_g=c1_g*k_factor;

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

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

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

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

%% Splitting mesh into 2 material groups

X=V(:,1); Y=V(:,2); Z=V(:,3);
VE=[mean(X(E),2) mean(Y(E),2) mean(Z(E),2)];

logicMaterial_1=VE(:,1)<0;

elementMaterialID=logicMaterial_1+1;

%Reoder E to cope with FEBio bug in relation to element ordering and
%multiple material sets
E=[E(elementMaterialID==1,:); E(elementMaterialID==2,:);];
elementMaterialID=[elementMaterialID(elementMaterialID==1,:); elementMaterialID(elementMaterialID==2,:);];

%Fix meshStruct to allow for meshView based visualization
meshStruct.elementMaterialID=elementMaterialID;
meshStruct.elements=E;

%%
% Plotting model boundary surfaces and a cut view

hFig=cFigure;

hs=subplot(1,2,1); hold on;
title('Model boundary surfaces and labels','FontSize',fontSize);
gpatch(Fb,V,Cb,'k',faceAlpha1); view(3)
colormap(gjet(12)); icolorbar;
axisGeom(gca,fontSize);

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

drawnow;

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

%Define supported node sets
logicFace=Cb==5; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList=unique(Fr(:)); %Node set part of selected face

%Prescribed displacement nodes
logicPrescribe=Cb==6; %Logic for current face set
Fr=Fb(logicPrescribe,:); %The current face set
bcPrescribeList=unique(Fr(:)); %Node set part of selected face

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

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

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

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

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

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

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

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

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

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

%Create control structure for use by all steps
stepStruct.Control.analysis='STATIC';
stepStruct.Control.time_steps=numTimeSteps;
stepStruct.Control.step_size=1/numTimeSteps;
stepStruct.Control.solver.max_refs=max_refs;
stepStruct.Control.solver.max_ups=max_ups;
stepStruct.Control.time_stepper.dtmin=dtmin;
stepStruct.Control.time_stepper.dtmax=dtmax;
stepStruct.Control.time_stepper.max_retries=max_retries;
stepStruct.Control.time_stepper.opt_iter=opt_iter;

%Add template based default settings to proposed control section
[stepStruct.Control]=structComplete(stepStruct.Control,febio_spec.Control,1); %Complement provided with default if missing

%Remove control field (part of template) since step specific control sections are used
febio_spec=rmfield(febio_spec,'Control');

febio_spec.Step.step{1}.Control=stepStruct.Control;
febio_spec.Step.step{1}.ATTR.id=1;
febio_spec.Step.step{2}.Control=stepStruct.Control;
febio_spec.Step.step{2}.ATTR.id=2;

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

materialName2='Multigen_material';
febio_spec.Material.material{2}.ATTR.name=materialName2;
febio_spec.Material.material{2}.ATTR.id=2;
febio_spec.Material.material{2}.ATTR.type='multigeneration';

febio_spec.Material.material{2}.generation{1}.ATTR.id=1;
febio_spec.Material.material{2}.generation{1}.start_time=0;
febio_spec.Material.material{2}.generation{1}.solid{1}.ATTR.type='Ogden unconstrained';
febio_spec.Material.material{2}.generation{1}.solid{1}.c1=c1_g(1);
febio_spec.Material.material{2}.generation{1}.solid{1}.m1=m1;
febio_spec.Material.material{2}.generation{1}.solid{1}.c2=c1_g(1);
febio_spec.Material.material{2}.generation{1}.solid{1}.m2=-m1;
febio_spec.Material.material{2}.generation{1}.solid{1}.cp=k_g(1);

febio_spec.Material.material{2}.generation{2}.ATTR.id=2;
febio_spec.Material.material{2}.generation{2}.start_time=1;
febio_spec.Material.material{2}.generation{2}.solid{1}.ATTR.type='Ogden unconstrained';
febio_spec.Material.material{2}.generation{2}.solid{1}.c1=c1_g(2);
febio_spec.Material.material{2}.generation{2}.solid{1}.m1=m1;
febio_spec.Material.material{2}.generation{2}.solid{1}.c2=c1_g(2);
febio_spec.Material.material{2}.generation{2}.solid{1}.m2=-m1;
febio_spec.Material.material{2}.generation{2}.solid{1}.cp=k_g(2);

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

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

partName2='Part2';
febio_spec.Mesh.Elements{2}.ATTR.name=partName2; %Name of this part
febio_spec.Mesh.Elements{2}.ATTR.type='hex8'; %Element type
febio_spec.Mesh.Elements{2}.elem.ATTR.id=(1+nnz(elementMaterialID==1):1:nnz(elementMaterialID==1)+nnz(elementMaterialID==2))'; %Element id's
febio_spec.Mesh.Elements{2}.elem.VAL=E(elementMaterialID==2,:);

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

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

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

febio_spec.MeshDomains.SolidDomain{2}.ATTR.name=partName2;
febio_spec.MeshDomains.SolidDomain{2}.ATTR.mat=materialName2;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x,y,z';

%STEP 1 Tension
febio_spec.Step.step{1}.Boundary.bc{1}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{1}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{1}.Boundary.bc{1}.dof='z';
febio_spec.Step.step{1}.Boundary.bc{1}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{1}.scale.VAL=displacementMagnitude;
febio_spec.Step.step{1}.Boundary.bc{1}.relative=1;

febio_spec.Step.step{1}.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Step.step{1}.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{1}.Boundary.bc{2}.dofs='x,y';

%STEP 2 Return form tension
febio_spec.Step.step{2}.Boundary.bc{1}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{1}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{2}.Boundary.bc{1}.dof='z';
febio_spec.Step.step{2}.Boundary.bc{1}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{1}.scale.VAL=-displacementMagnitude;
febio_spec.Step.step{2}.Boundary.bc{1}.relative=1;

febio_spec.Step.step{2}.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Step.step{2}.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{2}.Boundary.bc{2}.dofs='x,y';

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

febio_spec.LoadData.load_controller{2}.ATTR.id=2;
febio_spec.LoadData.load_controller{2}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{2}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{2}.points.point.VAL=[1 0; 2 1];

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

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_strainEnergy;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sed';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.element_data{1}.VAL=1:size(E,1);

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

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

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

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

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

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

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

%% Import FEBio results

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

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

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

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

    %%
    % Importing element stress from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_strainEnergy),1,1);

    %Access data
    E_energy=dataStruct.data;

    [F,CF]=element2patch(E,E_energy(:,:,1));

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

    c1_plot=c1*ones(size(timeVec));
    cg_plot=c1_g(1)*ones(size(timeVec));
    cg_plot(timeVec>=1)=c1_g(2);

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);

    subplot(1,2,1); hold on;
    title('Ogden parameter c_1');
    xlabel('Time'); ylabel('c_1');
    plot(timeVec,c1_plot,'b-','lineWidth',2);
    plot(timeVec,cg_plot,'r-','lineWidth',2);
    hp1=plot(timeVec(1),c1_plot(1),'b.','MarkerSize',50);
    hp2=plot(timeVec(1),cg_plot(1),'r.','MarkerSize',50);
    legend([hp1 hp2],'Material 1','Material 2');
    axis tight; axis square; set(gca,'fontsize',fontSize);
    grid on;

    subplot(1,2,2); hold on;
    hp3=gpatch(F,V_DEF(:,:,end),CF,'k',1); %Add graphics object to animate
    gpatch(Fb,V,0.5*ones(1,3),'k',0.25); %A static graphics object

    colormap(gjet(250)); hc=colorbar;
    caxis([0 max(E_energy(:))]);
    axisGeom(gca,fontSize);
    axis(axisLim(V_DEF)); %Set axis limits statically
    axis manual;
    camlight headlight;
    drawnow;

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

        [~,CF]=element2patch(E,E_energy(:,:,qt));

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp3 hp3 hp1 hp1 hp2 hp2]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData','XData','YData','XData','YData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),CF,timeVec(qt),c1_plot(qt),timeVec(qt),cg_plot(qt)}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

end

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

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function exact = p12_exact ( )

%*****************************************************************************80
%
%% P12_EXACT returns the estimated integral for problem 12.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the estimated value of the integral.
%
  exact = 0.5;

  return
end


text:
function Y=tensor2Mat(X,param2,C)
%%TENSOR2MAT Unfold a (real or complex) tensor into a matrix. This can
%            either be done as an n-way matricization using a standard
%            ordering of the modes, or one can explicitly specify the
%            ordering of the modes. The tensor is just a hypermatrix. The
%            unfolding produces a 2D matrix as the output. The
%            transformation of tensors into matrices plays a role in
%            numerous multilinear algorithms, such as the multilinear
%            singular value decomposition.
%
%INPUTS: X The hypermatrix (tensor) that is to be unfolded. This should
%          have 2 or more modes. For example, a 4-mode hyper matrix would
%          be addressable as X(i1,i2,i3,i4). The matrix can be real or
%          complex.
% param2, C If the algorithm is called with only two parameters, then
%          param2=n is the number n of the mode of the unfolding and an
%          n-mode matricization is performed. The mode number ranges
%          from 1 through N --the number of modes of the matrix X (the
%          number of indices needed to address an element in X). The
%          permutation of the modes to the columns is as in [1], which is
%          [1:1:(n-1),(n+1):1:N] and the number of rows of the matrix Y
%          returned is equal to the number of elements in the nth way. On
%          the other hand, some authors, such as in [2], use a different
%          permutation. In such an instance, the explicit ordering can be
%          specified using param2=R and also providing C as in [3]. R
%          contains the indices of the modes that are to be mapped to rows
%          and C contains the indices of modes that are to be mapped to
%          columns. For example, for a 4-mode matrix one might use R=[2;3],
%          C=[4;1], which makes the dimensionality of the rows of the
%          output equal to the sum of the dimensionalities of
%          modes 2 and 3.
%
%OUTPUT: Y A 2D matrix containing the elements of the hypermatrix X
%          arranged according to the inputs.
%
%The algorithm can be called as
%Y=tensor2Mat(X,n);
%for a standard n-mode matricization or as
%Y=tensor2Mat(X,R,C);
%for general ordering to be used. The standard n-mode matricization is
%consistent with the orderings used in the nModeProd function.
%
%As an example, consider the 3X2X3 hypermatrix
% A=zeros(3,2,3);
% A(1,1,1)=1;
% A(1,1,2)=1;
% A(2,1,1)=1;
% A(2,1,2)=-1;
% A(2,1,3)=2;
% A(3,1,1)=2;
% A(3,1,3)=2;
% A(1,2,1)=2;
% A(1,2,2)=2;
% A(2,2,1)=2;
% A(2,2,2)=-2;
% A(2,2,3)=4;
% A(3,2,1)=4;
% A(3,2,3)=4;
% %A standard 1-mode matrix unfolding is
% Y=tensor2Mat(A,1)
%Providing the result
% Y=[1     2     1     2     0     0;
%    1     2    -1    -2     2     4;
%    2     4     0     0     2     4];
%On the other hand, the example in [2] gives a different answer, because
%they use a different ordering. That is, they rearrange, the columns in the
%unfolding differently. Here, one can get the same answer using
% Y=tensor2Mat(A,1,[3,2])
%which returns
% Y=[1     1     0     2     2     0;
%    1    -1     2     2    -2     4;
%    2     0     2     4     0     4];
%
%Introductions to tensor operations are in [3] and [4].
%
%REFERENCES:
%[1] J. Salmi, A. Richter, and V. Koivunen, "Sequential unfolding SVD for
%    tensors with applications in array signal processing," IEEE
%    Transactions on Signal Processing, vol. 57, no. 12, pp. 4719-4733,
%    Dec. 2009.
%[2] L. de Lathauwer, B. de Moore, and J. Vandewalle, "A multilinear
%    singular value decomposition," SIAM Journal on Matrix Analysis and
%    Applications, vol. 21, no. 4, pp. 1253-1278, 2000.
%[3] T. G. Kolda, "Multilinear operators for higher-order decompositions,"
%    Sandia National Laboratories, Tech. Rep. SAND2006-2081, Apr. 2006.
%    [Online]. Available: http://www.sandia.gov/~tgkolda/pubs/pubfiles/SAND2006-2081.pdf
%[4] R. G. Kolda and B. W. Bader, "Tensor decompositions and applications,"
%    SIAM Review, vol. 51, no. 3, pp. 455-500, 2009.
%
%June 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

I=size(X);%The dimensionality of each way of X.

%If a standard n-mode matricization is being performed
if(nargin<3)
    d=length(I);
    C=[1:1:(param2-1),(param2+1):1:d];
end

R=param2;
J=prod(I(R));
K=prod(I(C));
Y=reshape(permute(X,[R(:);C(:)]),J,K);%Convert X to the matrix Y.

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function x = cs_pvec (b,p)                                                  %#ok
%CS_PVEC x=b(p)
%
% Example:
%   x = cs_pvec (b,p)
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_pvec mexFunction not found') ;


text:
%% DEMOLONG  Short demonstration of long numbers
%

%%
% The purpose of the long toolbox was to compute rigorous bounds for the
% value of certain standard functions. Those values are needed to initialize
% the INTLAB system. The long toolbox is slow, but fast enough to do the job.

setround(0)                 % set rounding to nearest
longprecision(0);           % default option

%% Definition of long numbers
%
% Long numbers are stored in midpoint radius representation.
% The midpoint is stored in an array of precision to be specified,
% the error is stored in one number. Long numbers or vectors are generated
% by the constructor long:

x = long(7)
V = long([-1;3])

%%
% Long vectors are always column vectors (or forced to column vectors).

%% Conversion
% Conversion of double to long by the constructor "long" is always into
% the current internal precision.
%
% The internal precision may be specified by "longprecision".
% The call without input parameter gives the current working precision
% in decimals, the call with input parameter sets working precision.

p = longprecision
longprecision(50)

%%
% Now, the statement "x = long(7)" generates a long number with
% approximately 50 decimal digits. This is only approximate because
% internal representation is to some base beta, a power of 2.

%% Output of long numbers
% Output is usually to a little more than double precision. If you want
% to see more digits, say k, use "display" with second parameter equal
% to k. To see all digits, use k=0.

longprecision
x = 1/long(7)
display(x,40)
display(x,0)

%%
% Output of long numbers is not rigorous. All but a few of the last
% digits are correct.

%% Arithmetic operations
% Long operations +,-,*,/ and ^ are supported. Note that operations on
% vectors are always performed elementwise.

x = [ long(3) ; -7 ]
x*x

%% Output of long intervals
% The display routine takes uncertainties into account. Only the
% correct digits plus some extra are displayed.

longprecision(50);
x = long(1)/37;
display(x,0)
for i=1:100
  x=x*x; x=x*37;
end
display(x,0)

%% Interval and non-interval operations
% Computing with uncertainties may be switched off by

longinit('WithoutErrorTerm')
longprecision(50);
x = long(1)/37;
display(x,0)
for i=1:100
  x=x*x; x=x*37;
end
display(x,0)

%%
% In this case all digits including incorrect ones are displayed.
% Computing without error term is a usual long precision arithmetic
% with specified precision. Note that scalar operations suffer from
% quite some interpretation overhead.

%% Conversion between long and double
% Conversion from long to double is approximately to nearest, conversion
% to interval is rigorous.
%
% For example, in the following the function "longpi" calculates
% "pi" to the specified longprecision, "IntPi" is a true inclusion of the
% transcendental number "pi".

longinit('WithErrorTerm');
longprecision(100);
Pi = longpi;
display(Pi,0)
flptPi = long2dble(Pi)
IntPi = long2intval(Pi)
format long
infsup(IntPi)

%% Long numbers with error term
% Long numbers may be specified with an explicit error term.
% For example,

longprecision(50);
x = long(-1.5)
display(x,0)
x = addlongerror(x,1e-40)
display(x,0)

%%
% defines x to be an interval with midpoint -1.5 and radius
% approximately 10^(-40). Only meaningful digits are stored and displayed.

%% Specifying extremely small errors
% For very small errors leaving the range double precision
% floating point numbers, the error may be specified by
% the mantissa and the exponent of the error:

longprecision(50);
x = long(2^-1000)^2;
x = addlongerror(x,1,-620)

%%
% The final x, which is 2^(-2000), is afflicted with an error
% of 10^(-620).

%% Taylor series: an example
% As an example, the following code computes the value of E = exp(x)
% by a Taylor series:

p = 100; longprecision(p);
x = -30;
t = 1; T = long(1); E = T; k = 0;
while abs(t)>10^(-p)
  k = k+1;
  t = t*x/k;
  T = T*x/k;
  E = E + T;
end
k
exp(x)
display(E,0)

%%
% Note that for large negative values of x there quite some
% cancellation. This can be seen by

x = 30;
t = 1; T = long(1); E = T; k = 0;
while abs(t)>10^(-p)
  k = k+1;
  t = t*x/k;
  T = T*x/k;
  E = E + T;
end
k
1/exp(x)
display(1/E,0)

%% Ill-conditioned polynomials
% Consider the following polynomial:

P = inline(' 4999*x.^6 - 200*x.^5 + 102*x.^4 - 2*x.^3 - 2500*x.^2 + 100*x - 1 ')

%%
% This is an example Bugeaud-Mignotte polynomial. The general form is
%
% ( X^n - aX + 1 )^k - 2X^(nk-k)(aX-1)^k
%
% where a>=10, n>=3 and k>=2.
%
% Those polynomials are constructed to have a pair of very close real roots near c=1/a+1/a^(n+1).
% A graph near c looks as follows:

e = 3e-8;
c = 1/50+1/50^4;
x = c*(1+linspace(-e,e));
close
plot(x,P(x),x,0*x)

%%
% From the graph it is not clear whether the polynomial has no, a double or two real roots
% in the interval c*[1-e,1+e]. An evaluation using the long package yields
% the following:

y = long2dble(P(long(x)));
close
plot(x,y,x,0*x)

%% Sample programs
% For sample programs using long numbers, see for example the
% source codes of long\longpi.m or long\@long\exp.m  .

%% Enjoy INTLAB
% INTLAB was designed and written by S.M. Rump, head of the Institute for Reliable Computing,
% Hamburg University of Technology. Suggestions are always welcome to rump (at) tuhh.de

repo: OpenMAS


text:
function inspect_issue509

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY

%%
% all defaults, should resemble design from Stefan

cfg = [];
data = ft_steadystatesimulation(cfg);

cfg = [];
cfg.viewmode = 'vertical';
ft_databrowser(cfg, data);

%%
% the following is an attempt to resemble the design from Cassia

cfg = [];
cfg.ntrials = 120;
cfg.level1.condition  = [1 2]; % attend low/high
cfg.level1.gain       = [1.1 1.2];
cfg.level2.condition  = [1 2]; % sequence goes up/down
cfg.level2.gain       = [1.1 1.2];
cfg.level3.condition  = [1 2 3]; % note 1 2 3
cfg.level3.gain       = [1.1 1.2 1.3];
cfg.stimulus1.mode = 'periodic';
cfg.stimulus1.isi = 1/20;
cfg.baseline = 0.1;
cfg.stimulus1.onset = 0;
cfg.stimulus1.onsetjitter = 0;
cfg.stimulus2.mode = 'off';

data = ft_steadystatesimulation(cfg);

cfg = [];
cfg.viewmode = 'vertical';
cfg.channel = [2 3];
ft_databrowser(cfg, data);


text:
function [bor,F0,F1] = spm_BMS_bor(L,posterior,priors,C)
% Compute Bayes Omnibus Risk
% FORMAT [bor,F0,F1] = spm_BMS_bor(L,posterior,priors,C)
%
% L         Log model evidence table (models x  subjects)
% posterior .a model counts, .r model-subject probs
% priors    .a model counts
% C         if this field is specified then BOR under family prior
%           is computed, otherwise BOR under model prior is computed.
%           C(k,f) = 1 if model k belongs to family f (0 otherwise)
%
% REFERENCES:
%
% Rigoux, L, Stephan, KE, Friston, KJ and Daunizeau, J. (2014)
% Bayesian model selection for group studies - Revisited.
% NeuroImage 84:971-85. doi: 10.1016/j.neuroimage.2013.08.065
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_BMS_bor.m 6444 2015-05-21 11:15:48Z guillaume $

if nargin < 4
    options.families = 0;
    % Evidence of null (equal model freqs)
    F0 = FE_null(L,options);
else
    options.families = 1;
    options.C = C;
    % Evidence of null (equal model freqs) under family prior
    [tmp,F0] = FE_null(L,options);
end

% Evidence of alternative
F1 = FE(L,posterior,priors);

% Implied by Eq 5 (see also p39) in Rigoux et al.
% See also, last equation in Appendix 2
bor = 1/(1+exp(F1-F0));

function [F,ELJ,Sqf,Sqm] = FE(L,posterior,priors)
% derives the free energy for the current approximate posterior
% This routine has been copied from the VBA_groupBMC function
% of the VBA toolbox http://code.google.com/p/mbb-vb-toolbox/
% and was written by Lionel Rigoux and J. Daunizeau
%
% See equation A.20 in Rigoux et al. (should be F1 on LHS)

[K,n] = size(L);
a0 = sum(posterior.a);
Elogr = psi(posterior.a) - psi(sum(posterior.a));
Sqf = sum(gammaln(posterior.a)) - gammaln(a0) - sum((posterior.a-1).*Elogr);
Sqm = 0;
for i=1:n
    Sqm = Sqm - sum(posterior.r(:,i).*log(posterior.r(:,i)+eps));
end
ELJ = gammaln(sum(priors.a)) - sum(gammaln(priors.a)) + sum((priors.a-1).*Elogr);
for i=1:n
    for k=1:K
        ELJ = ELJ + posterior.r(k,i).*(Elogr(k)+L(k,i));
    end
end
F = ELJ + Sqf + Sqm;

function [F0m,F0f] = FE_null (L,options)
% Free energy of the 'null' (H0: equal frequencies)
%
% F0m       Evidence for null (ie. equal probs) over models
% F0f       Evidence for null (ie. equal probs) over families
%
% This routine derives from the VBA_groupBMC function
% of the VBA toolbox http://code.google.com/p/mbb-vb-toolbox/
% written by Lionel Rigoux and J. Daunizeau
%
% See Equation A.17 in Rigoux et al.

[K,n] = size(L);
if options.families
    f0 = options.C*sum(options.C,1)'.^-1/size(options.C,2);
    F0f = 0;
else
    F0f = [];
end
F0m = 0;
for i=1:n
    tmp = L(:,i) - max(L(:,i));
    g = exp(tmp)./sum(exp(tmp));
    for k=1:K
        F0m = F0m + g(k).*(L(k,i)-log(K)-log(g(k)+eps));
        if options.families
            F0f = F0f + g(k).*(L(k,i)-log(g(k)+eps)+log(f0(k)));
        end
    end
end


text:
function EO = gaborconvolve_BW(im, nscale, norient, minWaveLength, sigmaOnf, mult)
%---------------------------------------------------
% This code gives the logGabor decomposition
% of input image 'im', 'nacale' is the number
% of scale, 'norient' is the number of
% orientations and 'EO' is the logGabor
% decomposed output. EO is a cell of
% dimension nscale times norient

% This code has used codes from author: Peter Kovesi
% Department of Computer Science & Software Engineering
% The University of Western Australia
% pk@cs.uwa.edu.au  www.cs.uwa.edu.au/~pk
% May 2001
%---------------------------------------------------

[rows cols] = size(im);

% Author: Peter Kovesi
% Department of Computer Science & Software Engineering
% The University of Western Australia
% pk@cs.uwa.edu.au  www.cs.uwa.edu.au/~pk
%
% May 2001

dThetaOnSigma   = 1.5;

imagefft = fft2(im);                 % Fourier transform of image
EO = cell(nscale, norient);          % Pre-allocate cell array

% Pre-compute some stuff to speed up filter construction

x = ones(rows,1) * (-cols/2 : (cols/2 - 1))/(cols/2);
y = (-rows/2 : (rows/2 - 1))' * ones(1,cols)/(rows/2);
radius = sqrt(x.^2 + y.^2);       % Matrix values contain *normalised* radius from centre.
radius(round(rows/2+1),round(cols/2+1)) = 1; % Get rid of the 0 radius value in the middle
                                             % so that taking the log of the radius will
                                             % not cause trouble.

% Precompute sine and cosine of the polar angle of all pixels about the
% centre point

theta = atan2(-y,x);              % Matrix values contain polar angle.
                                  % (note -ve y is used to give +ve
                                  % anti-clockwise angles)
sintheta = sin(theta);
costheta = cos(theta);
clear x; clear y; clear theta;      % save a little memory

thetaSigma = pi/norient/dThetaOnSigma;  % Calculate the standard deviation of the
                                        % angular Gaussian function used to
                                        % construct filters in the freq. plane.
% The low pass filter
lp = fftshift(lowpassfilter([rows,cols],.45,10));   % Radius .4, 'sharpness' 10

% The main loop...

for o = 1:norient,                   % For each orientation.
%   fprintf('Processing orientation %d \n', o);
  angl = (o-1)*pi/norient;           % Calculate filter angle.
  wavelength = minWaveLength;        % Initialize filter wavelength.

  % Pre-compute filter data specific to this orientation
  % For each point in the filter matrix calculate the angular distance from the
  % specified filter orientation.  To overcome the angular wrap-around problem
  % sine difference and cosine difference values are first computed and then
  % the atan2 function is used to determine angular distance.

  ds = sintheta * cos(angl) - costheta * sin(angl);     % Difference in sine.
  dc = costheta * cos(angl) + sintheta * sin(angl);     % Difference in cosine.
  dtheta = abs(atan2(ds,dc));                           % Absolute angular distance.
  spread = exp((-dtheta.^2) / (2 * thetaSigma^2));      % Calculate the angular filter component.

  for s = 1:nscale,                  % For each scale.

    % Construct the filter - first calculate the radial filter component.
    fo = 1.0/wavelength;                  % Centre frequency of filter.
    rfo = fo/0.5;                         % Normalised radius from centre of frequency plane
                                          % corresponding to fo.
    logGabor = exp((-(log(radius/rfo)).^2) / (2 * log(sigmaOnf)^2));
    logGabor(round(rows/2+1),round(cols/2+1)) = 0; % Set the value at the center of the filter
                                                   % back to zero (undo the radius fudge).

    filter = fftshift(logGabor .* lp .* spread); % Multiply by the angular spread to get the filter
                                           % and swap quadrants to move zero frequency
                                           % to the corners.

    % Do the convolution, back transform, and save the result in EO
    EO{s,o} = ifft2(imagefft .* filter);

    wavelength = wavelength * mult;       % Finally calculate Wavelength of next filter
  end                                     % ... and process the next scale

end  % For each orientation

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function Scan=Dicom2FourDAnalyze_MedCom_UCSF(Scan,baseDir);
% Scan=Dicom2FourDAnalyze_MedCom_UCSF(Scan);
% where Scan has to be a struct with the fields
% Scan.DirName
% Scan.Slices
% Scan.Volumes
% Scan.Filenames
% Scan.Filename4d
%
% the script will generate a 4danalyze file
% and put the name into the field Scan.Filename4d
% written 2006.01.12 by Mark Schira mark@ski.org

if ~exist('baseDir')
    baseDir = pwd;
end

Directory=[baseDir,filesep,Scan.DirName];
Filename=Scan.Filenames;
Slices=Scan.Slices;
Volumes=Scan.Volumes;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Stuff for fsl
fslBase='/raid/MRI/toolbox/FSL/fsl';
if (ispref('VISTA','fslBase'))
    disp('Setting fslBase to the one specified in the VISTA matlab preferences:');
    fslBase=getpref('VISTA','fslBase');
    disp(fslBase);
end

fslPath=fullfile(fslBase,'bin'); % This is where FSL lives - should also be able to get this from a Matlab pref
reconPath='/raid/MRI/toolbox/Recon'; % required for the recon program to convert .mag files into Analyze format

cd (Directory);

%******************************************************************
%now do the Recon

try
    In3dFileString=' ';
    for thisVol=1:Volumes;
        fistFile=1+(thisVol-1)*Slices;
        filenumbers=[fistFile:fistFile+Slices-1]; %these are the filenumbers of the slices we need
        filestring=' '; %we need a filestring to start with
        for thisSlice=1:Slices
            filestring=sprintf('%s %s%d.DCM',filestring,Filename,filenumbers(thisSlice));
        end
        Out3dVolStr=sprintf('%s%03d.hdr',Filename,thisVol);
        evalStr1=['!medcon -f',filestring,' -c anlz -stacks -o ',Out3dVolStr];
        eval(evalStr1);
        In3dFileString=sprintf('%s m000-stacks-%s%04d.DCM',In3dFileString,Out3dVolStr);
        In3dFileStringb=sprintf('%s m000-stacks-%s%04d.DCM',In3dFileString,Out3dVolStr);

    end
    OutVolStr=sprintf('%s_4d.hdr',Filename);
    %Usage: avwmerge <-x/y/z/t> <output> <file1 file2 .......>
    evalStr2=['!/raid/MRI/toolbox/FSL/fsl/bin/avwmerge -t ',OutVolStr,In3dFileString];
    eval(evalStr2);
    Scan.Filename4d=OutVolStr;
    display(['saved ',OutVolStr]);

    evalStr=(['! rm m000-stacks*']);
    eval(evalStr);
catch
    display('Warning!! Making4d file failed!');
end
cd (baseDir)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [coords] = GetSquare(x1,x2,x3,x4)

% This function draws a square defined by vertices

% Sanity check
assert(numel(x1) <= 3,"Coordinates can have upto three dimensions.");
assert(iscolumn(x1) == iscolumn(x2) == iscolumn(x3) == iscolumn(x3),"The coordinate vectors must be the same length.");

coords = zeros(5,3);
for i = 1:numel(x1)
    coords(:,i) = [x1(i);x2(i);x3(i);x4(i);x1(i)];
end

end


text:
% POP_TIMEF - Returns estimates and plots of event-related (log) spectral
%           perturbation (ERSP) and inter-trial coherence (ITC) changes
%           timelocked to a set of input events in one data channel.
%
% Usage:
%   >> pop_timef(EEG, typeplot); % pop_up window
%   >> pop_timef(EEG, typeplot, lastcom); % pop_up window
%   >> pop_timef(EEG, typeplot, channel); % do not pop-up
%   >> pop_timef(EEG, typeproc, num, tlimits,cycles,
%                        'key1',value1,'key2',value2, ... );
%
% Inputs:
%   INEEG    - input EEG dataset
%   typeproc - type of processing. 1 process the raw
%              data and 0 the ICA components
%   num      - component or channel number
%   tlimits  - [mintime maxtime] (ms) sub-epoch time limits
%   cycles   -  >0 -> Number of cycles in each analysis wavelet
%               0 -> Use FFTs (with constant window length)
%
% Optional inputs:
%    See the TIMEF function.
%
% Outputs: same as TIMEF, no outputs are returned when a
%          window pops-up to ask for additional arguments
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: TIMEF, EEGLAB

% Copyright (C) 2002 arno@salk.edu, Arnaud Delorme, CNL / Salk Institute
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad
% 03-08-02 add eeglab option & optimize variable sizes -ad
% 03-10-02 change timef call -ad
% 03-18-02 added title -ad & sm
% 04-04-02 added outputs -ad & sm

function varargout = pop_timef(EEG, typeproc, num, tlimits, cycles, varargin );

varargout{1} = '';
% display help if not enough arguments
% ------------------------------------
if nargin < 2
	help pop_timef;
	return;
end;
lastcom = [];
if nargin < 3
	popup = 1;
else
	popup = ischar(num) | isempty(num);
	if ischar(num)
		lastcom = num;
	end
end

% pop up window
% -------------
if popup
	[txt vars] = gethelpvar('timef.m');

	geometry = { [1 0.5 0.5] [1 0.5 0.5] [1 0.5 0.5] [0.92 0.1 0.78] [1 0.5 0.5] [1 0.8 0.2] [1] [1 1]};
    uilist = { ...
                           { 'Style', 'text', 'string', fastif(typeproc, 'Channel number', 'Component number'), 'fontweight', 'bold'  } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,3,[],'1') } {} ...
			   { 'Style', 'text', 'string', 'Epoch time range [min max] (msec)', 'fontweight', 'bold', ...
				 'tooltipstring', 'Sub epoch time limits' } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,4,[],[ int2str(EEG.xmin*1000) ' ' int2str(EEG.xmax*1000) ]) } {} ...
			   { 'Style', 'text', 'string', 'Wavelet cycles (0->FFT, see >> help timef)', 'fontweight', 'bold', ...
				 'tooltipstring', context('cycles',vars,txt) } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,5,[],'3 0.5') } {} ...
			   { 'Style', 'text', 'string',  '[set]->Linear coher / [unset]->Phase coher', 'fontweight', 'bold', ...
				 'tooltipstring', ['Compute linear inter-trial coherence (coher)' 10 ...
					'OR Amplitude-normalized inter-trial phase coherence (phasecoher)'] } ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'phasecoher','present',1) } { } ...
			   { 'Style', 'text', 'string', 'Bootstrap significance level (Ex: 0.01 -> 1%)', 'fontweight', 'bold', ...
				 'tooltipstring', context('alpha',vars,txt) } ...
			   { 'Style', 'edit', 'string', getkeyval(lastcom,'alpha') } {} ...
			   { 'Style', 'text', 'string', 'Optional timef() arguments (see Help)', 'fontweight', 'bold', ...
				 'tooltipstring', 'See timef() help via the Help button on the right...' } ...
			   { 'Style', 'edit', 'string', '''padratio'', 4, ''plotphase'',''off''' } ...
			   { 'Style', 'pushbutton', 'string', 'Help', 'callback',  'pophelp(''timef'');' } ...
			   {} ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotersp','present',0), 'string', ...
				 'Plot Event Related Spectral Power', 'tooltipstring', ...
				 'Plot log spectral perturbation image in the upper panel' } ...
			   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotitc','present',0), 'string', ...
				 'Plot Inter Trial Coherence', 'tooltipstring', ...
				 'Plot the inter-trial coherence image in the lower panel' } ...
			 };
      % { 'Style', 'edit', 'string', '''padratio'', 4, ''plotphase'', ''off''' } ...
			   %{ 'Style', 'text', 'string',  '[set] -> Plot ITC phase sign', 'fontweight', 'bold', ...
			%	 'tooltipstring', ['Plot the sign (+/-) of inter-trial coherence phase' 10 ...
			%		'as red (+) or blue (-)'] } ...
			%   { 'Style', 'checkbox', 'value', ~getkeyval(lastcom,'plotphase','present',1) } { } ...

	result = inputgui( geometry, uilist, 'pophelp(''pop_timef'');', ...
					   fastif(typeproc, 'Plot channel time frequency -- pop_timef()', ...
							  'Plot component time frequency -- pop_timef()'));
	if length( result ) == 0 return; end

	num	     = eval( [ '[' result{1} ']' ] );
	tlimits	 = eval( [ '[' result{2} ']' ] );
	cycles	 = eval( [ '[' result{3} ']' ] );
    if result{4}
    	options = [ ',''type'', ''coher''' ];
    else
		options = [',''type'', ''phasecoher''' ];
    end;

    % add topoplot
    % ------------
    if isfield(EEG.chanlocs, 'theta')
        if ~isfield(EEG, 'chaninfo'), EEG.chaninfo = []; end
        if typeproc == 1
            options = [options ', ''topovec'', ' int2str(num) ...
                        ', ''elocs'', EEG.chanlocs, ''chaninfo'', EEG.chaninfo' ];
        else
            options = [options ', ''topovec'', EEG.icawinv(:,' int2str(num) ...
                       '), ''elocs'', EEG.chanlocs, ''chaninfo'', EEG.chaninfo' ];
        end
    end

    % add title
    % ---------
	if isempty( findstr(  '''title''', result{6}))
        if ~isempty(EEG.chanlocs) && typeproc
            chanlabel = EEG.chanlocs(num).labels;
        else
            chanlabel = int2str(num);
        end
	    switch lower(result{4})
	       case 'coher', options = [options ', ''title'',' fastif(typeproc, '''Channel ', '''Component ') chanlabel ...
	       ' power and inter-trial coherence' fastif(~ isempty(EEG.setname), [' (' EEG.setname ')''' ], '''') ];
	       otherwise, options = [options ', ''title'',' fastif(typeproc, '''Channel ', '''Component ') chanlabel ...
	       ' power and inter-trial phase coherence' fastif(~ isempty(EEG.setname), [' (' EEG.setname ')''' ], '''') ];
	    end
	end
	if ~isempty( result{5} )
		options      = [ options ', ''alpha'',' result{5} ];
	end
	if ~isempty( result{6} )
		  options = [ options ',' result{6} ];
	end
	if ~result{7}
		options = [ options ', ''plotersp'', ''off''' ];
	end
	if ~result{8}
		options = [ options ', ''plotitc'', ''off''' ];
	end
	figure; try, icadefs; set(gcf, 'color', BACKCOLOR); catch, end
else
    options = [ ',' vararg2str(varargin) ];
end

% compute epoch limits
% --------------------
if isempty(tlimits)
	tlimits = [EEG.xmin, EEG.xmax]*1000;
end;
pointrange1 = round(max((tlimits(1)/1000-EEG.xmin)*EEG.srate, 1));
pointrange2 = round(min((tlimits(2)/1000-EEG.xmin)*EEG.srate, EEG.pnts));
pointrange = [pointrange1:pointrange2];

% call function sample either on raw data or ICA data
% ---------------------------------------------------
if typeproc == 1
	tmpsig = EEG.data(num,pointrange,:);
else
	if ~isempty( EEG.icasphere )
        tmpsig = eeg_getdatact(EEG, 'component', num, 'samples', pointrange);
	else
		error('You must run ICA first');
	end;
end;

% outputs
% -------
outstr = '';
if ~popup
    for io = 1:nargout, outstr = [outstr 'varargout{' int2str(io) '},' ]; end
    if ~isempty(outstr), outstr = [ '[' outstr(1:end-1) '] =' ]; end
end

% plot the datas and generate output command
% --------------------------------------------
if length( options ) < 2
    options = '';
end
if nargin < 4
    varargout{1} = sprintf('figure; pop_timef( EEG, %d, %d, %s, %s %s);', typeproc, num, ...
                           vararg2str({tlimits}), vararg2str({cycles}), options);
end
com = sprintf('%s timef( tmpsig(:, :), length(pointrange), [tlimits(1) tlimits(2)], EEG.srate, cycles %s);', outstr, options);
eval(com)

return;

% get contextual help
% -------------------
function txt = context(var, allvars, alltext);
	loc = strmatch( var, allvars);
	if ~isempty(loc)
		txt= alltext{loc(1)};
	else
		disp([ 'warning: variable ''' var ''' not found']);
		txt = '';
	end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
classdef nme_gen_acp < mp.nme_gen_ac & mp.form_acp

%   MATPOWER
%   Copyright (c) 2019-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end     %% properties

%     methods
%     end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function c = medalConstants(query);
%  c = medalConstants(query);%--------------------------------------------------------------------------
% Helper function for MEDAL.
%--------------------------------------------------------------------------
% DES

rootDir = which('medalConstants.m');
tmp = fileparts(rootDir);
rootDir = fileparts(tmp);

cTmp = load(fullfile(rootDir,'.config.mat'),query);
f = fields(cTmp);
c = cTmp.(f{:});

author: dustinstansbury

repo: medal


text:
function rez = learnAndSolve8b(rez, iorder)
% This is the main optimization. Takes the longest time and uses the GPU heavily.

Nbatches = rez.ops.Nbatch;

rng(iseed);
% if getOr(rez.ops, 'midpoint', 0)
%     rez.iorig = randperm(Nbatches);
% end
% rez.istart = ceil(Nbatches/2); % this doesn't really matter anymore

rez.iorig = iorder;

rez     = learnTemplates(rez, rez.iorig);

% if ~isfield(rez, 'W') || isempty(rez.W)
%     rez     = learnTemplates(rez, rez.iorig);
% else
%     rez     = learnTemplates2(rez, rez.iorig);
% end

rez.ops.fig = 0;
rez         = runTemplates(rez);

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
% MVNORMRND - Multivariate Normal - Random Number Generation
% Copyright (c) 1998, Harvard University. Full copyright in the file Copyright
%
% Y = mvnormrnd(mu, sigma, n)
%
%   mu = p by 1 mean column vector or n by p matrix of means
%   sigma = covariance matrix
%   n = number of observations to generate
%
%   Y = an n by p matrix of row vectors with mean mu and covariance sigma
%
% Note: works slightly different from Matlab builtin MVNRND.
%
%   if mu is a column vector, n rows will be returned, all with mean mu
%
%   if mu is a matrix, a matrix of the same size will be returned with
%   row Y(i,:) having mean mu(i,:) .
%
% See also: MVNORMPDF, MVNORMLPR

function [Y] = mvnormrnd (mu,sigma,n)

[d1,d2] = size(mu);
S = chol(sigma)';

if d2==1,
  % then mu is a column vector
  X = normrnd(0,1,n,d1);
  Y = X*S' + ones(n,1)*mu' ;
else
  X = normrnd(0,1,d1,d2);
  Y = X*S' + mu ;
end

author: Sable

repo: mcbench-benchmarks


text:
function objHandle = drawSegmentsObject(O,color,lineWidth)

% DRAWSEGMENTSOBJECT  Draw segments object.
%   DRAWSEGMENTSOBJECT(O)  plots all the segments in object O.
%
%   DRAWSEGMENTSOBJECT(O,C,L) uses segments of color C and width L.
%
%   H = DRAWSEGMENTSOBJECT(...) returns a vector H of line handles to each
%   segment in O.
%
%   See also HOUSE, DRAWSEGMENTSIMAGE.

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

X = O([1 4],:);
Y = O([2 5],:);
Z = O([3 6],:);

if nargin < 3
     lineWidth = 1;
     if nargout < 2
         color = 'b';
     end
end

if nargout == 1

    objHandle = line(X,Y,Z,'color',color,'lineStyle','-','linewidth',lineWidth);

else

    line(X,Y,Z,'color',color,'lineStyle','-','linewidth',lineWidth);

end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


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

%   Copyright (C) 2011 Jonathan Currie (I2C2)

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

%Names and Defaults
Names = {'algorithm','var_lin','cons_lin','allowable_fraction_gap','allowable_gap','cutoff','cutoff_decr','enable_dynamic_nlp','iteration_limit',...
         'nlp_failure_behavior','node_comparison','num_cut_passes','num_cut_passes_at_root',...
         'number_before_trust','number_strong_branch','solution_limit','sos_constraints','tree_search_strategy',...
         'variable_selection','feasibility_pump_objective_norm','heuristic_RINS','heuristic_dive_MIP_fractional',...
         'heuristic_dive_MIP_vectorLength','heuristic_dive_fractional','heuristic_dive_vectorLength',...
         'heuristic_feasibility_pump','pump_for_minlp','candidate_sort_criterion','maxmin_crit_have_sol',...
         'maxmin_crit_no_sol','min_number_strong_branch','number_before_trust_list','number_look_ahead',...
         'number_strong_branch_root','setup_pseudo_frac','coeff_var_threshold','dynamic_def_cutoff_decr','first_perc_for_cutoff_decr',...
         'max_consecutive_infeasible','num_resolve_at_infeasibles',...
         'num_resolve_at_node','num_resolve_at_root','second_perc_for_cutoff_decr','milp_solver','milp_strategy','bb_log_interval',...
         'bb_log_level','lp_log_level','milp_log_level','nlp_log_level','nlp_log_at_root','oa_log_frequency','oa_log_level',...
         'max_consecutive_failures','max_random_point_radius','num_retry_unsolved_random_point','random_point_perturbation_interval',...
         'random_point_type','warm_start'}';
Defaults = {'B-BB',[],[],0,0,1e100,1e-5,'no',2^31-1,'stop','best-bound',1,20,8,20,2^31-1,'enable','probed-dive','osi-strong',1,'no',...
            'yes','no','no','no','no','no','best-ps-cost',0.1,0.7,0,0,0,2^31-1,0.5,0.1,'no',-0.02,0,0,0,0,-0.05,'Cbc_D','find_good_sol',...
            100,1,0,0,0,0,100,1,10,100000,0,1,'Jon','none'}';
%Get IPOPT Defaults if required
if(nargin == 2 && strcmpi(varargin{2},'noipopt'))
    varargin = varargin(1);
else
    iset = ipoptset('bonmin');
    Names = [Names;fieldnames(iset)];
    Defaults = [Defaults;struct2cell(iset)];
end

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

function checkfield(field,value)
%Check a field contains correct data type
switch lower(field)
    %Scalar double
    case {'first_perc_for_cutoff_decr','second_perc_for_cutoff_decr'}
        err = opticheckval.checkScalarDbl(value,field);
    %Scalar non negative double
    case {'min_number_strong_branch','number_look_ahead','number_strong_branch_root','coeff_var_threshold'}
        err = opticheckval.checkScalarNonNeg(value,field);
    %Scalar non zero double
    case {'oa_log_frequency','max_random_point_radius','random_point_perturbation_interval'}
        err = opticheckval.checkScalarGrtZ(value,field);
    %Scalar non negative integer
    case {'iteration_limit','num_cut_passes','num_cut_passes_at_root','number_before_trust','number_strong_branch',...
          'solution_limit','max_consecutive_infeasible','num_resolve_at_infeasibles','num_resolve_at_node',...
          'num_resolve_at_root','bb_log_interval','max_consecutive_failures','num_retry_unsolved_random_point'}
        err = opticheckval.checkScalarIntNonNeg(value,field);
    %Scalar double with bounds
    case 'cutoff_decr'
        err = opticheckval.checkScalarBoundLELE(value,field,-1e10,1e10);
    case {'allowable_fraction_gap','allowable_gap'}
        err = opticheckval.checkScalarBoundLELE(value,field,-1e20,1e20);
    case 'cutoff'
        err = opticheckval.checkScalarBoundLELE(value,field,-1e100,1e100);
    case {'maxmin_crit_have_sol','maxmin_crit_no_sol','setup_pseudo_frac'}
        err = opticheckval.checkScalarBoundLELE(value,field,0,1);
    case 'number_before_trust_list'
        err = opticheckval.checkScalarBoundLEL(value,field,-1,inf);
    %Scalar integer with bounds
    case 'feasibility_pump_objective_norm'
        err = opticheckval.checkScalarIntBoundLELE(value,field,1,2);
    case 'bb_log_level'
        err = opticheckval.checkScalarIntBoundLELE(value,field,0,5);
    case {'lp_log_level','milp_log_level'}
        err = opticheckval.checkScalarIntBoundLELE(value,field,0,4);
    case {'nlp_log_level','oa_log_level'}
        err = opticheckval.checkScalarIntBoundLELE(value,field,0,2);
    case 'nlp_log_at_root'
        err = opticheckval.checkScalarIntBoundLELE(value,field,0,12);
    %Vector integer 0/1
    case {'var_lin','cons_lin'}
        err = opticheckval.checkVector01(value,field);
    %Yes / No
    case {'heuristic_rins','heuristic_dive_mip_fractional','heuristic_dive_mip_vectorlength','heuristic_dive_fractional',...
          'heuristic_dive_vectorlength','heuristic_feasibility_pump','pump_for_minlp','dynamic_def_cutoff_decr'}
        err = opticheckval.checkYesNo(value,field);
    %Yes / No with Yes Crash
    case 'enable_dynamic_nlp'
        err = opticheckval.checkYesCrashNo(value,field);
    %Enable / Disable
    case 'sos_constraints'
        err = opticheckval.checkEnDis(value,field);
    %Misc String methods
    case 'algorithm'
        err = opticheckval.checkValidString(value,field,{'B-BB','B-OA','B-QG','B-Hyb','B-Ecp','B-iFP'});
    case 'nlp_failure_behavior'
        err = opticheckval.checkValidString(value,field,{'stop','fathom'});
    case 'node_comparison'
        err = opticheckval.checkValidString(value,field,{'best-bound','depth-first','breadth-first','dynamic','best-guess'});
    case 'tree_search_strategy'
        err = opticheckval.checkValidString(value,field,{'top-node','dive','probed-dive','dfs-dive','dfs-dive-dynamic'});
    case 'variable_selection'
        err = opticheckval.checkValidString(value,field,{'most-fractional','strong-branching','reliability-branching','qp-strong-branching',...
                                                         'lp-strong-branching','nlp-strong-branching','osi-simple','osi-strong','random'});
    case 'candidate_sort_criterion'
        err = opticheckval.checkValidString(value,field,{'best-ps-cost','worst-ps-cost','most-fractional','least-fractional'});
    case 'milp_solver'
        err = opticheckval.checkValidString(value,field,{'cbc_d','cplex'});
    case 'milp_strategy'
        err = opticheckval.checkValidString(value,field,{'find_good_sol','solve_to_optimality'});
    case 'random_point_type'
        err = opticheckval.checkValidString(value,field,{'jon','andreas','claudia'});
    case 'warm_start'
        err = opticheckval.checkValidString(value,field,{'none','fake_basis','optimum','interior_point'});

    %IPOPT or Unknown Parameter
    otherwise
        err = [];
        try
            %try for ipoptset arg
            iset = ipoptset(field,value);
        catch ME
            throw(ME);
        end
        if(isempty(iset))
            err = MException('OPTI:SetFieldError','Unrecognized parameter name ''%s''.', field);
        end
end
if(~isempty(err)), throw(err); end

function printfields()
%Print out fields with defaults
fprintf(' ----- BONMIN (MIP SOLVER) SETTINGS ----- \n');
fprintf('  BONMIN GENERAL SETTINGS:\n');
fprintf('                          algorithm: [ MIP Algorithm {''B-BB''}, ''B-OA'', ''B-QG'', ''B-Hyb'', ''B-Ecp'', ''B-iFP'' ] \n');
fprintf('                            var_lin: [ Decision Variable Linearity (0 - Nonlinear, 1 - Linear) {zeros(n,1)} ] \n');
fprintf('                           cons_lin: [ Constraint Linearity (0 - Nonlinear, 1 - Linear) {zeros(m,1)} ] \n');

fprintf('\n  BONMIN BRANCH & BOUND SETTINGS:\n');
fprintf('             allowable_fraction_gap: [ Value of Relative Gap when algorithm stops {0} ] \n');
fprintf('                      allowable_gap: [ Value of Absolute Gap when algorithm stops {0} ] \n');
fprintf('                             cutoff: [ The algorithm will only look for values better than cutoff {1e100} ] \n');
fprintf('                        cutoff_decr: [ Amount by which cutoff is decremented below a new best upper bound {1e-5} ] \n');
fprintf('                 enable_dynamic_nlp: [ Enable dynamic linear and quadratic rows addition in NLP ''yes'', {''no''} ] \n');
fprintf('                    iteration_limit: [ Cumulated maximum number of iterations in the algorithm to process continuous relaxations {2^31-1} ] \n');
fprintf('               nlp_failure_behavior: [ Action to take when an NLP is unsolved by IPOPT {''stop''}, ''fathom'' ] \n');
fprintf('                    node_comparison: [ Node selection strategy {''best-bound''}, ''depth-first'', ''breadth-first'', ''dynamic'', ''best-guess'' ] \n');
fprintf('                     num_cut_passes: [ Maximum number of cut passes at regular nodes {1} ] \n');
fprintf('             num_cut_passes_at_root: [ Maximum number of cut passes at root node {20} ] \n');
fprintf('                number_before_trust: [ Number of branches on a variable before its pseudo costs are to be believed in dynamic strong branching {8} ] \n');
fprintf('               number_strong_branch: [ Maximum number of variables considered for strong branching {20} ] \n');
fprintf('                     solution_limit: [ Abort after this many integer feasible solutions have been found (or 0 to deactivate) {2^31-1} ] \n');
fprintf('                    sos_constraints: [ Activate SOS Type 1 constraints {''enable''}, ''disable'' ] \n');
fprintf('               tree_search_strategy: [ Tree traversing strategy {''probed-dive''}, ''top-node'', ''dive'', ''dfs-dive'', ''dfs-dive-dynamic'' ] \n');
fprintf('                 variable_selection: [ Variable selection strategy {''osi-strong''}, ''most-fractional'', ''reliability-branching'', ''qp-strong-branching'', ''lp-strong-branching'', ''nlp-strong-branching'', ''osi-simple'', ''strong-branching'', ''random'' ] \n');

fprintf('\n  BONMIN MINLP HEURISTIC SETTINGS:\n');
fprintf('    feasibility_pump_objective_norm: [ Norm of feasibility pump objective norm {1}, 2 ] \n');
fprintf('                     heuristic_RINS: [ Use RINS heuristic {''no''}, ''yes'' ] \n');
fprintf('      heuristic_dive_MIP_fractional: [ Use Dive MIP fractional heuristic {''no''}, ''yes'' ] \n');
fprintf('    heuristic_dive_MIP_vectorLength: [ Use Dive MIP vectorLength heuristic {''no''}, ''yes'' ] \n');
fprintf('          heuristic_dive_fractional: [ Use Dive fractional heuristic ''no'', {''yes''} ] \n');
fprintf('        heuristic_dive_vectorLength: [ Use Dive vectorLength heuristic ''no'', {''yes''} ] \n');
fprintf('         heuristic_feasibility_pump: [ Use feasibility pump heuristic {''no''}, ''yes'' ] \n');
fprintf('                     pump_for_minlp: [ Use FP for MINLP {''no''}, ''yes'' ] \n');

fprintf('\n  BONMIN NON-CONVEX PROBLEM SETTINGS:\n');
fprintf('                coeff_var_threshold: [ Coefficient of variation threshold (for dynamic definition of cutoff_decr) {0.1} ] \n');
fprintf('            dynamic_def_cutoff_decr: [ Define cutoff_decr dynamically? {''no''}, ''yes'' ] \n');
fprintf('         first_perc_for_cutoff_decr: [ Percentage used when the coeff of variance is smaller than the threshold, to compute cutoff_decr dynamically {-0.02} ] \n');
fprintf('         max_consecutive_infeasible: [ Maximum number of consecutive infeasible subproblems before aborting a branch {0} ] \n');
fprintf('         num_resolve_at_infeasibles: [ Number of tries to resolve an infeasible node (other than the root) of the tree with different starting points {0} ] \n');
fprintf('                num_resolve_at_node: [ Number of tries to resolve a node (other than the root) of the tree with different starting points (like multi-start) {0} ] \n');
fprintf('                num_resolve_at_root: [ Number of tries to resolve the root node with different starting points {0} ] \n');
fprintf('        second_perc_for_cutoff_decr: [ Percentage used when the coeff of variance is greater than the threshold, to compute cutoff_decr dynamically {-0.05} ] \n');

fprintf('\n  BONMIN NLP SOLUTION ROBUSTNESS SETTINGS:\n');
fprintf('           max_consecutive_failures: [ Number of consecutive unsolved problems before aborting a branch of the tree {10} ] \n');
fprintf('            max_random_point_radius: [ Max value for coordinate of a random point {100000} ] \n');
fprintf('    num_retry_unsolved_random_point: [ Number of times that the algorithm will try to resolve an unsolved NLP with a random starting point (with new point) {0} ] \n');
fprintf(' random_point_perturbation_interval: [ Amount by which starting point is perturbed when choosing a random point {1} ] \n');
fprintf('                  random_point_type: [ Random starting point method {''Jon''}, ''Andreas'', ''Claudia'' ] \n');
fprintf('                         warm_start: [ Warm start method {''none''}, ''fake_basis'', ''optimum'', ''interior'' ] \n');

fprintf('\n  BONMIN STRONG BRANCHING SETTINGS:\n');
fprintf('           candidate_sort_criterion: [ Criterion used to choose candidates in strong-branching {''best-ps-cost''}, ''worst-ps-cost'', ''most-fractional'', ''least-fractional'' ] \n')
fprintf('               maxmin_crit_have_sol: [ Weight towards minimum of lower and upper branching estimates when a solution has been found {0.1} ] \n')
fprintf('                 maxmin_crit_no_sol: [ Weight towards minimum of lower and upper branching estimates when no solution has been found {0.7} ] \n')
fprintf('           min_number_strong_branch: [ Minimum number of variables for strong branching (overriding trust) {0} ] \n')
fprintf('           number_before_trust_list: [ Number of branches on a variable vbefore its pseudo costs are to be believed during setup of strong branching list {0} ] \n')
fprintf('                  number_look_ahead: [ Limit of look-ahead strong-branching trials {0} ] \n')
fprintf('          number_strong_branch_root: [ Maximum number of variables considered for strong branching at root node {2^31-1} ] \n')
fprintf('                  setup_pseudo_frac: [ Proportion of strong branching list that has to be taken from most integer infeasible list {0.5} ] \n')

fprintf('\n  BONMIN MILP SOLVER SETTINGS:\n');
fprintf('                        milp_solver: [ Subsolver to solve MILP problems in OA decompositions {''Cbc_D''}, ''Cplex'' (must be installed on your PC) ] \n');
fprintf('                      milp_strategy: [ MILP solving strategy {''find_good_sol''}, ''solve_to_optimality'' ] \n');

fprintf('\n  BONMIN OUPUT AND LOG LEVEL SETTINGS:\n');
fprintf('                    bb_log_interval: [ Interval at which node level output is printed (number of nodes) {100} ] \n');
fprintf('                       bb_log_level: [ Level of branch and bound log detail (0-5) {1} ] \n');
fprintf('                       lp_log_level: [ Level of LP solver log detail (0-4) {0} ] \n');
fprintf('                     milp_log_level: [ Level of MILP solver log detail (0-4) {0} ] \n');
fprintf('                      nlp_log_level: [ Level of NLP solver log detail (0-2) {0} ] \n');
fprintf('                    nlp_log_at_root: [ Level of NLP solver log detail at root node (0-12) {0} ] \n');
fprintf('                   oa_log_frequency: [ Frequency (in seconds) of OA log messages  {100} ] \n');
fprintf('                       oa_log_level: [ Level of OA decomposition log detail (0-2) {1} ] \n');

fprintf('\n ----- IPOPT (RELAXED SOLVER) SETTINGS ----- \n');
ipoptset('bonmin');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function x = add_dc(y,dc,columns)
%ADD_DC Add DC channel to signals.
%   X = ADD_DC(Y,DC) adds the specified DC value to the (possibly
%   multi-dimensional) signal Y, returning the result as X. DC should be a
%   scalar value.
%
%   X = ADD_DC(Y,DC,'columns') where Y is a 2D matrix and DC is an array of
%   length size(Y,2), treats the columns of Y as individual 1D signals,
%   adding to each one the corresponding DC value from the DC array. X is
%   the same size as Y and contains the resulting signals.
%
%   See also REMOVE_DC.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  April 2009

if (nargin==3 && strcmpi(columns,'columns')), columns = 1;
else columns = 0;
end

if (columns)
  x = addtocols(y,dc);
else
  x = y + dc;
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [edata, adjlist, boundmap, perim] = mcmcGetDisjointEdgeData(imsegs, spdata, npertype, ignorelist)
% edata(nadj, nfeatures)
%   edata feature descriptions:
%      01 - 03: abs diff mean rgb
%      04 - 06: abs diff hsv conversion
%      07 - 07: chi-squared hue histogram
%      08 - 08: chi-squared sat histogram
%      09 - 23: abs diff mean texture response
%      24 - 24: chi-squared texture response histogram
%      25 - 26: abs diff mean x-y
%      27 - 27: smaller area / larger area

% spdata(nseg, nf)
%   spdata feature descriptions:
%      01 - 03: mean rgb
%      04 - 06: hsv conversion
%      07 - 11: hue histogram
%      12 - 14: sat histogram
%      15 - 29: mean texture response
%      30 - 44: texture response histogram
%      45 - 46: mean x-y
%      47 - 48: 10th, 90th perc. x
%      49 - 50: 10th, 90th perc. y
%      51 - 51: h / w
%      52 - 52: area

nfeatures = 27;

nseg = imsegs.nseg;

%[boundmap, perim] = mcmcGetSuperpixelBoundaries(imsegs);

nadj = nseg*(nseg-1)/2;

adjlist = zeros(nadj, 2);
edata = zeros(nadj, nfeatures);
c = 0;
for s1 = 1:nseg
    ns1 = nseg-s1;
    adjlist(c+1:c+ns1, 1) = s1;
    adjlist(c+1:c+ns1, 2) = [s1+1:nseg]';
    c = c + ns1;
end

[imh, imw] = size(imsegs.segimage);

for k = 1:nadj
    s1 = adjlist(k, 1);
    s2 = adjlist(k, 2);

    % differences of spdata
    %edata(k, 1:6) = abs(spdata(s1, 1:6) - spdata(s2, 1:6));
    edata(k, 7) = jensenshannon(spdata(s1, 7:11), spdata(s2, 7:11));
    %edata(k, 8) = jensenshannon(spdata(s1, 12:14), spdata(s2, 12:14));
    %edata(k, 9:23) = abs(spdata(s1, 15:29) - spdata(s2, 15:29));
    edata(k, 24) = jensenshannon(spdata(s1, 30:44), spdata(s2, 30:44));
    %edata(k, 25:26) = abs(spdata(s1, 45:46) - spdata(s2, 45:46));
    %edata(k, 27) = min(spdata([s1 s2], 52)) / max(spdata([s1 s2], 52));

end

% ignore adjacent sp (these are handeld in mcmcGetEdgeData)
keep = ones(nadj, 1);
if exist('ignorelist') && ~isempty(ignorelist)
    ignoremat = zeros(nseg);
    for k = 1:size(ignorelist, 1)
        ignoremat(ignorelist(k, 1), ignorelist(k, 2)) = 1;
    end
    for k = 1:nadj
        if ignoremat(adjlist(k, 1), adjlist(k, 2))
            keep(k) = 0;
        end
    end
    edata = edata(find(keep), :);
    adjlist = adjlist(find(keep), :);
end

% keep only top npertype edges for each type (color, texture)
keep = zeros(nadj, 1);
[val, ind] = sort(edata(:, 7), 'ascend'); % hue chi-square stat
keep(ind(1:npertype)) = 1;
[val, ind] = sort(edata(:, 24), 'ascend'); % texture chi-square stat
keep(ind(1:npertype)) = 1;
edata = edata(find(keep), :);
adjlist = adjlist(find(keep), :);

nadj = size(adjlist, 1);

for k = 1:nadj
    s1 = adjlist(k, 1);
    s2 = adjlist(k, 2);

    % differences of spdata
    edata(k, 1:6) = abs(spdata(s1, 1:6) - spdata(s2, 1:6));
    %edata(k, 7) = jensenshannon(spdata(s1, 7:11), spdata(s2, 7:11));
    edata(k, 8) = jensenshannon(spdata(s1, 12:14), spdata(s2, 12:14));
    edata(k, 9:23) = abs(spdata(s1, 15:29) - spdata(s2, 15:29));
    %edata(k, 24) = jensenshannon(spdata(s1, 30:44), spdata(s2, 30:44));
    edata(k, 25:26) = abs(spdata(s1, 45:46) - spdata(s2, 45:46));
    edata(k, 27) = min(spdata([s1 s2], 52)) / max(spdata([s1 s2], 52));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function x = chisquare(hist1, hist2)
x = sum((hist1-hist2).^2 ./ hist2);

function d = jensenshannon(hist1, hist2)
% Jensen-Shannon divergence (see wikipedia)
avehist = (hist1 + hist2)/2;
d = 0.5*(sum(hist1.*log(hist1./hist2)) + sum(hist2.*log(hist2./hist1)));

author: zouchuhang

repo: LayoutNet


text:
function addframe(vw,img)
%WORKED=ADDFRAME(VW)
%  Writes a video frame.  VW must be a videoWriter object and IMG a 2D
%  numeric array with 1 or 3 channels that is autoconverted to a uint8
%  image as follows:
%
%    type    assumed range
%    ----    -------------
%    uint8   0 to 255
%    double  0 to 1
%    logical 0 to 1
%
%SEE ALSO
%  videoWriter
%
%Copyright (c) 2007 Gerald Dalley
%See "MIT.txt" in the installation directory for licensing details (especially
%when using this library on GNU/Linux).

[h,w,d] = size(img);

if (isa(img, 'uint8'))
  if (h ~= vw.h || w ~= vw.w)
    img = uint8(255*imresize(double(img)/255, [vw.h vw.w]));
  else
    % no changes needed
  end

elseif (isa(img, 'double') || islogical(img))
  if (h ~= vw.h || w ~= vw.w)
    img = uint8(255*imresize(img, [vw.h vw.w]));
  else
    img = uint8(255*img);
  end

else
  error('Invalid image type.');
end

if (d == 1)
  img = repmat(img, [1 1 3]);
end

feval(vw.plugin, 'addframe', vw.handle, img);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode