text:
function [engine, loglik, niter] = enter_evidence(engine, evidence, varargin)
% ENTER_EVIDENCE Add the specified evidence to the network (pearl)
% [engine, loglik, num_iter] = enter_evidence(engine, evidence, ...)
% evidence{i} = [] if if X(i) is hidden, and otherwise contains its observed value (scalar or column vector)
%
% The following optional arguments can be specified in the form of name/value pa irs:
% [default value in brackets]
%
% maximize - if 1, does max-product instead of sum-product [0]
% 'filename' -  msgs will be printed to this file, so you can assess convergence while it runs [engine.filename]
%
% e.g., engine = enter_evidence(engine, ev, 'maximize', 1)
%
% For discrete nodes, loglik is the negative Bethe free energy evaluated at the final beliefs.
% For Gaussian nodes, loglik is currently always 0.
%
% 'num_iter' returns the number of iterations used.

maximize = 0;
filename = engine.filename;

% parse optional params
args = varargin;
nargs = length(args);
if nargs > 0
  for i=1:2:nargs
    switch args{i},
     case 'maximize', maximize = args{i+1};
     case 'filename', filename = args{i+1};
     otherwise,
      error(['invalid argument name ' args{i}]);
    end
  end
end

if maximize
  error('can''t handle max-prop yet')
end

engine.maximize = maximize;
engine.filename = filename;
engine.bel = []; % reset if necessary

bnet = bnet_from_engine(engine);
N = length(bnet.dag);
ns = bnet.node_sizes(:);

observed_bitv = ~isemptycell(evidence);
disconnected = find(engine.disconnected_nodes_bitv);
if ~all(observed_bitv(disconnected))
  error(['The following discrete nodes must be observed: ' num2str(disconnected)])
end
msg = init_pearl_msgs(engine.msg_type, engine.msg_dag, ns, evidence);

niter = 1;
switch engine.protocol
 case 'parallel', [msg, niter] = parallel_protocol(engine, evidence, msg);
 case 'tree', msg = tree_protocol(engine, evidence, msg);
 otherwise,
  error(['unrecognized protocol ' engine.protocol])
end
engine.niter = niter;

engine.marginal = cell(1,N);
nodes = find(~engine.disconnected_nodes_bitv);
for n=nodes(:)'
  engine.marginal{n} = compute_bel(engine.msg_type, msg{n}.pi, msg{n}.lambda);
end

engine.evidence = evidence; % needed by marginal_nodes and marginal_family
engine.msg = msg;  % needed by marginal_family

if (nargout >= 2)
  if (engine.msg_type == 'd')
    loglik = bethe_free_energy(engine, evidence);
  else
    loglik = 0;
  end
end

%%%%%%%%%%%

function msg =  init_pearl_msgs(msg_type, dag, ns, evidence)
% INIT_MSGS Initialize the lambda/pi message and state vectors
% msg =  init_msgs(dag, ns, evidence)
%

N = length(dag);
msg = cell(1,N);
observed = ~isemptycell(evidence);
lam_msg = 1;

for n=1:N
  ps = parents(dag, n);
  msg{n}.pi_from_parent = cell(1, length(ps));
  for i=1:length(ps)
    p = ps(i);
    msg{n}.pi_from_parent{i} = mk_msg(msg_type, ns(p));
  end

  cs = children(dag, n);
  msg{n}.lambda_from_child = cell(1, length(cs));
  for i=1:length(cs)
    c = cs(i);
    msg{n}.lambda_from_child{i} = mk_msg(msg_type, ns(n), lam_msg);
  end

  msg{n}.lambda = mk_msg(msg_type, ns(n), lam_msg);
  msg{n}.pi = mk_msg(msg_type, ns(n));

  if observed(n)
    msg{n}.lambda_from_self = mk_msg_with_evidence(msg_type, ns(n), evidence{n});
  else
    msg{n}.lambda_from_self = mk_msg(msg_type, ns(n), lam_msg);
  end
end

%%%%%%%%%

function msg =  mk_msg(msg_type, sz, is_lambda_msg)

if nargin < 3, is_lambda_msg = 0; end

switch msg_type
 case 'd', msg = ones(sz, 1);
 case 'g',
  if is_lambda_msg
    msg.precision = zeros(sz, sz);
    msg.info_state = zeros(sz, 1);
  else
    msg.Sigma = zeros(sz, sz);
    msg.mu = zeros(sz,1);
  end
end

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

function msg = mk_msg_with_evidence(msg_type, sz, val)

switch msg_type
 case 'd',
  msg = zeros(sz, 1);
  msg(val) = 1;
 case 'g',
  %msg.observed_val = val(:);
  msg.precision = inf;
  msg.mu = val(:);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function l = TSTL_append(l, argument)

%tstoolbox/@list/append
%   Syntax:
%     * list = TSTL_append(list, string)
%     * list = TSTL_append(list, list)
%
%   Add string(s) to existing list.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

l.data = l.data(:); 		% ensure compatibility with previous implementation

if ischar(argument)
    l.data = [l.data; {argument}];
	l.len = l.len + 1;
elseif isa(argument, 'list')
	l.len = l.len + argument.len;
	l.data = [l.data; argument.data];
elseif iscellstr(argument)
	l = TSTL_append(l, list(argument));
else
	error('Wrong type of argument(s) given');
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [W,D,E] = cubic_winding_number(C,V,ispoly)
  % CUBIC_WINDING_NUMBER Cubic the winding number of points in V with respect to
  % a cubic Bezier curve in C
  %
  % [W,D,E] = cubic_winding_number(C,V)
  %
  % Inputs:
  %   C  4 by 2 list of control points
  %   V  #V by 2 list of query points
  % Outputs:
  %   W  #V list of winding number values
  %   D  #V list of max-depths of recursive algorithm
  %   E  #V list of total evaluations of recursive algorithm
  %

  function I = inpolygon_convex(V,C);
    %I = inpolygon(V(:,1),V(:,2),C(:,1),C(:,2));
    %I = abs(winding_number(C,[1:size(C,1);2:size(C,1) 1]',V))>0.5;
    %% Geez, this is barely faster than winding_number
    %Nx = C([2:end 1],2)-C(:,2);
    %Ny = C(:,1)-C([2:end 1],1);
    %S = (V(:,1)-C(:,1)').*Nx' + (V(:,2)-C(:,2)').*Ny';
    %I = all(S<0,2);
    % Oh, matlab, why:
    I = all(((V(:,1)-C(:,1)').*(C([2:end 1],2)-C(:,2))' + (V(:,2)-C(:,2)').*(C(:,1)-C([2:end 1],1))')<0,2);
  end

  if nargin<3 || isempty(ispoly)
    ispoly = false;
  end

  max_depth = 40;
  W = zeros(size(V,1),1);
  D = zeros(size(V,1),1);
  E = zeros(size(V,1),1);
  k = 0;
  % "queue"
  Q = {{C,(1:size(V,1))',1}};
  while ~isempty(Q)
    % Pop off back: faster but uglier traversal for animation
    Qi = Q{end};
    Ci = Qi{1};
    Ji = Qi{2};
    depth = Qi{3};
    Q = Q(1:end-1);
    %% Pop off front: better looking traversal for animation, but slower
    %Qi = Q{1};
    %Ci = Qi{1};
    %Ji = Qi{2};
    %depth = Qi{3};
    %Q = Q(2:end);
    D(Ji) = max(D(Ji),depth);
    E(Ji) = E(Ji)+1;
    if depth >= max_depth
      %warning('spline_winding_number exceeded max depth (%d)\n',max_depth);
      continue;
    end
    if size(Ci,1)<3
      I = false(numel(Ji),1);
    else
      H = convhull(Ci);
      I = inpolygon_convex(V(Ji,:),Ci(H(1:end-1),:));
    end
    Wi = winding_number(Ci,[size(Ci,1) 1],V(Ji(~I),:));
    W(Ji(~I)) = W(Ji(~I)) - Wi;
    %WW = zeros(size(W));
    %WW(Ji(~I)) = -Wi;
    %ss = [915 938];
    %surf( ...
    %  reshape(V(:,1),ss), ...
    %  reshape(V(:,2),ss), ...
    %  reshape(0*V(:,2),ss), ...
    %  'CData', ...
    %  reshape(WW,ss), fphong, 'EdgeColor', 'none');
    %hold on
    %[pe,p] = plot_cubic(C);
    %set(p,'Color',0.75*[1 1 1]);
    %arrayfun(@(p) set(p,'Color',1+0.5*(get(p,'Color')-1)),pe);
    %plot_cubic(Ci);
    %plot_edges(Ci,[1 4],':','LineWidth',2,'Color',0.4*[1 1 1]);
    %hold off;
    %view(2);
    %axis equal;
    %axis tight;
    %caxis([-1 1])
    %colormap((flipud(cbrewer('RdBu',45))))
    %set(gca,'YDir','reverse');
    %set(gcf,'color','w');
    %set(gca,'Visible','off','Position',[0 0 1 1]);
    %figpng(sprintf('cubic-winding_number-%02d.png',k));
    %k=k+1;

    Ki = Ji(I);
    if ~isempty(Ki)
      if ispoly
        s = ceil(size(Ci,1)/2);
        C1 = Ci(1:s,:);
        C2 = Ci(s:end,:);
        Q{end+1} = {C1,Ki,Qi{3}+1};
        Q{end+1} = {C2,Ki,Qi{3}+1};
      else
        [C1,C2] = cubic_split(Ci,0.5);
        Q{end+1} = {C1,Ki,Qi{3}+1};
        Q{end+1} = {C2,Ki,Qi{3}+1};
      end
    end
  end
end

author: alecjacobson

repo: gptoolbox


text:
function [varargout]=amfStruct2xml(varargin)

% function [domNode]=amfStruct2xml(parseStruct,fileName,optionStruct)
% ------------------------------------------------------------------------
%
% This function provides the basis for coding AMF XML input file content
% and helps to generate the XML input file.
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2022/09/21 Created based on febioStruct2xml
%------------------------------------------------------------------------

%% Parse input

defaultOptionStruct.attributeKeyword='ATTR';
defaultOptionStruct.valueKeyword='VAL';
defaultOptionStruct.arrayLoopKeywords={''};

defaultOptionStruct.fieldOrder={'metadata','object','material'};

defaultOptionStruct.arrayParseMethod=1;

switch nargin
    case 1
        parseStruct=varargin{1};
        fileName=[];
        optionStruct=[];
    case 2
        parseStruct=varargin{1};
        fileName=varargin{2};
        optionStruct=[];
    case 3
        parseStruct=varargin{1};
        fileName=varargin{2};
        optionStruct=varargin{3};
    otherwise
        error('Wrong number of intput arguments, between 1 and 3 input arguments supported i.e.: parseStruct,fileName,optionStruct')
end

[optionStruct]=structComplete(optionStruct,defaultOptionStruct,1); %Complement provided with default if missing or empty

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

%Create temp folder if it does not exist
if ~exist(savePath,'file')
    mkdir(savePath);
end

%Create output file name if not provided
if isempty(fileName)
    fileName=fullfile(savePath,'temp.xml');
end

%% Initialize XML object

domNode = com.mathworks.xml.XMLUtils.createDocument('amf'); %Create the overall amf field

%% Add comment
% commentString = ['Created using GIBBON, ',char(datetime('now','TimeZone','local','Format','yyyy-dd-MM HH:mm:ss'))];
% commentNode = domNode.createComment(commentString);
% rootNode=domNode.getElementsByTagName('amf').item(0);
% rootNode.appendChild(commentNode);

rootNode=domNode.getElementsByTagName('amf').item(0);
%% Reorder fields
[parseStruct]=reorderStruct(parseStruct,optionStruct.fieldOrder);

%% Add content to XML
[domNode]=struct2xmlStep(domNode,rootNode,parseStruct,optionStruct,fileName);

%% Export XML file
xmlwrite_xerces(fileName,domNode); %Custom XML write function

%Import to cut out bad element names from text write mode
T=txtfile2cell(fileName);

logicKeep=~gcontains(T,'xxxxxxxxxxx_CUT');
cell2txtfile(fileName,T(logicKeep),0,0);

%% Collect output
if nargout>0
    varargout{1}=domNode;
end

end

function [domNode,rootNode]=struct2xmlStep(domNode,rootNode,parseStruct,optionStruct,fileName)

arrayParseMethod=optionStruct.arrayParseMethod;

%%
%Get keywords
attributeKeyword=optionStruct.attributeKeyword;
valueKeyword=optionStruct.valueKeyword;
arrayLoopKeywords=optionStruct.arrayLoopKeywords;

%Check if parseStruct if empty
if isempty(parseStruct) %throw error
    error('Empty structure component encountered. Note this may occur when index is skipped in amf structure')
end

%Get field names occuring in the input structure
fieldNameSet = fieldnames(parseStruct);
attributeStruct=[]; %Initialize as empty
for q_field=1:1:numel(fieldNameSet) %Loop for all field names

    %Get current field names and values
    currentFieldName=fieldNameSet{q_field}; %attribute name or element name
    currentFieldValue=parseStruct.(fieldNameSet{q_field}); %attribute value or element value

    %Switch behaviour depending on name, i.e. create attribute, assign
    %value (to current entity), or create an element with a value
    switch currentFieldName
        case attributeKeyword %ATTRIBUTE
            attributeStruct=currentFieldValue; %The attribute value should be a structure e.g. ATTR.attributeName.attributeValue
            if ~isempty(rootNode)
                %Add the current set of attribute for the element at rootNode
                [domNode]=addAttributeSetXML(domNode,rootNode,attributeStruct);
            else
                error('The rootNode is empty, the current attribute appears to lack a parent element');
            end
        case valueKeyword %VALUE
            %Add the current value to rootNode
            [domNode]=addElementValueXML(domNode,rootNode,currentFieldValue);
        otherwise %ELEMENT

            %Check if the current element tag correspods to one requiring
            %looping
            logicArray=any(strcmp(currentFieldName,arrayLoopKeywords));

            if logicArray %Check for looping across values
                parseStructSub=currentFieldValue; %The current element should be a structure
                if isfield(parseStructSub,valueKeyword) %if the structure contains a fielname which mathes the value keyword
                    arrayDataSet=parseStructSub.(valueKeyword); %Get the array of data
                    parseStructSub = rmfield(parseStructSub,valueKeyword); %Remove the value field from the array
                    numArray=size(arrayDataSet,1); %Measure the number of entries in the array
                else %Some other field
                    arrayDataSet=[];
                    numArray=[];
                end

                if isfield(parseStructSub,attributeKeyword) %if the structure contains a fielname which mathes the attribute keyword
                    attributeStruct=parseStructSub.(attributeKeyword);
                    parseStructSub = rmfield(parseStructSub,attributeKeyword); %Remove the value field from the array
                    fieldNameSetSub=fieldnames(attributeStruct);
                else
                    fieldNameSetSub=[];
                end

                if ~isempty(arrayDataSet) %if array entries exist
                    switch arrayParseMethod
                        case 0
                            for q_array=1:1:numArray %Loop over elements in array
                                %Add the element and its value
                                if iscell(arrayDataSet)
                                    [domNode,elementNode,rootNode]=addElementXML(domNode,rootNode,currentFieldName,arrayDataSet{q_array});
                                else
                                    [domNode,elementNode,rootNode]=addElementXML(domNode,rootNode,currentFieldName,arrayDataSet(q_array,:));
                                end
                                if ~isempty(fieldNameSetSub) %If attributes exist
                                    for q_fieldSub=1:1:numel(fieldNameSetSub) %loop over all attributes
                                        currentFieldNameSub=fieldNameSetSub{q_fieldSub}; %Get current attribute name
                                        %Create structure for attribute creation
                                        if iscell(attributeStruct.(currentFieldNameSub))
                                            attributeStructSub.(currentFieldNameSub)=attributeStruct.(currentFieldNameSub){q_array};
                                        else
                                            attributeStructSub.(currentFieldNameSub)=attributeStruct.(currentFieldNameSub)(q_array,:);
                                        end
                                        [domNode]=addAttributeSetXML(domNode,elementNode,attributeStructSub); %Add attribute
                                    end
                                end
                            end
                        otherwise
                            [domNode,rootNode]=textModeElementAttributeParse(domNode,rootNode,currentFieldName,attributeStruct,arrayDataSet,numArray,fieldNameSetSub,fileName,arrayParseMethod);
                    end

                elseif ~isempty(fieldNameSetSub) %If this is not empty there are attributes to add despite empty array

                    %Get length from first
                    currentFieldNameSub=fieldNameSetSub{1};
                    numArray=numel(attributeStruct.(currentFieldNameSub));

                    switch arrayParseMethod
                        case 0
                            for q_array=1:1:numArray %Loop over all entries
                                [domNode,elementNode,rootNode]=addElementXML(domNode,rootNode,currentFieldName,[]); %Create element without value
                                for q_fieldSub=1:1:numel(fieldNameSetSub) %Loop over and add all attributes
                                    currentFieldNameSub=fieldNameSetSub{q_fieldSub};
                                    %Create structure for attribute creation
                                    if iscell(attributeStruct.(currentFieldNameSub))
                                        attributeStructSub.(currentFieldNameSub)=attributeStruct.(currentFieldNameSub){q_array};
                                    else
                                        attributeStructSub.(currentFieldNameSub)=attributeStruct.(currentFieldNameSub)(q_array,:);
                                    end
                                    [domNode]=addAttributeSetXML(domNode,elementNode,attributeStructSub); %Add attribute
                                end
                            end
                        otherwise
                            [domNode,rootNode]=textModeElementAttributeParse(domNode,rootNode,currentFieldName,attributeStruct,arrayDataSet,numArray,fieldNameSetSub,fileName,arrayParseMethod);
                    end
                elseif ~isempty(parseStructSub) %Despite arrayLoopKeywords trigger no values or attributes are to be added, but other fields might exist
                    [domNode,rootNode]=recursiveElementParse(domNode,rootNode,currentFieldName,currentFieldValue,optionStruct,fileName);
                end

            else
                [domNode,rootNode]=recursiveElementParse(domNode,rootNode,currentFieldName,currentFieldValue,optionStruct,fileName);
            end
    end
end

end

function [domNode,rootNode]=recursiveElementParse(domNode,rootNode,currentFieldName,currentFieldValue,optionStruct,fileName)

if isnumeric(currentFieldValue) || ischar(currentFieldValue)
    [domNode,~,rootNode]=addElementXML(domNode,rootNode,currentFieldName,currentFieldValue);
elseif isstruct(currentFieldValue)
    [domNode,elementNode,rootNode]=addElementXML(domNode,rootNode,currentFieldName,[]); %Create the current element
    [domNode,elementNode]=struct2xmlStep(domNode,elementNode,currentFieldValue,optionStruct,fileName); %Add whatever is nested in here
elseif iscell(currentFieldValue)
    for q_cell=1:1:numel(currentFieldValue)
        parseStructSub=currentFieldValue{q_cell};
        [domNode,elementNode,rootNode]=addElementXML(domNode,rootNode,currentFieldName,[]);
        [domNode,elementNode]=struct2xmlStep(domNode,elementNode,parseStructSub,optionStruct,fileName); %Add whatever is nested in here
    end
end

end

%%

function [domNode,rootNode]=textModeElementAttributeParse(domNode,rootNode,currentFieldName,attributeStruct,arrayDataSet,numArray,fieldNameSetSub,fileName,arrayParseMethod)

switch arrayParseMethod
    case 1 %Non-flexible, assumes equal attribute type and size (fastest)
        textFormat=['<',currentFieldName,' '];
        if ~isempty(fieldNameSetSub) %If attributes exist
            D=[];
            for q_fieldSub=1:1:numel(fieldNameSetSub) %loop over all attributes
                attributeName=fieldNameSetSub{q_fieldSub}; %Get current attribute name
                if iscell(attributeStruct.(attributeName))
                    attributeNameValue=attributeStruct.(attributeName);
                else
                    attributeNameValue=attributeStruct.(attributeName);
                end
                D=[D attributeNameValue];
                if all(isrounded(attributeNameValue))
                    t_form=repmat(['%d',', '],1,size(attributeNameValue,2));
                else
                    t_form=repmat(['%6.7e',', '],1,size(attributeNameValue,2));
                end
                t_form=t_form(1:end-2); %Take away last space and comma
                textFormat=[textFormat,' ',attributeName,'="',t_form,'"'];
            end
        else
            D=[];
        end

        if all(isrounded(arrayDataSet))
            t_form=repmat(['%d',', '],1,size(arrayDataSet,2));
        else
            t_form=repmat(['%6.7e',', '],1,size(arrayDataSet,2));
        end
        t_form=t_form(1:end-2); %Take away last space and comma
        textFormat=[textFormat,'>',t_form,'</',currentFieldName,'>','\n'];
        D=[D arrayDataSet];

        %Convert to string
        t=sprintf(textFormat,D');

        T_add={t};
    case 2 %Flexible, handles non-uniform attributes, uses cellfun to "loop"
        T_add=cell(numArray,1);
        tquote=repmat({'"'},numArray,1);
        if ~isempty(fieldNameSetSub) %If attributes exist
            for q_fieldSub=1:1:numel(fieldNameSetSub) %loop over all attributes
                attributeName=fieldNameSetSub{q_fieldSub}; %Get current attribute name
                attributeNameValue=attributeStruct.(attributeName);
                T_attributeValue=var2cellstr(attributeNameValue);
                T_add=strcat(T_add,repmat({[' ',attributeName,'=']},numArray,1),tquote,T_attributeValue,tquote);
            end
        end
        T_val=var2cellstr(arrayDataSet);
        T_add=strcat(repmat({['<',currentFieldName,' ']},numArray,1),T_add,repmat({'>'},numArray,1),T_val,repmat({['</',currentFieldName,'>']},numArray,1));
    case 3 %Flexible, handles non-uniform attributes, loop through all entries (slowest)
        T_add=cell(numArray,1);
        for q_array=1:1:numArray %Loop over elements in array
            attributeText=[];
            if ~isempty(fieldNameSetSub) %If attributes exist
                for q_fieldSub=1:1:numel(fieldNameSetSub) %loop over all attributes
                    attributeName=fieldNameSetSub{q_fieldSub}; %Get current attribute name
                    if iscell(attributeStruct.(attributeName))
                        attributeNameValue=attributeStruct.(attributeName){q_array};
                    else
                        attributeNameValue=attributeStruct.(attributeName)(q_array,:);
                    end
                    [attributeNameValue]=vec2strIntDouble(attributeNameValue,'%6.7e');
                    attributeText=[attributeText,' ',attributeName,'="',attributeNameValue,'"'];
                end
            end
            if ~isempty(arrayDataSet)
                [currentValue]=vec2strIntDouble(arrayDataSet(q_array,:),'%6.7e');
            else
                currentValue=[];
            end
            T_add{q_array}=['<',currentFieldName,' ',attributeText,'>',currentValue,'</',currentFieldName,'>'];
        end
end

fieldNameCut=[currentFieldName,'_xxxxxxxxxxx_CUT'];
T_top{1}=['<',fieldNameCut,'>'];
totalTextCell=T_top;
totalTextCell(end+1:end+numel(T_add))=T_add;
totalTextCell(end+1)={['</',fieldNameCut,'>']};

% Export text cell to .feb file
cell2txtfile(fileName,totalTextCell,0,0); %Export to text file

% Reimport XML type
domNode2 = xmlread(fileName); %Reimport docNode

nodeList=domNode2.getElementsByTagName(fieldNameCut);
rootNode2=nodeList.item(0);

firstDocImportedNode = domNode.importNode(rootNode2.getChildNodes,true);
rootNode.appendChild(firstDocImportedNode);

%Delete file
delete(fileName);

end

%%
function [B]=reorderStruct(A,fieldOrder)

%Start with fieldOrder fields
fieldNameSet=fieldnames(A);
for q=1:1:numel(fieldOrder)
    fieldNameNow=fieldOrder{q};
    logicFind=strcmp(fieldNameNow,fieldNameSet);
    if any(logicFind)
        B.(fieldNameNow)=A.(fieldNameNow);
        A=rmfield(A,fieldNameNow);
    end
end

%Add remaining fields
fieldNameSet=fieldnames(A);
for q=1:1:numel(fieldNameSet)
    fieldNameNow=fieldNameSet{q};
    B.(fieldNameNow)=A.(fieldNameNow);
end

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 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 [ n_data, j1, j2, j3, j4, j5, j6, j7, j8, j9, fx ] = nine_j_values ...
  ( n_data )

%*****************************************************************************80
%
%% NINE_J_VALUES returns some values of the Wigner 9J function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real J1, J2, J3, J4, J5, J6, J7, J8, J9,
%    the arguments of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 9;

  fx_vec = [ ...
     0.0004270039294528318, ...
    -0.001228915451058514, ...
    -0.0001944260688400887, ...
     0.003338419923885592, ...
    -0.0007958936865080434, ...
    -0.004338208690251972, ...
     0.05379143536399187, ...
     0.006211299937499411, ...
     0.03042903097250921 ];
  j1_vec = [ ...
    1.0, ...
    1.5, ...
    2.0, ...
    1.0, ...
    1.5, ...
    2.0, ...
    0.5, ...
    1.0, ...
    1.5  ];
  j2_vec = [ ...
    8.0, ...
    8.0, ...
    8.0, ...
    3.0, ...
    3.0, ...
    3.0, ...
    0.5, ...
    0.5, ...
    0.5  ];
  j3_vec = [ ...
    7.0, ...
    7.0, ...
    7.0, ...
    2.0, ...
    2.0, ...
    2.0, ...
    1.0, ...
    1.0, ...
    1.0 ];
  j4_vec = [ ...
    6.5, ...
    6.5, ...
    6.5, ...
    4.0, ...
    4.0, ...
    4.0, ...
    2.0, ...
    2.0, ...
    2.0 ];
  j5_vec = [ ...
    7.5, ...
    7.5, ...
    7.5, ...
    1.5, ...
    1.5, ...
    1.5, ...
    1.0, ...
    1.0, ...
    1.0 ];
  j6_vec = [ ...
    7.5, ...
    7.5, ...
    7.5, ...
    3.0, ...
    3.0, ...
    3.0, ...
    1.5, ...
    1.5, ...
    1.5 ];
  j7_vec = [ ...
    6.0, ...
    6.0, ...
    6.0, ...
    3.5, ...
    3.5, ...
    3.5, ...
    1.5, ...
    1.5, ...
    1.5 ];
  j8_vec = [ ...
    10.0, ...
    10.0, ...
    10.0, ...
     2.0, ...
     2.0, ...
     2.0, ...
     0.5, ...
     0.5, ...
     0.5 ];
  j9_vec = [ ...
    6.0, ...
    6.0, ...
    6.0, ...
    2.0, ...
    2.0, ...
    2.0, ...
    1.5, ...
    1.5, ...
    1.5 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    j1 = 0.0;
    j2 = 0.0;
    j3 = 0.0;
    j4 = 0.0;
    j5 = 0.0;
    j6 = 0.0;
    j7 = 0.0;
    j8 = 0.0;
    j9 = 0.0;
    fx = 0.0;
  else
    j1 = j1_vec(n_data);
    j2 = j2_vec(n_data);
    j3 = j3_vec(n_data);
    j4 = j4_vec(n_data);
    j5 = j5_vec(n_data);
    j6 = j6_vec(n_data);
    j7 = j7_vec(n_data);
    j8 = j8_vec(n_data);
    j9 = j9_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function badChannels = getBadChannelStructure()
  badChannels = struct( ...
    'all', [], ...
    'badChannelsFromNaNs', [], ...
    'badChannelsFromNoData', [], ...
    'badChannelsFromLowSNR', [], ...
    'badChannelsFromHFNoise', [], ...
    'badChannelsFromCorrelation', [], ...
    'badChannelsFromDeviation', [], ...
    'badChannelsFromRansac', [], ...
    'badChannelsFromDropOuts', []);

author: VisLab

repo: EEG-Clean-Tools


text:
function [Gbest,CrowdDis] = UpdateGbest(Gbest,N)
% Update the global best set

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    Gbest    = Gbest(NDSort(Gbest.objs,1)==1);
    CrowdDis = CrowdingDistance(Gbest.objs);
    [~,rank] = sort(CrowdDis,'descend');
    Gbest    = Gbest(rank(1:min(N,length(Gbest))));
    CrowdDis = CrowdDis(rank(1:min(N,length(Gbest))));
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%
%
%  Setup_Nesta
%
%  sets up the path for Nesta
%
%
%
% Written by: Jerome Bobin, Caltech
% Email: bobin@acm.caltech.edu
% Created: May 2009
%
% NESTA NESTA Version 1.1
%   See also NESTA

% This adds a relative path, so it only works
% when you're in the NESTA directory

%addpath RecPF_v1.1/solver/
%addpath Misc/

% This adds an absolute path, so it works
% when the user is in any directory

p = mfilename('fullpath'); % the location of this file
path = fileparts(p);
addpath( fullfile(path,'RecPF_v1.1','solver' ) );
addpath( fullfile(path,'Misc') );

repo: REST


text:
function Obs = observationInnovation(Obs,innType)

% OBSERVATIONINNOVATION  Observation innovation.
%   Obs = OBSERVATIONINNOVATION(Obs,innType) updates the structure Obs.inn with the
%   innovation. The used fields are Obs.meas and Obs.exp.
%
%   For line landmarks, ie. when Obs.ltype = '???Lin', the function uses
%   the innovation space defined in innType.
%
%   See also INNOVATION, OBSERVEKNOWNLMKS.

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

switch Obs.ltype(4:6)

    case 'Pnt'   % Use regular Euclidean innovation z = y - h(x)
        [Obs.inn.z, Obs.inn.Z, Obs.inn.iZ, Obs.inn.MD2, Z_e] = innovation(...
            Obs.meas.y, Obs.meas.R,...
            Obs.exp.e, Obs.exp.E);

    case 'Lin'   % Lines don't admit Euclidean innovation. Select one:
        switch innType % innovation type for segments
            case 'ortDst'
                [Obs.inn.z, Obs.inn.Z, Obs.inn.iZ, Obs.inn.MD2, Z_e] = innovation(...
                    Obs.meas.y, Obs.meas.R, ...
                    Obs.exp.e, Obs.exp.E,   ...
                    @hms2hh);

            case 'rhoThe'
                error('??? Line''s ''%s'' innovation not yet implemented.',innType)

            otherwise
                error('??? Unknown innovation type ''%s''.',innType);
        end

    otherwise
        error('??? Unknown landmark type ''%s''.',Obs.ltype);
end

% Jacobians of innovation - the chain rule
Obs.Jac.Z_r = Z_e * Obs.Jac.E_r;
Obs.Jac.Z_s = Z_e * Obs.Jac.E_s;
Obs.Jac.Z_l = Z_e * Obs.Jac.E_l;

% ========== 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 out_fem_knw(FemMat, CondTensor, OutputFile)
% OUT_FEM_KNW Write Cauchy conductivity file (.knw)
%
% INPUT:
%    - TessMat    : Brainstorm FEM head model mesh (tetrahedral or hexahedral)
%    - CondTensor : [nLayers x 6] or [xx, yy, zz, xy, yz, zx]
%    - OutputFile : Full path to output file

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

% Indices of all elements
index_elem = (1:size(FemMat.Tissue, 1))';

% Anisotropic: Tensor value given for each element
if (length(CondTensor) == length(index_elem))
    values = [index_elem, CondTensor];
% Isotropic: Same value for all the elements of a give tissues, or same value for all the elements
else
    if (numel(CondTensor) == 6)
        CondTensor = repmat(CondTensor(:)', length(FemMat.TissueLabels), 1);
    else
        CondTensor = reshape(CondTensor,[],6);
    end
    % Check dimensions
    if size(CondTensor,1) ~= length(FemMat.TissueLabels)
        error(['Input mismatch: ' num2str(size(CondTensor,1)) ' conductivity tensors for  '  num2str(length(FemMat.TissueLabels)) ' tissues.']);
    end
    values = [index_elem, CondTensor(FemMat.Tissue,:)];
end

% Open file
[fid, message] = fopen(OutputFile, 'wt');
if (fid < 0)
    error(['Could not create file: ' message]);
end
% Write file
% Note: the spaces between values are important
fprintf(fid, 'BOI - TENSORVALUEFILE\n');
fprintf(fid, '========================================================\n');
fprintf(fid, '========================================================\n');
fprintf(fid,'BOI - TENSOR\n');
fprintf(fid,'         %d   %1.5f       %1.5f       %1.5f\n             %1.5f       %1.5f       %1.5f\n', values');
fprintf(fid,'EOI - TENSOR\n');
fprintf(fid,'========================================================\n');
fprintf(fid,'========================================================\n');
fprintf(fid,'EOI - TENSORVALUEFILE\n');
% Close file
fclose(fid);

author: brainstorm-tools


text:
% REMOVE_KS2_DUPLICATE_SPIKES2 Double-counted spikes are hard to avoid with
% Kilosort's template matching algorithm since the overall fit can be
% improved by having multiple templates jointly account for a single variable waveform.
%
% This function takes the kilosort2 output rez and identifies pair of
% spikes that are close together in time and space. The temporal threshold
% is give by the parameter OVERLAP_S which is 0.5ms by default and
% the spatial threshold (applied to the template primary sites) is given by
% CHANNEL_SEPARATION_UM and is 100 by default.
%
% From these spike pairs, it identifies the pair with the larger template as
% being the "main" or "reference" cluster and the duplicate spikes from the
% other cluster are removed.
%
% All spike pairs are considered, not just those from CCG-contaminated
% pairs, as in REMOVE_KS2_DUPLICATE_SPIKES2.
%
% Adrian Bondy, 2020
%
%=INPUT
%
%   rez structure
%
%=OPTIONAL INPUT, NAME-VALUE PAIRS
%
%   overlap_s
%       the time interval, in second, within which a sequence of spikes are
%       vetted for duplicates.
%
%   channel_separation_um
%       When the primay channels of two spikes are within this distance, in
%       microns, then the two spikes are vetted for duplicate.
%
%=EXAMPLE
%
%   >> rez = remove_ks2_duplicate_spikes(rez)
function rez = remove_ks2_duplicate_spikes(rez, varargin)
    input_parser = inputParser;
    addParameter(input_parser, 'overlap_s', 5e-4, @(x) (isnumeric(x))) % the temporal window within which pairs of spikes will be considered duplicates (if they are also within the spatial window)
    addParameter(input_parser, 'channel_separation_um', 100, @(x) (ischar(x))) % the spatial window within which pairs of spikes will be considered duplicates (if they are also within the temporal window)
    parse(input_parser, varargin{:});
    P = input_parser.Results;

    spike_times = uint64(rez.st3(:,1));
    spike_templates = uint32(rez.st3(:,2));

    rez.U=gather(rez.U);
    rez.W = gather(rez.W);
    templates = zeros(rez.ops.Nchan, size(rez.W,1), size(rez.W,2), 'single');
    for iNN = 1:size(templates,3)
       templates(:,:,iNN) = squeeze(rez.U(:,iNN,:)) * squeeze(rez.W(:,iNN,:))';
    end
    templates = permute(templates, [3 2 1]); % now it's nTemplates x nSamples x nChannels
    n_spikes=numel(spike_times);
    %% Make sure that the spike times are sorted
    if ~issorted(spike_times)
        [spike_times, spike_idx] = sort(spike_times);
        spike_templates = spike_templates(spike_idx);
    else
        spike_idx=(1:n_spikes)';
    end
    %% deal with cluster 0
    if any(spike_templates==0)
        error('Currently this function can''t deal with existence of cluster 0. Should be OK since it ought to be run first in the post-processing.');
    end
    %% Determine the channel where each spike had that largest amplitude (i.e., the primary) and determine the template amplitude of each cluster
    whiteningMatrix = rez.Wrot/rez.ops.scaleproc;
    whiteningMatrixInv = whiteningMatrix^-1;

    % here we compute the amplitude of every template...
    % unwhiten all the templates
    tempsUnW = zeros(size(templates));
    for t = 1:size(templates,1)
        tempsUnW(t,:,:) = squeeze(templates(t,:,:))*whiteningMatrixInv;
    end

    % The amplitude on each channel is the positive peak minus the negative
    tempChanAmps = squeeze(max(tempsUnW,[],2))-squeeze(min(tempsUnW,[],2));

    % The template amplitude is the amplitude of its largest channel
    [tempAmpsUnscaled,template_primary] = max(tempChanAmps,[],2);
    %without undoing the whitening
    %template_amplitude = squeeze(max(templates, [], 2) - min(templates, [], 2));
    %[~, template_primary] = max(template_amplitude, [], 2);

    template_primary = cast(template_primary, class(spike_templates));
    spike_primary = template_primary(spike_templates);

    %% Number of samples in the overlap
    n_samples_overlap = round(P.overlap_s * rez.ops.fs);
    n_samples_overlap = cast(n_samples_overlap, class(spike_times));
    %% Distance between each channel
    chan_dist = ((rez.xcoords - rez.xcoords').^2 + (rez.ycoords - rez.ycoords').^2).^0.5;
    n_spikes=numel(spike_times);
    remove_idx = [];
    reference_idx = [];
    % check pairs of spikes in the time-ordered list for being close together in space and time.
    % Check pairs that are separated by N other spikes,
    % starting with N=0. Increasing N until there are no spikes within the temporal overlap window.
    % This means only ever computing a vector operation (i.e. diff(spike_times))
    % rather than a matrix one (i.e. spike_times - spike_times').
    diff_order=0;
    while 1==1
        diff_order=diff_order+1;
        fprintf('Now comparing spikes separated by %g other spikes.\n',diff_order-1);
        isis=spike_times(1+diff_order:end) - spike_times(1:end-diff_order);
        simultaneous = isis<n_samples_overlap;
        if any(isis<0)
            error('ISIs less than zero? Something is wrong because spike times should be sorted.');
        end
        if ~any(simultaneous)
            fprintf('No remaining simultaneous spikes.\n');
            break
        end
        nearby = chan_dist(sub2ind(size(chan_dist),spike_primary(1:end-diff_order),spike_primary(1+diff_order:end)))<P.channel_separation_um;
        first_duplicate = find(simultaneous & nearby); % indexes the first (earliest in time) member of the pair
        n_duplicates = length(first_duplicate);
        if ~isempty(first_duplicate)
            fprintf('On iteration %g, %g duplicate spike pairs were identified.\n',diff_order,n_duplicates);
            amps_to_compare=tempAmpsUnscaled(spike_templates([first_duplicate first_duplicate(:)+diff_order]));
            if length(first_duplicate)==1
                amps_to_compare = amps_to_compare(:)'; % special case requiring a dimension change
            end
            first_is_bigger =  diff(amps_to_compare,[],2)<=0;
            remove_idx = [remove_idx ; spike_idx([first_duplicate(~first_is_bigger);(first_duplicate(first_is_bigger)+diff_order)])];
            reference_idx = [reference_idx ; spike_idx([(first_duplicate(~first_is_bigger)+diff_order);first_duplicate(first_is_bigger)])];
        end
    end
    [remove_idx,idx] = unique(remove_idx);
    reference_idx = reference_idx(idx);
    logical_remove_idx = ismember((1:n_spikes)',remove_idx);
    rez = remove_spikes(rez,logical_remove_idx,'duplicate','reference_time',spike_times(reference_idx),...
        'reference_cluster',spike_templates(reference_idx));
end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function p02_title ( )

%*****************************************************************************80
%
%% P02_TITLE prints a title for problem 02.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Problem 02\n' );
  fprintf ( 1, '  Name:       QuadSum\n' );
  fprintf ( 1, '              Davis, Rabinowitz, page 370, #1.\n' );
  fprintf ( 1, '  Region:     0 <= X(i) <= 1\n' );
  fprintf ( 1, '  Integrand:  F(X) = ( sum ( 2 * X(i) - 1 ) )^4\n' );

  return
end


text:
function dicomrt_createmclinks = dicomrt_createmclinks(study,beam_module,int)
% dicomrt_createmclinks(study,beam_module,int)
%
% Create links from original MC phase space file directory to a user specified directory
% this is to avoid PHSPFILENAME being longer than 80 character supported by BEAM/DOSXYZ
%
% rtplanform contains plan data coded (warmed) by dicomrt_mcwarm.
% beam_module is the MC model to use
% int is a parameter which determines an interactive session (int~=0). During interactive session
%     it is possible to specify the name of directories where files will be linked. (OPTIONAL)
%
% Example:
%
% dicomrt_createmclinks(A,'BEAM_Clinac21CD')
%
% will create:
%
% A.createlinks     when executed links phsp files from original BEAM 'BEAM_Clinac21CD' dir
%                   to user defined dir this may be useful when path to original BEAM dir
%                   is too long to fit in the BEAM/DOSXYZ supported 80 characters
%
%
% See also dicomrt_mcwarm, dicomrt_BEAMexport, dicomrt_DOSXYZexport, dicomrt_createmcscripts
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check number of argument
error(nargchk(2,3,nargin))

% Check input
[study,type]=dicomrt_checkinput(study);
if isequal(type,'RTPLAN')==1
    rtplanformc=dicomrt_mcwarm(study);
else
    error('dicomrt_createmcscripts: invalid input data format');
end

if exist('int')==0
    int=0;
end

% Default directories
LOCAL_EGS=['/user/phexternal/phx1pl/cardiff/egs4'];
PHSPFILEDIR=['/user/phexternal/phx1pl/cardiff'];

% Almost ready to write link file: ask some question if session is interactive
if int~=0
    disp(' ');
    disp('Current default for local egs4 directory is: ');
    disp(LOCAL_EGS);
    egs4option = input('Do you want to change this default path Y/N [N]:','s');
    if egs4option == 'Y' | egs4option == 'y';
        LOCAL_EGS = input('Input the filename complete of full path : ','s');
    end

    % PHSPFILEDIR is the directory where phsp files are suppose to be linked
    % this is to avoid PHSPFILENAME being longer than 80 character supported by BEAM/DOSXYZ
    disp(' ');
    disp('Current directory where phsp files will be linked is: ');
    disp(PHSPFILEDIR);
    linkoption = input('Do you want to change this default path Y/N [N]:','s');
    if linkoption == 'Y' | linkoption == 'y';
        PHSPFILEDIR = input('Input the full path : ','s');
    end
end

ORIGDIR=[LOCAL_EGS,'/',beam_module,'/'];
PHSPFILEDIR = [PHSPFILEDIR,'/'];

% Prepare title and strings for internal file identification
MAINtitle=['#MCTP - BEAM link phsp files utility. Case study name: ',inputname(1)];
comma=[','];
space=[' '];

% Write link file
linkfilename=[inputname(1),'.createlinks'];
linkfilenameid=fopen(linkfilename,'w');

if linkfilenameid<0
    error('dicomrt_createmclinks: Error writing script file. Exit now!');
    return
end

time=fix(clock);
fprintf(linkfilenameid,'#!/bin/sh');
fprintf(linkfilenameid,'\n');
fprintf(linkfilenameid,'#****************************************************************************************');
fprintf(linkfilenameid,'\n');
fprintf(linkfilenameid,'#Session started: ');
fprintf(linkfilenameid,date);fprintf(linkfilenameid,'%c',comma);
fprintf(linkfilenameid,'%2i',time(4));
fprintf(linkfilenameid,'%c',':');
fprintf(linkfilenameid,'%2i',time(5));
fprintf(linkfilenameid,'\n');
fprintf(linkfilenameid,MAINtitle);
fprintf(linkfilenameid,'\n');
fprintf(linkfilenameid,'#****************************************************************************************');
fprintf(linkfilenameid,'\n');
fprintf(linkfilenameid,'echo ----------');
fprintf(linkfilenameid,'\n');
for i=1:size(rtplanformc,1) % loop over # beams
    for j=1:size(rtplanformc{i,3},1)
        filename=[inputname(1),'_b',int2str(i),'s',int2str(j),'.egs4phsp1'];
        fprintf(linkfilenameid,'ln -s ');
        fprintf(linkfilenameid,ORIGDIR);
        fprintf(linkfilenameid,filename);
        fprintf(linkfilenameid,space);
        fprintf(linkfilenameid,PHSPFILEDIR);
        fprintf(linkfilenameid,filename);
        fprintf(linkfilenameid,'\n');
    end
end
fclose(linkfilenameid);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ iters, endpt ] = hooke ( nvars, startpt, rho, eps, itermax, f )

%*****************************************************************************80
%
%% HOOKE seeks a minimizer of a scalar function of several variables.
%
%  Discussion:
%
%    This routine find a point X where the nonlinear objective function
%    F(X) has a local minimum.  X is an N-vector and F(X) is a scalar.
%    The objective function F(X) is not required to be differentiable
%    or even continuous.  The program does not use or require derivatives
%    of the objective function.
%
%    The user supplies three things:
%    1) a subroutine that computes F(X),
%    2) an initial "starting guess" of the minimum point X,
%    3) values for the algorithm convergence parameters.
%
%    The program searches for a local minimum, beginning from the
%    starting guess, using the Direct Search algorithm of Hooke and
%    Jeeves.
%
%    This program is adapted from the Algol pseudocode found in the
%    paper by Kaupe, and includes improvements suggested by Bell and Pike,
%    and by Tomlin and Smith.
%
%    The algorithm works by taking "steps" from one estimate of
%    a minimum, to another (hopefully better) estimate.  Taking
%    big steps gets to the minimum more quickly, at the risk of
%    "stepping right over" an excellent point.  The stepsize is
%    controlled by a user supplied parameter called RHO.  At each
%    iteration, the stepsize is multiplied by RHO  (0 < RHO < 1),
%    so the stepsize is successively reduced.
%
%    Small values of rho correspond to big stepsize changes,
%    which make the algorithm run more quickly.  However, there
%    is a chance (especially with highly nonlinear functions)
%    that these big changes will accidentally overlook a
%    promising search vector, leading to nonconvergence.
%
%    Large values of RHO correspond to small stepsize changes,
%    which force the algorithm to carefully examine nearby points
%    instead of optimistically forging ahead.  This improves the
%    probability of convergence.
%
%    The stepsize is reduced until it is equal to (or smaller
%    than) EPS.  So the number of iterations performed by
%    Hooke-Jeeves is determined by RHO and EPS:
%
%      RHO^(number_of_iterations) = EPS
%
%    In general it is a good idea to set RHO to an aggressively
%    small value like 0.5 (hoping for fast convergence).  Then,
%    if the user suspects that the reported minimum is incorrect
%    (or perhaps not accurate enough), the program can be run
%    again with a larger value of RHO such as 0.85, using the
%    result of the first minimization as the starting guess to
%    begin the second minimization.
%
%    Normal use:
%    (1) Code your function F() in the C language;
%    (2) Install your starting guess;
%    (3) Run the program.
%
%    If there are doubts about the result, the computed minimizer
%    can be used as the starting point for a second minimization attempt.
%
%    To apply this method to data fitting, code your function F() to be
%    the sum of the squares of the errors (differences) between the
%    computed values and the measured values.  Then minimize F()
%    using Hooke-Jeeves.
%
%    For example, you have 20 datapoints (T(i), Y(i)) and you want to
%    find A, B and C so that:
%
%      A*t*t + B*exp(t) + C*tan(t)
%
%    fits the data as closely as possible.  Then the objective function
%    F() to be minimized is just
%
%      F(A,B,C) = sum ( 1 <= i <= 20 )
%        ( y(i) - A*t(i)*t(i) - B*exp(t(i)) - C*tan(t(i)) )^2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2008
%
%  Author:
%
%    ALGOL original by Arthur Kaupe.
%    C version by Mark Johnson.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    M Bell, Malcolm Pike,
%    Remark on Algorithm 178: Direct Search,
%    Communications of the ACM,
%    Volume 9, Number 9, September 1966, page 684.
%
%    Robert Hooke, Terry Jeeves,
%    Direct Search Solution of Numerical and Statistical Problems,
%    Journal of the ACM,
%    Volume 8, Number 2, April 1961, pages 212-229.
%
%    Arthur Kaupe,
%    Algorithm 178:
%    Direct Search,
%    Communications of the ACM,
%    Volume 6, Number 6, June 1963, page 313.
%
%    FK Tomlin, LB Smith,
%    Remark on Algorithm 178: Direct Search,
%    Communications of the ACM,
%    Volume 12, Number 11, November 1969, page 637-638.
%
%  Parameters:
%
%    Input, integer NVARS, the number of spatial dimensions.
%
%    Input, real STARTPT(NVARS), the user-supplied
%    initial estimate for the minimizer.
%
%    Input, real RHO, a user-supplied convergence parameter
%    which should be set to a value between 0.0 and 1.0.  Larger values
%    of RHO give greater probability of convergence on highly nonlinear
%    functions, at a cost of more function evaluations.  Smaller
%    values of RHO reduce the number of evaluations and the program
%    running time, but increases the risk of nonconvergence.
%
%    Input, real EPS, the criterion for halting
%    the search for a minimum.  When the algorithm
%    begins to make less and less progress on each
%    iteration, it checks the halting criterion: if
%    the stepsize is below EPS, terminate the
%    iteration and return the current best estimate
%    of the minimum.  Larger values of EPS (such
%    as 1.0e-4) give quicker running time, but a
%    less accurate estimate of the minimum.  Smaller
%    values of EPS (such as 1.0e-7) give longer
%    running time, but a more accurate estimate of
%    the minimum.
%
%    Input, integer ITERMAX, a limit on the number of iterations.
%
%    Input, function handle F, the name of the function routine,
%    which should have the form:
%      function value = f ( x, n )
%
%    Output, integer ITERS, the number of iterations taken.
%
%    Output, real ENDPT(NVARS), the estimate for the
%    minimizer, as calculated by the program.
%
  verbose = 0;

  for i = 1 : nvars
    newx(i) = startpt(i);
  end

  for i = 1 : nvars
    xbefore(i) = startpt(i);
  end

  for i = 1 : nvars
    if ( startpt(i) == 0.0 )
      delta(i) = rho;
    else
      delta(i) = rho * abs ( startpt(i) );
    end
  end

  funevals = 0;
  steplength = rho;
  iters = 0;
  fbefore = f ( newx, nvars );
  funevals = funevals + 1;
  newf = fbefore;

  while ( iters < itermax & eps < steplength )

    iters = iters + 1;

    if ( verbose )

      fprintf ( 1, '\n' );
      fprintf ( 1, '  FUNEVALS = %d, F(X) = %e\n', funevals, fbefore );
      for i = 1 : nvars
        fprintf ( 1, '  %8d  %e\n', i, xbefore(i) );
      end

    end
%
%  Find best new point, one coordinate at a time.
%
    for i = 1 : nvars
      newx(i) = xbefore(i);
    end

    [ newf, newx, funevals ] = best_nearby ( delta, newx, fbefore, nvars, ...
    f, funevals );
%
%  If we made some improvements, pursue that direction.
%
    keep = 1;

    while ( newf < fbefore & keep == 1 )

      for i = 1 : nvars
%
%  Arrange the sign of DELTA.
%
        if ( newx(i) <= xbefore(i) )
          delta(i) = - abs ( delta(i) );
        else
          delta(i) = abs ( delta(i) );
        end
%
%  Now, move further in this direction.
%
        tmp = xbefore(i);
        xbefore(i) = newx(i);
        newx(i) = newx(i) + newx(i) - tmp;
      end

      fbefore = newf;
      [ newf, newx, funevals ] = best_nearby ( delta, newx, fbefore, nvars, ...
        f, funevals );
%
%  If the further (optimistic) move was bad...
%
      if ( fbefore <= newf )
        break;
      end
%
%  Make sure that the differences between the new and the old points
%  are due to actual displacements; beware of roundoff errors that
%  might cause NEWF < FBEFORE.
%
      keep = 0;

      for i = 1 : nvars
        if ( 0.5 * abs ( delta(i) ) < abs ( newx(i) - xbefore(i) ) )
          keep = 1;
          break
        end
      end

    end

    if ( eps <= steplength & fbefore <= newf )
      steplength = steplength * rho;
      for i = 1 : nvars
        delta(i) = delta(i) * rho;
      end
    end

  end

  for i = 1 : nvars
    endpt(i) = xbefore(i);
  end

  return
end


text:
function trans = scaling(varargin)
%SCALING return 3*3 matrix of a scaling in 2 dimensions
%
%   TRANS = scaling(SX, SY);
%   return the matrix corresponding to scaling by SX and SY in the 2
%   main directions.
%   The returned matrix has the form:
%   [SX  0  0]
%   [0  SY  0]
%   [0   0  1]
%
%   TRANS = scaling(SX);
%   Assume SX and SY are equals.
%
%   See also:
%   transforms2d, transformPoint, createTranslation, createRotation
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 07/04/2004.

%   HISTORY
%   04/01/2007: rename as scaling
%   22/04/2009: copy to createScaling and deprecate

% deprecation warning
warning('geom2d:deprecated', ...
    '''scaling'' is deprecated, use ''createScaling'' instead');

% call current implementation
trans = createScaling(varargin{:});

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%% Ellipsoid parameters
%
% Silson et al. have a very different set of ellipses, far more oriented
% than anyone else.  This was pointed out by a few of us at Stanford, and
% in discussing with Reynolds (the AFNI guy) it turned out that they solved
% the ellipse using
%
%    ax^2 + bxy + cy^2
%
% That's the wrong equation.  The correction equation is
%
%    ax^2 + 2bxy + cy^2
%
% So when the true value is b, they calculate a value of 2b.  This script
% analyses what we expect to find in terms of the ratios of major and minor
% axes if we solve erroneously, as per the AFNI calculation published by
% Silson et al.
%
% The formulae for the ellipsoid calculations are taken from Wolfram in
%
%   http://mathworld.wolfram.com/Ellipse.html
%
% The calculation shows that in many cases the ratio of the lengths of the
% major and minor axes differs quite significantly when one makes the AFNI
% error.
%
% The significance of this calculation is that Silson et al. deny that
% there is a significant difference.  Hence, there is a dispute based on
% mathematics.  Let's check this code and see who is right.
%
% Wandell, September 24, 2018

%% General quadratic notes
%
%  ax^2 + 2bxy + cy^2 + 2dx + 2fy + g=0
%
% For the centered ellipse we have d = f = 0.
% We set g = -1 for this case so
%
%   ax^2 + 2bxy + cy^2 = 1
%
% J = det([ a b; b c]) > 0 and
% D = det([ a b 0; b c 0; 0 0 1])
% I = a + c
% We require D .n.e. 0, J > 0 and D/I < 0
%

%% Set parameters

% These are the three parameters that can be changed
a = 4;
b = 1; originalB = b;
c = 3;

%{
a = 2;  b = 2; c = 3;  % 3.225 and 1.618
%}
g = -1;  % We fix g to -1 because we can always scale a,b and c to make it so

%%  The formulae for the axis lengths from the Wolfram web-page
%
% a' = sqrt((2(af^2+cd^2+gb^2-2bdf-acg))/((b^2-ac)[sqrt((a-c)^2+4b^2)-(a+c)]))
% b' = sqrt((2(af^2+cd^2+gb^2-2bdf-acg))/((b^2-ac)[-sqrt((a-c)^2+4b^2)-(a+c)])).
%
% The formulae below simplify because d=f=0 and g = -1;

b = originalB;
if ellipseValidate(a,b,c,g)

    top    = (2*(g*b^2 - a*c*g));
    bottom = (b^2 - a*c)*(sqrt((a - c)^2 + 4*b^2) - (a + c));
    aLength = sqrt(top / bottom);

    top    = (2*(g*b^2 - a*c*g));
    bottom = (b^2 - a*c)*(-1*sqrt((a - c)^2 + 4*b^2) - (a + c));
    bLength = sqrt(top / bottom);

    mx = max(aLength, bLength);
    mn = min(aLength, bLength);
    fprintf('Ratio 1: %f (b = %.2f)\n',mx/mn,b);
else
    fprintf('Not an ellipse when b = %f\n',b);
end

%% If we incorrectly solve using b2 = 2b
%
% It doesn't really matter whether we use b2 = 2b or b2 = b/2
% The b values are related by a factor of 2 and one is right and the other
% is wrong.  Reynolds says that this has no impact on the ratio of the
% lengths of the two axes.  That is bogus, as running most cases show, such
% as the ones above

b = originalB/2;
if ellipseValidate(a,b,c,g)

    top    = (2*(g*b^2 - a*c*g));
    bottom = (b^2 - a*c)*(sqrt((a - c)^2 + 4*b^2) - (a + c));
    aLength = sqrt(top / bottom);

    top    = (2*(g*b^2 - a*c*g));
    bottom = (b^2 - a*c)*(-1*sqrt((a - c)^2 + 4*b^2) - (a + c));
    bLength = sqrt(top / bottom);

    mx = max(aLength, bLength);
    mn = min(aLength, bLength);
    fprintf('Ratio 2: %f (b = %.2f)\n',mx/mn,b);
else
    fprintf('Not an ellipse when b = %f\n',b);
end

%%  Or if we divide by 2

b = originalB*2;
if ellipseValidate(a,b,c,g)

    top    = (2*(g*b^2 - a*c*g));
    bottom = (b^2 - a*c)*(sqrt((a - c)^2 + 4*b^2) - (a + c));
    aLength = sqrt(top / bottom);

    top    = (2*(g*b^2 - a*c*g));
    bottom = (b^2 - a*c)*(-1*sqrt((a - c)^2 + 4*b^2) - (a + c));
    bLength = sqrt(top / bottom);

    mx = max(aLength, bLength);
    mn = min(aLength, bLength);
    fprintf('Ratio 3: %f (b = %.2f)\n',mx/mn,b);
else
    fprintf('Not an ellipse when b = %f\n',b);
end

%%

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
Network MobileNetV2 {
	Layer CONV1 {
		Type: CONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 32, C: 3, R: 1, S: 1, Y:224, X:224 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

//Bottleneck 1
//BottleneckID_repeat_operatorID

	Layer Bottleneck1_1_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 32, C: 32, R: 1, S: 1, Y:112, X:112 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck1_1_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 32, R: 3, S: 3, Y:112, X:112 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck1_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 16, C: 32, R: 1, S: 1, Y:112, X:112 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck1_1_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 16, R: 1, S: 1, Y:112, X:112 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

// Bottleneck 2
	Layer Bottleneck2_1_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 96, C: 16, R: 1, S: 1, Y:112, X:112 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck2_1_2 {
		Type: DSCONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 1, C: 96, R: 3, S: 3, Y:112, X:112 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck2_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 96, C: 96, R: 1, S: 1, Y:56, X:56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck2_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 96, C: 96, R: 1, S: 1, Y:56, X:56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck2_2_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 96, R: 3, S: 3, Y:56, X:56 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck2_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 24, C: 96, R: 1, S: 1, Y:56, X:56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

//Bottleneck 3

	Layer Bottleneck3_1_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 144, C: 24, R: 1, S: 1, Y:56, X:56 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck3_1_2 {
		Type: DSCONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 1, C: 144, R: 3, S: 3, Y:56, X:56 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck3_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 144, C: 144, R: 1, S: 1, Y:28, X:28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck3_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 144, C: 144, R: 1, S: 1, Y:28, X:28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck3_2_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 144, R: 3, S: 3, Y:28, X:28 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck3_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 144, C: 144, R: 1, S: 1, Y:28, X:28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck3_2_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 144, R: 1, S: 1, Y:28, X:28 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck3_3_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 144, C: 144, R: 1, S: 1, Y:28, X:28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck3_3_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 144, R: 3, S: 3, Y:28, X:28 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck3_3_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 32, C: 144, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

//Bottleneck 4

	Layer Bottleneck4_1_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 192, C: 32, R: 1, S: 1, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_1_2 {
		Type: DSCONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 1, C: 192, R: 3, S: 3, Y: 28, X: 28 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 192, C: 192 R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 192, C: 192, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_2_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 192, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 192, C: 192 R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_2_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 192, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_3_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 192, C: 192, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_3_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 192, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_3_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 192, C: 192 R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_3_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 192, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_4_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 192, C: 192, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_4_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 192, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck4_4_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 64, C: 192 R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

// Bottleneck 5

	Layer Bottleneck5_1_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 384, C: 64, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck5_1_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 384, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck5_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 384, C: 384 R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck5_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 384, C: 384, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck5_2_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 384, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck5_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 96, C: 384 R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

// Bottleneck 6

	Layer Bottleneck6_1_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 576, C: 96, R: 1, S: 1, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck6_1_2 {
		Type: DSCONV
		Stride { X: 2, Y: 2 }
		Dimensions { K: 1, C: 576, R: 3, S: 3, Y: 14, X: 14 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck5_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 576, C: 576 R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck6_2_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 576, C: 576, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck6_2_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 576, R: 3, S: 3, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck6_2_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 576, C: 576 R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck6_2_Residual {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 576 R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck6_3_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 576, C: 576, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck6_3_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 576, R: 3, S: 3, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck6_3_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 160, C: 576 R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

// Bottleneck 7

	Layer Bottleneck7_1_1 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 960, C: 160, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck7_1_2 {
		Type: DSCONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1, C: 960, R: 3, S: 3, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer Bottleneck7_1_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 320, C: 960, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

// CONV2D

	Layer CONV2D_2 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1280, C: 320, R: 1, S: 1, Y: 7, X: 7 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

	Layer CONV2D_3 {
		Type: CONV
		Stride { X: 1, Y: 1 }
		Dimensions { K: 1000, C: 1280, R: 1, S: 1, Y: 1, X: 1 }
		Dataflow {
			TemporalMap (1,1) K;
			TemporalMap (1,1) C;
			TemporalMap (Sz(R),1) Y;
			SpatialMap (Sz(S),1) X;
			TemporalMap (Sz(R),Sz(R)) R;
			TemporalMap (Sz(S),Sz(S)) S;
		}
	}

}

author: maestro-project

repo: maestro


text:
function [ y ] = ill_si( t, x, lambda)
%ILL_SIR Summary of this function goes here
%   输入参数 x 包含一个分量，为：Infective 的比例
% lambda   日接触率

y = [lambda * x(1) * (1 - x(1))]';

end

repo: NumericalAnalysis


text:
% pop_readbdf() - Read Biosemi 24-bit BDF file
%
% Usage:
%   >> EEG = pop_readbdf;             % an interactive window pops up
%   >> EEG = pop_readbdf( filename ); % no pop-up window
%   >> EEG = pop_readbdf( filename, range, eventchans, ref, delchan);
%
% Graphical interface:
%   "Data block range to read" - {edit box] see command line
%                    parameter 'range' below.
%   "Event channel index(s)" - {edit box] see command line
%                    parameter 'eventchans' below.
%   "Index(s) of reference channel(s)" - {edit box] see command line
%                    parameter 'ref' below.
% Inputs:
%   filename       - Biosemi 24-bit BDF file name
%
% Optional input:
%   range          - [min max] integer range of data blocks to import.
%                    Default is empty -> import all data blocks.
%   eventchans     - [integer] event channel index(s). Default: none.
%   ref            - [integer] channel index or index(s) for the reference.
%                    Reference channels are not removed from the data,
%                    allowing easy re-referencing. If more than one
%                    channel, data are referenced to the average of the
%                    indexed channels. WARNING! Biosemi Active II data
%                    are recorded reference-free, but LOSE 40 dB of SNR
%                    if no reference is used!. If you do not know which
%                    channel to use, pick one and then re-reference after
%                    the channel locations are read in. {default: none}
%   delchan        - ['on'|'off'] delete event channel. Default if 'on'\
%
% Outputs:
%   EEG            - EEGLAB data structure
%
% Author: Arnaud Delorme, CNL / Salk Institute, 13 March 2002
%
% See also: openbdf(), readbdf()

% Copyright (C) 13 March 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%adbdf( filename ); % no pop-up window
%   >> EEG = pop_readbdf( filename, range, eventchans, ref, delchan);
%
% Graphical interface:
%   "Data block range to read" - {edit box] see command line
%                    parameter 'range' below.
%   "Event channel index(s)" - {edit box] see command line
%                    parameter 'eventchans' below.
%   "Index(s) of reference channel(s)" - {edit box] see command line
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% programmed from pop_readedf() version 1.15

function [EEG, command] = pop_readbdf(filename, blockrange, eventchans, ref, delchan);
EEG = [];
command = '';

if nargin < 2
    blockrange = [];
end;
if nargin < 3
    eventchans = [];
end;
if nargin < 4
    ref = [];
end;
if nargin < 5
    delchan = 'on';
end;

if nargin < 1
	% ask user
	[filename, filepath] = uigetfile('*.BDF;*.bdf', 'Choose an BDF file -- pop_readbdf()');
    drawnow;

	if filename == 0 return; end;
	filename = [filepath filename];
    disp('We highly recommend that you choose a reference channel IF these are Biosemi data');
    disp('(e.g., a mastoid or other channel). Otherwise the data will lose 40 dB of SNR!');

    dat = openbdf(filename);
    promptstr    = { [ 'Data block range to read (default all [1 ' int2str(dat.Head.NRec) '])' ]
                     [ 'Event channel number(s) (default:none [1-' int2str(dat.Head.NS)  '])' ]
                     'If Biosemi data, reference chan(s) number(s)'
                     'Delete event channel' };
    inistr       = { '' int2str(dat.Head.NS) '' 'on' };
    result       = inputdlg2( promptstr, 'Import BDF file -- pop_readbdf()', 1,  inistr, 'pop_readbdf');
    if length(result) == 0 return; end;
    blockrange   = eval( [ '[' result{1} ']' ] );
    eventchans   = eval( [ '[' result{2} ']' ] );
    ref          = eval( [ '[' result{3} ']' ] );
    delchan      = result{4};
end;

% load datas
% ----------
EEG = eeg_emptyset;
fprintf('Reading BDF data in 24-bit format...\n');
dat = openbdf(filename);
%dat.Head = sopen(filename);
if isempty(dat.Head.NRec)
    dat.Head.NRec = 100;
end;
if isempty(blockrange)
    blockrange = [1 dat.Head.NRec];
end;
vectrange = [blockrange(1):min(blockrange(2), dat.Head.NRec)];
DAT=readbdf(dat,vectrange);
EEG.nbchan          = size(DAT.Record,1);
EEG.srate           = dat.Head.SampleRate(1);
EEG.data            = DAT.Record;
EEG.pnts            = size(EEG.data,2);
EEG.trials          = 1;
EEG.setname 		= 'BDF file';
disp('Event information might be encoded in the last channel');
disp('To extract these events, use menu File > Import event info > From data channel');
EEG.filename        = filename;
EEG.filepath        = '';
EEG.xmin            = 0;
EEG.chanlocs        = struct('labels', cellstr(dat.Head.Label));
EEG = eeg_checkset(EEG);

% extract events
% --------------
disp('Extracting events...');
warning off;
if ~isempty(eventchans)
    if any(eventchans > EEG.nbchan) | any(eventchans < 1)
        error('Event channel index or indices out of range');
    end;
    disp('Applying 8-bit masks to event channels (only for BIOSEMI files).');
    %EEG.data(eventchans,:) = bitand(uint32(EEG.data(eventchans,:)), 32767);

    %thiscode = 0;lastcode=0;
    %for p = 1:size(EEG.data,2)-1
    %    prevcode = thiscode;
    %    thiscode = mod(EEG.data(end,p),256*256);   % andrey's code - 16 bits
    %    if (thiscode ~= 0) && (thiscode~=prevcode) && (thiscode~=lastcode) % fix to avoid repeated codes (per Ying's demand)
    %        EEG.event(end+1).latency =  p;
    %        EEG.event(end).type = thiscode;
    %        lastcode = thiscode;
    %    end;
    %end;
    thiscode = 0;
    for p = 1:size(EEG.data,2)-1
        prevcode = thiscode;
        thiscode = mod(EEG.data(end,p),256*256);   % andrey's code - 16 bits
        if (thiscode ~= 0) && (thiscode~=prevcode)
            EEG.event(end+1).latency =  p;
            EEG.event(end).type = thiscode;
        end;
    end;
    if strcmpi(delchan, 'on')
        EEG = pop_select(EEG, 'nochannel', size(EEG.data,1));
    end;
    EEG = eeg_checkset(EEG, 'makeur');

    %EEG = pop_chanevent(EEG, eventchans, 'edge', 'leading', 'delchan', delchan);
end;
warning on;

% rerefencing
% -----------
if ~isempty(ref)
    disp('Re-referencing...');
    if ~isempty(eventchans) & any(ref > eventchans(1))
        error('Event channel index cannot be before ref channel index, contact arno@salk.edu');
    end;
    EEG.data = EEG.data - repmat(mean(EEG.data(ref,:),1), [size(EEG.data,1) 1]);
    if length(ref) == 1
        disp([ 'Warning: channel ' int2str(ref) ' is now zeroed (but still present in the data)' ]);
    else
        disp([ 'Warning: data matrix rank has decreased through re-referencing' ]);
    end;
end;

% history
% -------
if ~isempty(blockrange) | ~isempty(eventchans)
    command = sprintf('EEG = pop_readbdf(''%s'', %s);', filename, vararg2str({blockrange eventchans ref}));
else
    command = sprintf('EEG = pop_readbdf(''%s'');', filename);
end;

return;

author: PatternRecognition

repo: OpenBMI


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QPSK demonstration packet-based transceiver for Chilipepper
% Toyon Research Corp.
% http://www.toyon.com/chilipepper.php
% Created 10/17/2012
% embedded@toyon.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% There are two major goals with this core. The first is to find the peak
% of the training sequence and then to subsequently pull out and pack the
% bits. The number of bytes transmitted is in the packet so we extract this
% to determine how many bytes to pull out.
% The second goal is to send these bytes off to the Microblaze processor.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%#codegen
function [byte_out, en_out, s_out, o_out] = ...
    qpsk_rx_correlator(s_i_in, s_q_in)

    persistent counter
    persistent sBuf_i sBuf_q
    persistent oLatch sLatch
    persistent q detPacket
    persistent ip op
    persistent bits symCount byteCount numBytes
    t_i = TB_i;
    t_q = TB_q;
    OS_RATE = 8;
    BIT_TO_BYTE = [1 2 4 8 16 32 64 128]';

    if isempty(counter)
        counter = 0;
        sBuf_i = zeros(1,65);
        sBuf_q = zeros(1,65);
        sLatch = 0;
        oLatch = 0;
        q = 0;
        detPacket = 0;
        ip = 0; op = 0;
        bits = zeros(1,8);
        symCount = 0;
        byteCount = 0;
        numBytes = 1000;
    end
    byte_out = 0;
    en_out = 0;

    % found a packet, now we're ready to write the data out
    if counter == 0 && detPacket == 1
        if s_i_in < 0
            sHard_i_t = -1;
        else
            sHard_i_t = 1;
        end
        if s_q_in < 0
            sHard_q_t = -1;
        else
            sHard_q_t = 1;
        end
        sHard_i = 0; sHard_q = 0;
        switch q
            case 0
                sHard_i = sHard_i_t;
                sHard_q = sHard_q_t;
            case 1
                sHard_i = sHard_q_t;
                sHard_q = -sHard_i_t;
            case 2
                sHard_i = -sHard_i_t;
                sHard_q = -sHard_q_t;
            case 3
                sHard_i = -sHard_q_t;
                sHard_q = sHard_i_t;
        end
        sLatch = sHard_i;
        oLatch = 1;
        bits(symCount*2+1) = (sHard_i+1)/2;
        bits(symCount*2+2) = (sHard_q+1)/2;

        symCount = symCount + 1;
        if symCount >= 4
            byteCount = byteCount + 1;
            symCount = 0;
            byte_out = bits*BIT_TO_BYTE;
            en_out = 1;
            % first byte is number of bytes in payload
            if byteCount == 1
                numBytes = byte_out;
            end
            % if we exceed the packet ID
            if byteCount > 3
                % exit if we've written all the bytes or above reasonable
                % threshold
                if byteCount == numBytes+6 || byteCount > 256
                    detPacket = 0;
                    counter = 1;
                end
            end
        end
    end

    % let's see if we can find a packet. only do so if MCU is ok to rcv packet
    if counter == 0 && detPacket == 0
        sLatch = 0;
        if s_i_in < 0
            ss_i = -1;
        else
            ss_i = 1;
        end
        if s_q_in < 0
            ss_q = -1;
        else
            ss_q = 1;
        end

        sBuf_i = [sBuf_i(2:end) ss_i];
        sBuf_q = [sBuf_q(2:end) ss_q];

        sc_iWithi = sBuf_i*t_i;
        sc_iWithq = sBuf_i*t_q;
        sc_qWithi = sBuf_q*t_i;
        sc_qWithq = sBuf_q*t_q;

        ip = abs(sc_iWithi)+abs(sc_qWithq);
        op = abs(sc_iWithq)+abs(sc_qWithi);

        % we found a packet. While we have frequency offset lock we don't
        % know the phase offset. Here we use the inphase and quadrature
        % phasing to determine how to rotate around the circle
        if ip > 100 % 0 or 180 angle
            if sc_iWithi > 10 && sc_qWithq > 10
                q = 0; % 0 degrees
            else
                q = 2; % 180 degrees;
            end
            detPacket = 1;
        end
        if op > 100
            if sc_iWithq > 10 && sc_qWithi < 10
                q = 3; % 90 degrees
            else
                q = 1; % 270 degrees;
            end
            detPacket = 1;
        end
        oLatch = ip+op;
        symCount = 0;
        byteCount = 0;
        numBytes = 1000;
    end

    s_out = sLatch;
    o_out = oLatch;

    % only pull data once every OS_RATE clocks
    counter = counter + 1;
    if counter >= OS_RATE
        counter = 0;
    end
end

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Bayesian Linear Regression 1D Example  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%          1) Generate 1D Regression Datasets                %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Generate some data
clear all;close all;clc;
nbSamples = 50;
X         = linspace(0,4,nbSamples);
w         = [0.5 25];
y         = X * w(1) + normrnd(0,0.2,1,nbSamples) + w(2);

options             = [];
options.points_size = 20;
options.title       = 'Training data';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1      = ml_plot_data([X(:),y(:)],options);

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                 2) Train BLR Model                       %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Train Baysian Gaussian Linear regressin
var         = 0.2;        % variance in the noise of the data
Sigma_p     = eye(2,2);
Xb          = [X];
yb          = y;
[w,invA]    = train_blr(Xb',y',var,Sigma_p);

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  3) Test BLR Model                       %%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Prediction

x_test  = [linspace(0,4,100);ones(1,100)];
y_test  = x_test' * w;
Sigma   = x_test' * invA * x_test;

if exist('h2','var') && isvalid(h2), delete(h2);end
if exist('h3','var') && isvalid(h3), delete(h3);end

h2 = figure;

subplot(1,2,1);
box on;
plot_gaussian_contour(gca,zeros(2,1),Sigma_p);
xlabel('b','FontSize',11);
ylabel('a','FontSize',11);
title('Prior weights','FontSize',11);
grid on;
axis square;
axis_limits = axis;

subplot(1,2,2);
box on;
plot_gaussian_contour(gca,w,invA);
xlabel('b','FontSize',11);
ylabel('a','FontSize',11);
title('Posterior weights','FontSize',11);
grid on;
axis square;

h3 = figure;
hold on;box on;
scatter(X,y,20,'filled');
hp = plot(x_test(1,:),y_test,'--r','LineWidth',2);
grid on;
axis square;
legend(hp,'regression line');
title('y = b * x + a');

author: epfl-lasa


text:
function obj = tide_fac(obj,t_s,t_e,incnstit)
% obj = tide_fac(obj,t_s,t_e,incnstit)
% Input a msh class obj, start and time simulation times, and constituent
% names and return populated f15 tidal potential options in the msh class
% obj
%
% Created by William Pringle March 15 2018,
% Copied functions from U_tide to calculate the tide factors, available at:
% https://www.mathworks.com/matlabcentral/fileexchange/
% 46523--utide--unified-tidal-analysis-and-prediction-functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Get the times and lat entries
ts = datenum(t_s);
te = datenum(t_e);

tref = 0.5*(ts + te);

% make a fake t with 1000 intervals
t = linspace(ts,te,1000)';

lat = mean(obj.p(:,2));

%% Allow for an entry just to be major8 to get the major eight constituents
jj = strcmp(incnstit,'major8');
if sum(jj) > 0
   incnstit = [{'M2','S2','N2','K2','K1','O1','Q1','P1'} incnstit(~jj)];
   disp('Using the major eight harmonic tidal constituents')
   disp(incnstit)
end

%% Now check and get the indices of the incnstits (this part came from U_tide)
load('tide_fac_constants.mat','const');
cnstit.NR.lind = nan*ones(size(incnstit,1),1);
for j = 1:length(incnstit)
    lind1 = strmatch(incnstit{j},const.name);
    if isempty(lind1)
        error(['tide_fac: unrecognized non-reference constituent: '...
               incnstit{j} '.']);
    else
        cnstit.NR.lind(j) = lind1;
    end
end
[~,seq] = sort(const.freq(cnstit.NR.lind));
cnstit.NR.lind = cnstit.NR.lind(seq);
if isempty(cnstit.NR.lind)
    error('tide_fac: no constituents specified');
end
cnstit.NR.frq = const.freq(cnstit.NR.lind);
cnstit.NR.name = cellstr(const.name(cnstit.NR.lind,:));
% Doodson amp multiplied by Doodson constant/g:
% (See : https://www.whoi.edu/fileserver.do?id=21351&pt=10&p=17272)
cnstit.NR.amp  = abs(const.doodsonamp(cnstit.NR.lind))*0.2675;
cnstit.NR.red  = const.earthreduc(cnstit.NR.lind);

%% Now get the F, U, and V (using the U_tide function ut_FUV)
[F,U,V] = ut_FUV(t,tref,cnstit.NR.lind,lat,zeros(1,4));

%% Final factors
F = mean(F)'; % This is the average nodal factor
phs = (mean(U)+V(1,:))'*360;  % This is the average nodal correction +
                              % astronomical argument at beginning of
                              % simulation

% Make sure phase between 0 and 360 for aesthetic purposes
while any(phs < 0)
    phs(phs<0) = phs(phs<0) + 360;
end

obj.f15.ntif = length(cnstit.NR.name);
for k = 1:obj.f15.ntif
    obj.f15.tipotag(k).name   = cnstit.NR.name{k}; % name in frq order
    obj.f15.tipotag(k).val(1) = cnstit.NR.amp(k); % potential amplitude of species
    obj.f15.tipotag(k).val(2) = cnstit.NR.frq(k)*2*pi/3600; % frq in rad/s format
    obj.f15.tipotag(k).val(3) = cnstit.NR.red(k); % earth rigidity reduction factor
    obj.f15.tipotag(k).val(4) = F(k);   % average nodal factor
    obj.f15.tipotag(k).val(5) = phs(k); % average nodal correction +
                         % astronomical argument at beginning of simulation
end
%EOF
end

% Stuff below has been copied and slightly adjusted from U_tide
% (which was copied from T_tide)
function [F,U,V] = ut_FUV(t,tref,lind,lat,ngflgs)
% UT_FUV()
% compute nodal/satellite correction factors and astronomical argument
% inputs
%   t = times [datenum UTC] (nt x 1)
%   tref = reference time [datenum UTC] (1 x 1)
%   lind = list indices of constituents in ut_constants.mat (nc x 1)
%   lat = latitude [deg N] (1 x 1)
%   ngflgs = [NodsatLint NodsatNone GwchLint GwchNone] each 0/1
% output
%   F = real nodsat correction to amplitude [unitless] (nt x nc)
%   U = nodsat correction to phase [cycles] (nt x nc)
%   V = astronomical argument [cycles] (nt x nc)
% UTide v1p0 9/2011 d.codiga@gso.uri.edu
% (uses parts of t_vuf.m from t_tide, Pawlowicz et al 2002)

nt = length(t);
nc = length(lind);
%% nodsat
if ngflgs(2) % none
    F = ones(nt,nc);
    U = zeros(nt,nc);
else
    if ngflgs(1) % linearized times
        tt = tref;
    else         % exact times
        tt = t;
    end
    ntt = length(tt);
    load('tide_fac_constants.mat');
    [astro,~]=ut_astron(tt');
    if abs(lat) < 0.5 % to make sure no close to zero values
        lat = sign(lat)*0.5;
        if sign(lat) == 0
            lat = 0.5;
        end
    end
    slat=sin(pi*lat/180);
    rr=sat.amprat;
    j=find(sat.ilatfac==1);
    rr(j)=rr(j).*0.36309.*(1.0-5.0.*slat.*slat)./slat;
    j=find(sat.ilatfac==2);
    rr(j)=rr(j).*2.59808.*slat;
    uu=rem( sat.deldood*astro(4:6,:)+sat.phcorr(:,ones(1,ntt)), 1);
    nfreq=length(const.isat); %#ok
    mat = rr(:,ones(1,ntt)).*exp(1i*2*pi*uu);
    F = ones(nfreq,ntt);
    ind = unique(sat.iconst);
    for i = 1:length(ind)
        F(ind(i),:) = 1+sum(mat(sat.iconst==ind(i),:),1);
    end
    U = imag(log(F))/(2*pi); % faster than angle(F)
    F=abs(F);
    for k=find(isfinite(const.ishallow))'
        ik=const.ishallow(k)+(0:const.nshallow(k)-1);
        j = shallow.iname(ik);
        exp1 = shallow.coef(ik);
        exp2 = abs(exp1);
        F(k,:)=prod(F(j,:).^exp2(:,ones(ntt,1)),1);
        U(k,:)=sum(U(j,:).*exp1(:,ones(ntt,1)),1);
    end
    F=F(lind,:)';
    U=U(lind,:)';
    if ngflgs(1) % nodal/satellite with linearized times
        F = F(ones(nt,1),:);
        U = U(ones(nt,1),:);
    end
end
%% gwch (astron arg)
if ngflgs(4) % none (raw phase lags not greenwich phase lags)
    if ~exist('const','var')
        load('tide_fac_constants.mat','const');
    end
    [~,ader] = ut_astron(tref);
    ii=isfinite(const.ishallow);
    const.freq(~ii) = (const.doodson(~ii,:)*ader)/(24);
    for k=find(ii)'
        ik=const.ishallow(k)+(0:const.nshallow(k)-1);
        const.freq(k)=sum(const.freq(shallow.iname(ik)).*shallow.coef(ik));
    end
    V = 24*(t-tref)*const.freq(lind)';
else
    if ngflgs(3)  % linearized times
        tt = tref;
    else
        tt = t;   % exact times
    end
    ntt = length(tt);
    if exist('astro','var')
        if ~isequal(size(astro,2),ntt)
            [astro,~]=ut_astron(tt');
        end
    else
        [astro,~]=ut_astron(tt');
    end
    if ~exist('const','var')
        load('tide_fac_constants.mat');
    end
    V=rem( const.doodson*astro+const.semi(:,ones(1,ntt)), 1);
    for k=find(isfinite(const.ishallow))'
        ik=const.ishallow(k)+(0:const.nshallow(k)-1);
        j = shallow.iname(ik);
        exp1 = shallow.coef(ik);
        V(k,:) = sum(V(j,:).*exp1(:,ones(ntt,1)),1);
    end
    V=V(lind,:)';
    if ngflgs(3)    % linearized times
        [~,ader] = ut_astron(tref);
        ii=isfinite(const.ishallow);
        const.freq(~ii) = (const.doodson(~ii,:)*ader)/(24);
        for k=find(ii)'
            ik=const.ishallow(k)+(0:const.nshallow(k)-1);
            const.freq(k)=sum( const.freq(shallow.iname(ik)).* ...
                shallow.coef(ik) );
        end
        V = V(ones(1,nt),:) + 24*(t-tref)*const.freq(lind)';
    end
end
%EOF
end

function [astro,ader] = ut_astron(jd)
% UT_ASTRON()
% calculate astronomical constants
% input
%   jd = time [datenum UTC] (1 x nt)
% outputs
%   astro = matrix [tau s h p np pp]T, units are [cycles] (6 x nt)
%   ader = matrix of derivatives of astro [cycles/day] (6 x nt)
% UTide v1p0 9/2011 d.codiga@gso.uri.edu
% (copy of t_astron.m from t_tide, Pawlowicz et al 2002)

d=jd(:)'-datenum(1899,12,31,12,0,0);
D=d/10000;
args=[ones(size(jd));
      d;
      D.*D;
      D.^3];
sc= [ 270.434164,13.1763965268,-0.0000850, 0.000000039];
hc= [ 279.696678, 0.9856473354, 0.00002267,0.000000000];
pc= [ 334.329556, 0.1114040803,-0.0007739,-0.00000026];
npc=[-259.183275, 0.0529539222,-0.0001557,-0.000000050];
ppc=[ 281.220844, 0.0000470684, 0.0000339, 0.000000070];
astro=rem( [sc;hc;pc;npc;ppc]*args./360.0 ,1);
tau=rem(jd(:)',1)+astro(2,:)-astro(1,:);
astro=[tau;astro];
dargs=[zeros(size(jd));
       ones(size(jd));
       2.0e-4.*D;
       3.0e-4.*D.*D];
ader=[sc;hc;pc;npc;ppc]*dargs./360.0;
dtau=1.0+ader(2,:)-ader(1,:);
ader=[dtau;ader];
end

author: CHLNDDEV