text:
classdef (InferiorClasses = {?rotation,?quaternion}) orientation < rotation
%
% The class *orientation* represents orientations and misorientations.
%
% Syntax
%   ori = orientation(rot)
%   ori = orientation.byEuler(phi1,Phi,phi2,cs,ss)
%   ori = orientation.byEuler(alpha,beta,gamma,'ZYZ',cs,ss)
%   ori = orientation.byMiller([h k l],[u v w],cs,ss)
%   ori = orientation.byAxisAngle(v,omega,cs,ss)
%   ori = orientation.byMatrix(A,cs)
%   ori = orientation.map(u1,v1,u2,v2,cs)
%   ori = orientation.goss(cs)
%   ori = orientation.cube(cs)
%
% Input
%  rot    - @rotation
%  cs, ss - @crystalSymmetry / @specimenSymmetry
%
% Output
%  ori    - @orientation
%
% Class Properties
%  CS, SS          - @crystalSymmetry / @specimenSymmetry
%  antipodal       - grain exchange symmetry for misorientations
%  phi1, Phi, phi2 - Euler angles
%  i               - inversion
%  a, b, c, d      - quaternion components
%
% See also
% DefinitionAsCoordinateTransform CrystalOperations CrystalReferenceSystem

properties

  CS = crystalSymmetry('1');   % crystal symmetry
  SS = specimenSymmetry('1');  % specimen symmetry or crystal symmetry
  antipodal = false

end

methods

  function o = orientation(varargin)

    % find and remove symmetries
    args  = cellfun(@(s) isa(s,'symmetry'),varargin,'uniformoutput',true);
    sym = varargin(args);
    varargin(args) = [];

    % call rotation constructor
    o = o@rotation(varargin{:});

    if nargin == 0, return;end

    % set symmetry
    if ~isempty(varargin) && isa(varargin{1},'orientation')
      o.CS = varargin{1}.CS;
      o.SS = varargin{1}.SS;
      o.antipodal = varargin{1}.antipodal;
    elseif ~isempty(varargin) && ischar(varargin{1}) && strcmpi(varargin{1},'map')
      if isa(varargin{2},'Miller'), o.CS = varargin{2}.CS; end
      if isa(varargin{3},'Miller'), o.SS = varargin{3}.CS; end
    else
      try %#ok<TRYNC>
        a = get_option(varargin,'axis');
        o.CS = a.CS;
        o.SS = a.CS;
      end
    end
    if ~isempty(sym), o.CS = sym{1};end
    if length(sym) > 1, o.SS = sym{2};end

    % empty constructor -> done
    if isempty(varargin), return; end

    % copy constructor
    switch class(varargin{1})

      case 'char'

        switch lower(varargin{1})

          case 'miller'

            o = orientation.byMiller(varargin{:});

          otherwise

            if exist([varargin{1},'Orientation'],'file')

              % there is a file defining this specific orientation
              o = eval([varargin{1},'Orientation(o.CS,o.SS)']);

            end
        end
    end
    o.antipodal = o.antipodal | check_option(varargin,'antipodal');
    if o.antipodal && o.CS ~= o.SS
      warning('antipodal symmetry is only meaningfull for misorientations between the same phase.')
    end
  end
end

methods (Static = true)

  function ori = nan(varargin)
    s = varargin(cellfun(@isnumeric,varargin));
    q = quaternion.nan(s{:});
    ori = orientation(q,varargin{:});
  end

  function ori = id(varargin)
    id = find(~cellfun(@isnumeric,varargin),1)-1;
    q = quaternion.id(varargin{1:id});
    ori = orientation(q,varargin{id+1:end});
  end

  function ori = rand(varargin)
    q = quaternion.rand(varargin{:});
    ori = orientation(q,varargin{:});
  end

  ori = byMiller(m1,m2,varargin);
  ori = byEuler(phi1,Phi,phi2,varargin);
  ori = byAxisAngle(v,omega,varargin);
  ori = byMatrix(M,varargin);
  ori = map(varargin);
  ori = fit(varargin);
  [ori,interface,options] = load(fname,varargin);

  function ori = cube(varargin)
    ori = orientation.byEuler(0,0,0,varargin{:});
  end

  function ori = cubeND22(varargin)
    ori = orientation.byEuler(22*degree,0,0,varargin{:});
  end

  function ori = cubeND45(varargin)
    ori = orientation.byEuler(45*degree,0,0,varargin{:});
  end

  function ori = cubeRD(varargin)
    ori = orientation.byEuler(0,22*degree,0,varargin{:});
  end

  function ori = goss(varargin)
    ori = orientation.byEuler(0,45*degree,0,varargin{:});
  end

  function ori = copper(varargin)
    %ori = orientation.byEuler(90*degree,35*degree,45*degree,varargin{:});
    ori = orientation.byMiller([1 1 2],[-1 1 -1],varargin{:});
  end

  function ori = copper2(varargin)
    %ori = orientation.byEuler(270*degree,35*degree,45*degree,varargin{:});
    ori = orientation.byMiller([1 1 2],[1 -1 1],varargin{:});
  end

  function ori = SR(varargin)
    ori = orientation.byEuler(53*degree,35*degree,63*degree,varargin{:});
  end

  function ori = SR2(varargin)
    ori = orientation.byEuler(233*degree,35*degree,63*degree,varargin{:});
  end

  function ori = SR3(varargin)
    ori = orientation.byEuler(307*degree,35*degree,27*degree,varargin{:});
  end

  function ori = SR4(varargin)
    ori = orientation.byEuler(127*degree,35*degree,27*degree,varargin{:});
  end

  function ori = brass(varargin)
    ori = orientation.byEuler(35*degree,45*degree,0*degree,varargin{:});
  end

  function ori = brass2(varargin)
    ori = orientation.byEuler(325*degree,45*degree,0*degree,varargin{:});
  end

  function ori = PLage(varargin)
    ori = orientation.byEuler(65*degree,45*degree,0*degree,varargin{:});
  end

  function ori = PLage2(varargin)
    ori = orientation.byEuler(245*degree,45*degree,0*degree,varargin{:});
  end

  function ori = QLage(varargin)
    ori = orientation.byEuler(65*degree,20*degree,0*degree,varargin{:});
  end

  function ori = QLage2(varargin)
    ori = orientation.byEuler(245*degree,20*degree,0*degree,varargin{:});
  end

  function ori = QLage3(varargin)
    ori = orientation.byEuler(115*degree,160*degree,0*degree,varargin{:});
  end

  function ori = QLage4(varargin)
    ori = orientation.byEuler(295*degree,160*degree,0*degree,varargin{:});
  end

  function ori = invGoss(varargin)
    ori = orientation.byEuler(90*degree,45*degree,0*degree,varargin{:});
  end

  function mori = Bain(csGamma,csAlpha)
    %
    % Syntax:
    %   mori = Bain(csGamma,csAlpha)
    %
    % Input
    %  csGamma - parent @crystalSymmetry (cubic fcc)
    %  csAlpha - child @crystalSymmetry (cubic bcc)
    %

    mori = orientation.map(Miller(1,0,0,csGamma),Miller(1,0,0,csAlpha),...
      Miller(0,1,0,csGamma,'uvw'),Miller(0,1,1,csAlpha,'uvw'));
  end

  function mori = KurdjumovSachs(csGamma,csAlpha)
    %
    % Syntax:
    %   mori = KurdjumovSachs(csGamma,csAlpha)
    %
    % Input
    %  csGamma - parent @crystalSymmetry (cubic fcc)
    %  csAlpha - child @crystalSymmetry (cubic bcc)
    %

    mori = orientation.map(Miller(1,1,1,csGamma),Miller(0,1,1,csAlpha),...
      Miller(-1,0,1,csGamma,'uvw'),Miller(-1,-1,1,csAlpha,'uvw'));
  end

  function mori = NishiyamaWassermann(csGamma,csAlpha)
    %
    % Syntax:
    %   mori = NishiyamaWassermann(csGamma,csAlpha)
    %
    % Input
    %  csGamma - parent @crystalSymmetry (cubic fcc)
    %  csAlpha - child @crystalSymmetry (cubic bcc)
    %

    mori = orientation.map(Miller(1,1,1,csGamma),Miller(0,1,1,csAlpha),...
      Miller(1,1,-2,csGamma,'uvw'),Miller(0,-1,1,csAlpha,'uvw'));
  end

  function mori = Pitsch(csGamma,csAlpha)
    %
    % Syntax:
    %   mori = Pitch(csGamma,csAlpha)
    %
    % Input
    %  csGamma - parent @crystalSymmetry (cubic fcc)
    %  csAlpha - child @crystalSymmetry (cubic bcc)
    %

    mori = orientation.map(Miller(0,1,0,csGamma),Miller(1,0,1,csAlpha),...
      Miller(1,0,1,csGamma,'uvw'),Miller(-1,1,1,csAlpha,'uvw'));

    %mori = orientation.map(Miller(1,1,0,csGamma),Miller(1,1,1,csAlpha),...
    %  Miller(0,0,1,csGamma,'uvw'),Miller(-1,1,0,csAlpha,'uvw'));

  end

  function mori = GreningerTrojano(csGamma,csAlpha)
    %
    % Syntax:
    %   mori = orientation.GreningerTrojano(csGamma,csAlpha)
    %
    % Input
    %  csGamma - parent @crystalSymmetry (cubic fcc)
    %  csAlpha - child @crystalSymmetry (cubic bcc)
    %

    mori = inv(orientation.byEuler(2.7*degree,46.6*degree,7.5*degree,csAlpha,csGamma));

    %mori = orientation.map(Miller(1,1,1,csGamma),Miller(1,1,0,csAlpha),...
    %  Miller(5,12,17,csGamma,'uvw'),Miller(17,17,7,csAlpha,'uvw'));

  end

  function mori = Burgers(csParent,csChild)
    %
    % Syntax:
    %   mori = orientation.Burgers(csParent,csChild)
    %
    % Input
    %  csParent - parent @crystalSymmetry
    %  csChild  - child @crystalSymmetry
    %

    if csParent.lattice.isTriHex
      mori = inv(orientation.Burgers(csChild,csParent));
    else
      mori = orientation.map(Miller(1,1,0,csParent),Miller(0,0,0,1,csChild),...
        Miller(-1,1,-1,csParent),Miller(2,-1,-1,0,csChild));
    end

  end

  function mori = Burger(varargin)

    warning('orientation.Burger has been renamed to orientation.Burgers')

    mori = orientation.Burgers(varargin{:});

  end

  function mori = ShojiNishiyama(csGamma,csEps)
    %
    % Syntax:
    %   mori = orientation.ShojiNishiyama(csGamma,csEps)
    %
    % Input
    %  csGamma - parent @crystalSymmetry
    %  csEps  - child @crystalSymmetry
    %

    mori = orientation.map(Miller(1, 1, 1,csGamma),Miller(0, 0, 0, 1,csEps),...
      Miller(1, 0, -1,csGamma,'uvw'),Miller(2, -1, -1, 0,csEps,'uvw'));

  end

end

methods (Static = true, Hidden = true)

  function check_dot(cs)

    % first setting
    cs = crystalSymmetry('m-3m');
    n = 10000000;
    ori1 = orientation.rand(n,cs);
    ori2 = orientation.rand(n,cs);

    tic
    mori = inv(ori1) .* ori2;

    d1 = max(dot_outer(mori,cs.rot,'noSymmetry'),[],2);
    toc

    tic
    d2 = dot(ori1,ori2);
    toc

    norm(d1-d2)

    %second setting
    n = 1000000;
    ori1 = orientation.rand(n,cs,cs);
    ori2 = orientation.rand(cs,cs);

    tic
    d1 = dot(ori1,ori2);
    toc

    tic
    d2 = dot(ori2,ori1);
    toc

    norm(d1-d2)

  end

end

end

author: mtex-toolbox

repo: mtex

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


text:
function [mmps2] = uGal2mmps2(uGal)
% Convert acceleration from microgals to millimeters per square second.
% Chad A. Greene 2012
mmps2 = uGal*1e-5;
end

author: Sable

repo: mcbench-benchmarks


text:
function [version] = hasyokogawa(desired)

% HASYOKOGAWA tests whether the data input toolbox for MEG systems by
% Yokogawa (www.yokogawa.com, designed by KIT/EagleTechnology) is
% installed. Only the newest version of the toolbox is accepted.
%
% Use as
%   [string]  = hasyokogawa;
% which returns a string describing the toolbox version, e.g. "12bitBeta3",
% "16bitBeta3", or "16bitBeta6". An empty string is returned if the toolbox
% is not installed. The string "unknown" is returned if it is installed but
% the version is unknown.
%
% Alternatively you can use it as
%   [boolean] = hasyokogawa(desired);
% where desired is a string with the desired version.
%
% See also READ_YOKOGAWA_HEADER, READ_YOKOGAWA_DATA, READ_YOKOGAWA_EVENT,
% YOKOGAWA2GRAD

% Copyright (C) 2010, Tilmann Sander-Thoemmes
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: hasyokogawa.m 2885 2011-02-16 09:41:58Z roboos $

% return empty if not present
version = [];

try
  warning('off', 'MATLAB:pfileOlderThanMfile');
  % Call some functions with input argument "Inf": If
  % the functions are present they return their revision number.
  % Call first GetMeg160ADbitInfoM as this is not present in
  % the 12bit library, in case of error the "catch" part will take over.
  % The code below is intentionally very literal for easy of reading.
  res = textscan(evalc('GetMeg160ADbitInfoM(Inf);'),'%s %s %c %s %s %d');
  rev_ADbitInfoM = res{6};
  res = textscan(evalc('GetMeg160ChannelInfoM(Inf);'),'%s %s %c %s %s %d');
  rev_ChannelInfoM = res{6};
  res = textscan(evalc('GetMeg160AmpGainM(Inf);'),'%s %s %c %s %s %d');
  rev_AmpGainM = res{6};
  res = textscan(evalc('GetMeg160MatchingInfoM(Inf);'),'%s %s %c %s %s %d');
  rev_MatchingInfoM = res{6};
  if [0 2 1 5] == [rev_ADbitInfoM rev_ChannelInfoM rev_AmpGainM rev_MatchingInfoM]
    version='16bitBeta3';
  elseif [0 2 2 5] == [rev_ADbitInfoM rev_ChannelInfoM rev_AmpGainM rev_MatchingInfoM]
    version='16bitBeta6';
  else
    warning('Yokogawa toolbox is installed, but the version cannot be determined.');
    version = 'unknown';
  end
  if nargin>0
    version = strcmpi(version, desired);
    if ~version
        warning('The required version of the Yokogawa input toolbox (%s) is not installed.', desired);
    end
  end
  warning('on', 'MATLAB:pfileOlderThanMfile');
catch
  warning('on', 'MATLAB:pfileOlderThanMfile');
  m = lasterror;
  m.identifier;
  if strcmp(m.identifier, 'MATLAB:UndefinedFunction') || strcmp(m.identifier, 'MATLAB:FileIO:InvalidFid')
    if (exist('GetMeg160ChannelInfoM') && exist('GetMeg160ContinuousRawDataM'));
      version = '12bitBeta3';
    else
      version = 'unknown';
   end
  end
  if nargin>0
    version = 0; % logical output
  end

end

author: PatternRecognition

repo: OpenBMI


text:
%%  The MAP algorithm
%---input---------------------------------------------------------
%   X: initial 2D labels
%   Y: image
%   Z: 2D constraints
%   mu: vector of means
%   sigma: vector of standard deviations
%   k: number of labels
%   MAP_iter: maximum number of iterations of the MAP algorithm
%   show_plot: 1 for showing a plot of energy in each iteration
%       and 0 for not showing
%---output--------------------------------------------------------
%   X: final 2D labels
%   sum_U: final energy

%   Copyright by Quan Wang, 2012/04/25
%   Please cite: Quan Wang. HMRF-EM-image: Implementation of the
%   Hidden Markov Random Field Model and its Expectation-Maximization
%   Algorithm. arXiv:1207.3510 [cs.CV], 2012.

function [X sum_U]=MRF_MAP(X,Y,Z,mu,sigma,k,MAP_iter,show_plot)

[m n]=size(Y);
x=X(:);
y=Y(:);
U=zeros(m*n,k);
sum_U_MAP=zeros(1,MAP_iter);
for it=1:MAP_iter % iterations
    fprintf('  Inner iteration: %d\n',it);
    U1=U;
    U2=U;

    for l=1:k % all labels
        yi=y-mu(l);
        temp1=yi.*yi/sigma(l)^2/2;
        temp1=temp1+log(sigma(l));
        U1(:,l)=U1(:,l)+temp1;

        for ind=1:m*n % all pixels
            [i j]=ind2ij(ind,m);
            u2=0;
            if i-1>=1 && Z(i-1,j)==0
                u2=u2+(l ~= X(i-1,j))/2;
            end
            if i+1<=m && Z(i+1,j)==0
                u2=u2+(l ~= X(i+1,j))/2;
            end
            if j-1>=1 && Z(i,j-1)==0
                u2=u2+(l ~= X(i,j-1))/2;
            end
            if j+1<=n && Z(i,j+1)==0
                u2=u2+(l ~= X(i,j+1))/2;
            end
            U2(ind,l)=u2;
        end
    end
    U=U1+U2;
    [temp x]=min(U,[],2);
    sum_U_MAP(it)=sum(temp(:));

    X=reshape(x,[m n]);
    if it>=3 && std(sum_U_MAP(it-2:it))/sum_U_MAP(it)<0.0001
        break;
    end
end

sum_U=0;
for ind=1:m*n % all pixels
    sum_U=sum_U+U(ind,x(ind));
end
if show_plot==1
    figure;
    plot(1:it,sum_U_MAP(1:it),'r');
    title('sum U MAP');
    xlabel('MAP iteration');
    ylabel('sum U MAP');
    drawnow;
end


text:
function varargout = process_extract_pthresh2( varargin )
% PROCESS_EXTRACT_PTHRESH2 Apply a statistical threshold from FilesA(stat) to a FilesB(regular file).
%
% USAGE:  OutputFiles = process_extract_pthresh2('Run', sProcess, sInput, sInput2)

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

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Apply statistic threshold';
    sProcess.Category    = 'File2';
    sProcess.SubGroup    = 'Test';
    sProcess.Index       = 140;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Statistics#Directionality:_Difference_of_absolute_values';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'pdata', 'presults', 'ptimefreq', 'pmatrix'};
    sProcess.OutputTypes = {'data',  'results',  'timefreq',  'matrix'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 1;
    sProcess.isPaired    = 1;
    % Define options
    sProcess = process_extract_pthresh('DefineOptions', sProcess);
end

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

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInput, sInput2) %#ok<DEFNU>
    % Get options
    [StatThreshOptions, strCorrect] = process_extract_pthresh('GetOptions', sProcess);
    % Load stat result (FileA)
    [sFileA, matNameA] = in_bst(sInput.FileName);
    % Load file to threshold (FileB)
    [sFileB, matNameB] = in_bst(sInput2.FileName);
    % Check that the two files have the same type and dimensions
    sizeMat = size(sFileA.(matNameA));
    if ~isequal(sizeMat, size(sFileB.(matNameB)))
        bst_report('Error', sProcess, [sInput, sInput2], 'File sizes do not match.');
        return;
    end

    % Process separately the types of files
    switch (sInput.FileType)
        case 'pdata'
            % Load channel file
            ChannelMat = in_bst_channel(sInput.ChannelFile);
            % Get only relevant sensors as multiple tests
            iChannels = good_channel(ChannelMat.Channel, sFileA.ChannelFlag, {'MEG', 'EEG', 'SEEG', 'ECOG', 'NIRS'});
            if isfield(sFileA, 'pmap') && ~isempty(sFileA.pmap)
                sFileA.pmap = sFileA.pmap(iChannels,:,:);
            end
            if isfield(sFileA, 'tmap') && ~isempty(sFileA.tmap)
                sFileA.tmap = sFileA.tmap(iChannels,:,:);
            end
            % Create a new data file structure
            pmask = zeros(sizeMat);
            pmask(iChannels,:) = (process_extract_pthresh('Compute', sFileA, StatThreshOptions) ~= 0);

        case {'presults', 'ptimefreq', 'pmatrix'}
            pmask = (process_extract_pthresh('Compute', sFileA, StatThreshOptions) ~= 0);
    end

    % Apply threshold mask to the data file
    sFileB.(matNameB) = sFileB.(matNameB) .* pmask;
    % Update file comment
    sFileB.Comment = [sFileB.Comment ' | ' strCorrect];
    % Add history entry
    sFileB = bst_history('add', sFileB, 'pthresh', ['Setting the stat threshold: ' strCorrect]);
    sFileB = bst_history('add', sFileB, 'pthresh', ['Stat file: ' sInput.FileName]);
    sFileB = bst_history('add', sFileB, 'pthresh', ['Data file: ' sInput2.FileName]);
    % Output filename
    [fPath, fBase, fExt] = bst_fileparts(file_fullpath(sInput2.FileName));
    DataFile = file_unique(fullfile(fPath, [fBase, '_pthresh', fExt]));
    % Save on disk
    bst_save(DataFile, sFileB, 'v6');
    % Register in database
    db_add_data(sInput.iStudy, DataFile, sFileB);
    % Return data file
    OutputFiles{1} = DataFile;
end

author: brainstorm-tools


text:
% Defines a parameter
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 11/17/2009
classdef Parameter < edu.stanford.covert.cell.kb.KnowledgeBaseObject
    properties
        process          = edu.stanford.covert.cell.kb.Process.empty(0, 0);
        state           = edu.stanford.covert.cell.kb.State.empty(0, 0);
        reactions       = edu.stanford.covert.cell.kb.Reaction.empty(0, 0);
        proteinMonomers = edu.stanford.covert.cell.kb.ProteinMonomer.empty(0, 0);
        proteinComplexs = edu.stanford.covert.cell.kb.ProteinComplex.empty(0, 0);
    end

    properties %(SetAccess = protected)
        index
        defaultValue
        units
        experimentallyConstrained
    end

    methods
        function this = Parameter(knowledgeBase, wid, wholeCellModelID, name, ...
                index, defaultValue, units, experimentallyConstrained, ...
                comments, crossReferences)

            if nargin == 0; return; end;

            this = edu.stanford.covert.cell.kb.Parameter.empty(size(wid, 1), 0);
            this(size(wid, 1), 1) = edu.stanford.covert.cell.kb.Parameter;
            for i = 1:size(wid, 1)
                this(i, 1).idx = i;
                this(i, 1).knowledgeBase = knowledgeBase;
                this(i, 1).wid = wid(i);
                this(i, 1).wholeCellModelID = wholeCellModelID{i};
                this(i, 1).name = name{i};
                this(i, 1).index = index{i};
                this(i, 1).defaultValue = this.parseDefaultValue(defaultValue{i});
                this(i, 1).units = units{i};
                this(i, 1).experimentallyConstrained = experimentallyConstrained(i);
                if exist('comments','var') && ~isempty(comments); this(i,1).comments = comments{i}; end;
                if exist('crossReferences', 'var')
                    if size(crossReferences, 1) > 1
                        this(i, 1).crossReferences = crossReferences(i);
                    else
                        this(i, 1).crossReferences = struct;
                        fields = fieldnames(crossReferences);
                        for j = 1:size(fields,1)
                            values = crossReferences.(fields{j});
                            this(i, 1).crossReferences.(fields{j}) = values(i);
                        end
                    end
                end
            end
        end

        function serializeLinks(this)
            for i = 1:numel(this)
                this(i).process         = this.serializeLinksHelper(this(i).process);
                this(i).state           = this.serializeLinksHelper(this(i).state);
                this(i).reactions       = this.serializeLinksHelper(this(i).reactions);
                this(i).proteinMonomers = this.serializeLinksHelper(this(i).proteinMonomers);
                this(i).proteinComplexs = this.serializeLinksHelper(this(i).proteinComplexs);
                serializeLinks@edu.stanford.covert.cell.kb.KnowledgeBaseObject(this(i))
            end
        end

        function deserializeLinks(this, kb)
            for i = 1:numel(this)
                this(i).process         = this.deserializeLinksHelper(this(i).process, kb.processes);
                this(i).state           = this.deserializeLinksHelper(this(i).state, kb.states);
                this(i).reactions       = this.deserializeLinksHelper(this(i).reactions, kb.reactions);
                this(i).proteinMonomers = this.deserializeLinksHelper(this(i).proteinMonomers, kb.proteinMonomers);
                this(i).proteinComplexs = this.deserializeLinksHelper(this(i).proteinComplexs, kb.proteinComplexs);
                deserializeLinks@edu.stanford.covert.cell.kb.KnowledgeBaseObject(this(i), kb);
            end
        end

        function deleteLinks(this)
            for i = 1:numel(this)
                this(i).process         = [];
                this(i).state           = [];
                this(i).reactions       = [];
                this(i).proteinMonomers = [];
                this(i).proteinComplexs = [];
                deleteLinks@edu.stanford.covert.cell.kb.KnowledgeBaseObject(this(i))
            end
        end
    end

    methods (Static = true)
        function parsedValue = parseDefaultValue(rawValue)
            if isempty(rawValue)
                parsedValue = {};
                return;
            end

            if iscell(rawValue) && numel(rawValue) == 1
                rawValue = rawValue{1};
            end

            rawValue = strsplit(';', rawValue)';
            parsedValue = zeros(size(rawValue));
            for i = 1:length(rawValue)
                parsedValue(i) = str2double(rawValue{i});
            end
            if any(isnan(parsedValue))
                parsedValue = rawValue;
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function title = p08_title ( )

%*****************************************************************************80
%
%% P08_TITLE returns the title of problem P08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Problem 8, Enright and Pryce #B3';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = drawAngleBetweenVectors3d(o, v1, v2, r, varargin)
%DRAWANGLEBETWEENVECTORS3D Draw an arc between 2 vectors.
%
%   drawAngleBetweenVectors3d(ORIGIN, VECTOR1, VECTOR2, RADIUS)
%   draws the arc between VECTOR1 and VECTOR2.
%
%   drawAngleBetweenVectors3d(...,'ConjugateAngle',1) draws the conjugate
%   angle instead of the small angle. Default is false.
%
%   H = drawAngleBetweenVectors3d(...)
%   returns the handle of the created LINE object
%
%   Example
%     o=-100 + 200*rand(1,3);
%     v1=normalizeVector3d(-1 + 2*rand(1,3));
%     v2=normalizeVector3d(-1 + 2*rand(1,3));
%     r = rand;
%     figure('color','w'); view(3)
%     hold on; axis equal tight; xlabel X; ylabel Y; zlabel Z;
%     drawVector3d(o, v1, 'r')
%     drawVector3d(o, v2, 'g')
%     drawAngleBetweenVectors3d(o, v1, v2, r,'Color','m','LineWidth', 3)
%
%   See also
%     drawCircleArc3d

% ------
% Author: oqilipo
% E-mail: N/A
% Created: 2020-02-02
% Copyright 2020-2022

% parse axis handle
hAx = gca;
if isAxisHandle(o)
    hAx = o;
    o = v1;
    v1 = v2;
    v2 = r;
    r = varargin{1};
    varargin(1) = [];
end

p = inputParser;
p.KeepUnmatched = true;
logParValidFunc=@(x) (islogical(x) || isequal(x,1) || isequal(x,0));
addParameter(p,'ConjugateAngle',false,logParValidFunc);
parse(p, varargin{:});
conjugate=p.Results.ConjugateAngle;
drawOptions=p.Unmatched;

% Normal of the two vectors
normal=normalizeVector3d(crossProduct3d(v1, v2));
% Align normal with the z axis.
ROT = createRotationVector3d(normal,[0 0 1]);
% Align first vector with x axis
ROTv1 = createRotationVector3d(transformVector3d(v1,ROT),[1 0 0]);
% Get Euler angles of the arc.
% The arc is an flat object. Hence, use the 'ZYZ' convention.
[PHI, THETA, PSI] = rotation3dToEulerAngles((ROTv1*ROT)','ZYZ');
% Get angle between the vectors
angle=rad2deg(vectorAngle3d(v1, v2));
% Draw the arc
if ~conjugate
    h = drawCircleArc3d(hAx, [o r [THETA PHI PSI] 0 angle],drawOptions);
else
    h = drawCircleArc3d(hAx, [o r [THETA PHI PSI] 0 angle-360],drawOptions);
end

% Format output
if nargout > 0
    varargout{1} = h;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ cx, cy ] = csrot ( n, cx, incx, cy, incy, c, s )

%*****************************************************************************80
%
%% CSROT applies a complex plane rotation.
%
%  Discussion:
%
%    The cosine C and sine S are real and the vectors CX and CY are complex.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vectors.
%
%    Input, complex CX(*), one of the vectors to be rotated.
%
%    Input, integer INCX, the increment between successive entries of CX.
%
%    Input, complex CY(*), one of the vectors to be rotated.
%
%    Input, integer INCY, the increment between successive elements of CY.
%
%    Input, real C, S, parameters (presumably the cosine and sine of
%    some angle) that define a plane rotation.
%
%    Output, complex CX(*), one of the vectors to be rotated.
%
%    Output, complex CY(*), one of the vectors to be rotated.
%
  temp(1:n)               =   c * cx(1:incx:1+(n-1)*incx) ...
                            + s * cy(1:incy:1+(n-1)*incy);

  cy(1:incy:1+(n-1)*incy) = - s * cx(1:incx:1+(n-1)*incx) ...
                            + c * cy(1:incy:1+(n-1)*incy);

  cx(1:incx:1+(n-1)*incx) = temp(1:n);

  return
end


text:
function y=PhaseDemod(x, bits)
%function y=PhaseDemod(x, bits)
%
% This function demodulates M-PSK with M=2^bits symbols
% This function should be used with PhaseMod.m

M = 2^bits;

if size(x,1) ==1
   x = x.';
end

ss = 0:M-1;
symbols = exp(j*2*pi*ss/M);
b = dec2bin(0:M-1);

for i=1:M
   e(i,:) = symbols(i)-x.';
end

[min,index] = min(abs(e));

y = b(index,:);

y = reshape(y, bits*size(y,1),1);

y = str2num(y).';


text:
function model = vargplvmModelInit(model, globalOpt)

% VARGPLVMMODELINIT Initialise a vargplvm model given global demo options
% COPYRIGHT: Andreas C. Damianou, 2012
% SEEALSO: vargplvmCreate
% VARGPLVM

%if globalOpt.DgtN
   % model.mOrig = model.m;
%    model = vargplvmParamInit(model, model.mOrig, model.X);
%else
    model = vargplvmParamInit(model, model.y, model.X);
%end

% Lengthscales
if strcmp(model.kern.type, 'rbfardjit')
    model.kern.inputScales = globalOpt.invWidthMult./(((max(model.X)-min(model.X))).^2); % Default 5
else
    model.kern.comp{1}.inputScales = globalOpt.invWidthMult./(((max(model.X)-min(model.X))).^2); % Default 5
end

model.vardist.covars = 0.5*ones(size(model.vardist.covars)) + 0.001*randn(size(model.vardist.covars));
%model.kern.comp{1}.variance = max(var(Y)); %%%

if model.N > 50 && globalOpt.enableParallelism
    fprintf('# Parallel computations w.r.t the datapoints!\n');
    model.vardist.parallel = 1;
end

%%%
if ~isfield(globalOpt, 'beta') || isempty(globalOpt.betaInit)
    if model.DgtN
        model.beta = 1/((1/globalOpt.initSNR * var(model.mOrig(:))));
    else
        model.beta = 1/((1/globalOpt.initSNR * var(model.m(:))));
    end
else
    model.beta = globalOpt.betaInit;
end

if model.beta < 1e-7
    warning('Beta was too small... Setting beta to 1e-7')
    model.beta = 1e-7;
elseif model.beta > 1e+7
    warning('Beta was too big... Setting beta to 1e+7')
    model.beta = 1e+7;
end

% Check for inconsistencies
if isfield(model, 'dynamics') && ~isempty(model.dynamics)
    if isfield(model.dynamics, 'reoptimise') && model.dynamics.reoptimise
        if isfield(model, 'fixInducing') && model.fixInducing
            msg = sprintf('Reoptimising inducing points for test and fixing ind. points are incompatible!\nSetting reoptimising to false.');
            warning(msg);
            model.dynamics.reoptimise = false;
            globalOpt.testReoptimise = false;
        end
    end
end

model.dataSetInfo.dataSetName = globalOpt.dataSetName;

params = vargplvmExtractParam(model);
model = vargplvmExpandParam(model, params);

model.date = date;
model.globalOpt = globalOpt;

author: SheffieldML

repo: vargplvm


text:
% Focussed 2D Array with Directional Elements
%
% This example demonstrates the use of k-Wave to compute the outputs from a
% curved detector array which consists of several elements, each of which
% consists of a number of grid points.
%
% author: Ben Cox
% date: 29th October 2010
% last update: 22nd February 2012
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

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

clear all;

% =========================================================================
% SETUP THE GRID
% =========================================================================

% create the computational grid
Nx = 180;           % number of grid points in the x (row) direction
Ny = 180;           % number of grid points in the y (column) direction
dx = 0.1e-3;        % grid point spacing in the x direction [m]
dy = 0.1e-3;        % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500;    % [m/s]

% define the array of time points [s]
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
Nt = round(0.75*length(kgrid.t_array));
kgrid.t_array = (0:Nt-1)*dt;

% =========================================================================
% DEFINE A FOCUSSED ARRAY OF DIRECTIONAL ELEMENTS
% =========================================================================

% define a semicircular sensor centered on the grid
semicircle_radius = 65; % [grid points]
arc = makeCircle(Nx, Ny, Nx/2, Ny/2, semicircle_radius, pi);

% find total number and indices of the grid points constituting the
% semicircle
arc_indices = find(arc == 1);
Nv = length(arc_indices);

% calculate angles between grid points in the arc and the centre of the
% grid
arc_angles = atan((kgrid.y(arc_indices))./kgrid.x(arc_indices));

% sort the angles into ascending order, and adjust the indices accordingly
[sorted_arc_angles,sorted_index] = sort(arc_angles);
sorted_arc_indices = arc_indices(sorted_index);

% divide the semicircle into Ne separate sensor elements
Ne = 13;
sensor.mask = zeros(Nx,Ny);
for loop = 1:Ne

    % the indices of the grid points belonging to one element.
    % (There is a two grid point gap between the elements.)
    voxel_indices = sorted_arc_indices(floor((loop-1)*Nv/Ne)+2:floor(loop*Nv/Ne)-1);

    % add the element to the sensor.mask
    sensor.mask(voxel_indices) = 1;

end

% =========================================================================
% SIMULATION
% =========================================================================

% Define an infinitely wide plane wave source. (This is achieved by turning
% off the PML.)
source.p_mask = zeros(Nx,Ny);
source.p_mask(140,:) = 1;

% set the display mask for the simulation
display_mask = source.p_mask + sensor.mask;

% define a time varying sinusoidal source
source_freq = 1e6;
source_mag = 0.5;
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);
source.p = filterTimeSeries(kgrid, medium, source.p);

% run the simulation with the PML switched off on the sides, so that the
% source continues up to the edge of the domain (and from there infinitely,
% because of the periodic assumption implicit in pseudospectral methods).
input_args = {'PMLAlpha', [2 0], 'DisplayMask', display_mask, 'PlotScale', [-0.75 0.75]};
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% split up the data, recorded on all the grid points, between the elements
element_data = zeros(Ne,Nt);
for loop = 1:Ne

    % the indices of the sensor grid points in the sensor mask
    sensor_indices = find(sensor.mask==1);

    % the indices of the grid points belonging to one element.
    voxel_indices = sorted_arc_indices(floor((loop-1)*Nv/Ne)+2:floor(loop*Nv/Ne)-1);

    % indices of sensor_data that refer to the data for this element
    data_indices = zeros(length(voxel_indices),1);
    for loop2 = 1:length(voxel_indices)
        data_indices(loop2) = find(sensor_indices == voxel_indices(loop2));
    end

    % for one element per loop, average the time series from each of the
    % element's grid points to give one time series for the whole element
    element_data(loop,:) = mean(sensor_data(data_indices,:),1);

end

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the elements of the array, and the source mask
figure;
imagesc(display_mask);
axis image;
colormap(flipud(gray));

% plot the time series recorded at each array element
figure;
stackedPlot(kgrid.t_array*1e6, element_data);
xlabel('time [\mus]');
ylabel('time series, one for each element');

repo: Aero-matlab


text:
function [TRI,xout,yout,uout,interp] = PlotField2D(Nout, xin, yin, uin)

% function [TRI,xout,yout,uout,interp] = PlotField2D(Nout, xin, yin, uin)
% Purpose: filled contour plot of solution data

Globals2D;

% build equally spaced grid on reference triangle
Npout = (Nout+1)*(Nout+2)/2;
rout = zeros(Npout,1); sout = zeros(Npout,1);
sk = 1;
for n=1:Nout+1
  for m=1:Nout+2-n
    rout(sk) = -1 + 2*(m-1)/Nout;
    sout(sk) = -1 + 2*(n-1)/Nout;
    counter(n,m) = sk; sk = sk+1;
  end
end

% build matrix to interpolate field data to equally spaced nodes
interp = InterpMatrix2D(rout, sout);

% build triangulation of equally spaced nodes on reference triangle
tri = [];
for n=1:Nout+1
  for m=1:Nout+1-n,
    v1 = counter(n,m);   v2 = counter(n,m+1);
    v3 = counter(n+1,m); v4 = counter(n+1,m+1);
    if(v4)
      tri = [tri;[[v1 v2 v3];[v2 v4 v3]]];
    else
      tri = [tri;[[v1 v2 v3]]];
    end
  end
end

% build triangulation for all equally spaced nodes on all elements
TRI = [];
for k=1:K
  TRI = [TRI; tri+(k-1)*Npout];
end

% interpolate node coordinates and field to equally spaced nodes
xout = interp*xin; yout = interp*yin; uout = interp*uin;

% render and format solution field
trisurf(TRI, xout(:), yout(:), uout(:));
shading interp,    material shiny,    lighting gouraud
camlight headlight
return


text:
% test for propagation and geodesic extraction on 2D planar shape

path(path, 'toolbox/');
path(path, 'data/');

name = 'chicken';
name = 'apple';
name = 'cavern';
name = 'camel';
name = 'giraffe';

rep = 'results/shape-geodesics/';
if not(exist(rep))
    mkdir(rep);
end

n = 128;
M = rescale( load_image(name,n), 0,1 );
M = perform_blurring(M,5);
M = double(M>0.5);

% make sure pixels on the boundary are black
if M(1)==1
    M = 1-M;
end

warning off;
imwrite(1-M, [rep name '-shape.png'], 'png');
warning off;

% compute geodesic distance
clf;
imagesc(M); axis image; axis off;
title('click on a point inside the shape');
[y,x] = ginput(1);
start_points = round([x y]');
W = ones(n);
L = zeros(n)-Inf; L(M==1) = +Inf;
options.constraint_map = L;
disp('Compute distance function');
[D,S,Q] = perform_fast_marching(W, start_points, options);

bound = compute_shape_boundary(M);
nbound = size(bound,1);
npaths = 30;
sel = round(linspace(1,nbound+1,npaths+1)); sel(end) = [];
end_points = bound(sel,:);

disp('Extract paths');
paths = {};
D1 = D; D1(M==0) = 1e9;
for i=1:npaths
    paths{i} = compute_geodesic(D1,end_points(i,:)');
%    paths{i} = compute_discrete_geodesic(D1,end_points(i,:)')';
end

ms = 30; lw = 3;
% display
A = convert_distance_color(D);
clf; hold on;
imageplot(A); axis image; axis off;
for i=1:npaths
    end_point = end_points(i,:);
    h = plot( paths{i}(2,:), paths{i}(1,:), 'k' );
    set(h, 'LineWidth', lw);
    h = plot(end_point(2),end_point(1), '.b');
    set(h, 'MarkerSize', ms);
end
h = plot(start_points(2),start_points(1), '.r');
set(h, 'MarkerSize', ms);
hold off;
colormap jet(256);
axis ij;
saveas(gcf, [rep name '-geodesics.png'], 'png');

author: alecjacobson

repo: gptoolbox


text:
slice = function mrThickOblVol(thick,obX,obY,volume,sagSize,numSlices,reflections,aTheta,cTheta,curSag);
% OBSOLETE
%
% mrThickOblVol(thick,obX,obY,volume,sagSize,numSlices,reflections
%                                    ,aTheta,cTheta,curSag);
%
% Places an averaged oblique image in the oblique window.
% POIRSON 08.09.96 Added reflection logic
% SPG 12.15.96 Replaced old oblique slice function calls. removed
%   reflection calls. they are covered by mrRotSagVol.

global obwin volslimin2 volslimax2

numextracts = 10;

if(obX(1) <= sagSize(2) & obY(1) <= sagSize(1) & obX(1) > 0 & obY(1) > 0 & ...
   obX(2) <= sagSize(2) & obY(2) <= sagSize(1) & obX(2) > 0 & obY(2) > 0 )

  d = sqrt((obY(2)-obY(1)).^2 + (obX(2)-obX(1)).^2);
  unitv = [(obX(2)-obX(1))/d, (obY(2)-obY(1))/d];
  perp = [-unitv(2), unitv(1)];
  obSize = mrFindObSize(obX,obY,sagSize,numSlices);

  newX = obX-perp(1)*(thick/2);
  newY = obY-perp(2)*(thick/2);

  for i = 1:numextracts
	newX = newX+perp(1)*(thick/numextracts);
	newY = newY+perp(2)*(thick/numextracts);
	[sagSlice,sagPts,temp,obPts] = mrRotSagVol(volume,newX,newY,obSize,sagSize,cTheta,aTheta,curSag,reflections,0);
        obSlices(i,:) = temp;
  end
	slice = mean(obSlices);

else

  error('Oblique points out of range. Re-Clip inplanes grid.');
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [prob_temp, t_templ,t_temph] = problTemperature(BandBT,idclr,l_pt,h_pt)
%PROBTEMPERATURE calculate temperature probability for land respectively.

    % [Temperature test (over land)]
    F_temp=BandBT(idclr); % get clear temperature
    clear idclr;
    t_buffer=4*100;
    % 0.175 percentile background temperature (low)
    t_templ=prctile(F_temp,100*l_pt);
    % 0.825 percentile background temperature (high)
    t_temph=prctile(F_temp,100*h_pt);
    clear F_temp l_pt h_pt;

    t_tempL=t_templ-t_buffer;
    t_tempH=t_temph+t_buffer;
    clear t_buffer;
    Temp_l=t_tempH-t_tempL;
    clear t_tempL;
    prob_temp=(t_tempH-BandBT)/Temp_l;
    clear BandBT t_tempH Temp_l;
    % Temperature can have prob > 1
    prob_temp(prob_temp<0)=0;
%     prob_temp(prob_temp>1)=1;

end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
function warning_flexible(varargin)
% Function allowing to have better control over the warnings
% that might not be necessary at some point
% _______________________________________________________________________
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: warning_flexible.m 5467 2013-05-05 20:03:50Z vladimir $

warning off backtrace
warning(varargin{:});
warning on backtrace

author: spm


text:
function [grad_scores] = normLayerBackprop(grad_alignWeights, alignWeights) %, maskedIds, params)
  %% from grad_alignWeights -> grad_scores
  % alignWeights a = softmax(scores)
  % Let's derive per indices grad align weight w.r.t scores
  %   der a_i / der s_j = der exp(s_i) / sum_k (exp(s_k)) / der s_j =
  %     (1/sum) * (der exp(s_i) / der s_j) - (exp(s_i)/sum^2)*exp(s_j) =
  %      a_i*I{i==j} - a_i*_a_j
  %
  % Now let's try to optimize the vector grad for a single example i:
  %   grad_score_i = (diag(a_i) - a_i*a_i')*grad_a_i
  %                = a_i.*grad_a_i - a_i*(a_i'*grad_a_i)
  %                = a_i.*grad_a_i - a_i*alpha_i
  % multiple examples: alpha = sum(a.*grad_a, 1) % 1*curBatchSize
  %     grad_scores = a.*grad - bsxfun(@times, a, alpha)
  % tmpResult = alignWeights.*grad_alignWeights; % numAttnPositions * curBatchSize

  tmpResult = alignWeights.*grad_alignWeights; % numAttnPositions * curBatchSize
  grad_scores = tmpResult - bsxfun(@times, alignWeights, sum(tmpResult, 1));

%   % assert
%   if params.assert
%     % compute grad_scores in a different way
%     grad_scores1 = zeroMatrix(size(grad_scores), params.isGPU, params.dataType);
%     for ii=1:params.curBatchSize
%       grad_scores1(:, ii) = (diag(alignWeights(:, ii))-alignWeights(:, ii)*alignWeights(:, ii)')*grad_alignWeights(:, ii);
%     end
%     assert(computeSum(grad_scores-grad_scores1, params.isGPU)<1e-5);
%
%     assert(computeSum(alignWeights(maskedIds), params.isGPU)==0);
%     assert(computeSum(grad_scores(maskedIds), params.isGPU)==0);
%     assert(computeSum(tmpResult(maskedIds), params.isGPU)==0);
%     tmpResult = bsxfun(@times, alignWeights, sum(tmpResult, 1));
%     assert(computeSum(tmpResult(maskedIds), params.isGPU)==0);
%   end
end

author: lmthang


text:
function test_qrsol
%TEST_QRSOL test cs_qrsol
%
% Example:
%   test_qrsol
% See also: testall

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

index = UFget ;
[ignore f] = sort (max (index.nrows, index.ncols)) ;

k = 0 ;
rs1 = zeros (1,0) ;
rs2 = zeros (1,0) ;

for i = f
    Prob = UFget (i,index) ;
    A = Prob.A ;
    if (~isreal (A))
        continue ;
    end

    [m n] = size (A) ;                                                  %#ok
    b = rand (m,1) ;

    x1 = A\b ;
    x2 = cs_qrsol (A,b) ;

    x1 (~isfinite (x1)) = 0 ;
    x2 (~isfinite (x2)) = 0 ;

    r1 = norm (A*x1-b) ;
    r2 = norm (A*x2-b) ;

    k = k + 1 ;
    rs1 (k) = r1 ;
    rs2 (k) = r2 ;

    fprintf ('%30s  MATLAB: %6.2e CS: %6.2e\n', Prob.name, r1, r2) ;

    loglog (rs1, rs2, 'o') ;
    drawnow

    clear A b x1 x2
    % pack

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function test11
%TEST11 test cs_rowcnt
%
% Example:
%   test11
% See also: testall

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

clear functions
index = UFget ;
[ignore f] = sort (max (index.nrows, index.ncols)) ;
f = f (1:200) ;

for i = f
    Prob = UFget (i, index) ;
    disp (Prob) ;
    A = Prob.A ;
    [m n] = size (A) ;
    if (~isreal (A) | m ~= n)                                               %#ok
        continue
    end

    A = spones (A) ;
    A = A+A' + speye(n) ;

    [cc h pa po R] = symbfact (A) ;
    rc1 = full (sum (R)) ;
    rc2 = cs_rowcnt (A, pa, po) ;
    if (any (rc1 ~= rc2))
        error ('!') ;
    end

    try
        p = amd (A) ;
    catch
        p = symamd (A) ;
    end
    A = A (p,p) ;

    [cc h pa po R] = symbfact (A) ;
    rc1 = full (sum (R)) ;
    rc2 = cs_rowcnt (A, pa, po) ;
    if (any (rc1 ~= rc2))
        error ('!') ;
    end

end


text:
function Offspring = Operator(Problem,Population)
% Differential evolution in FROFI

%------------------------------- 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".
%--------------------------------------------------------------------------

    %% Parameter setting
    PopDec = Population.decs;
    [N,D]  = size(PopDec);
    CR = [0.1 0.2 1]';
    CR = repmat(CR(randi(end,N,1)),1,D);
    F  = [0.6 0.8 1]';
    F  = repmat(F(randi(end,N,1)),1,D);

    %% Parents
    [~,P] = sort(rand(N),2);
    P1    = PopDec(P(:,1),:);
    P2    = PopDec(P(:,2),:);
    P3    = PopDec(P(:,3),:);
    [~,B] = min(Population.objs);
    PB    = repmat(PopDec(B,:),N,1);

    %% Offspring generation
    Rand   = rand(N,D);
    k1     = repmat(rand(N,1)<0.5,1,D);
    k2     = ~k1 & rand(N,D)<CR;
    OffDec = PopDec;
    OffDec(k1) = PopDec(k1) + Rand(k1).*(P1(k1)-PopDec(k1)) + F(k1).*(P2(k1)-P3(k1));
    OffDec(k2) = P1(k2) + Rand(k2).*(PB(k2)-P1(k2)) + F(k2).*(P2(k2)-P3(k2));
    Offspring  = Problem.Evaluation(OffDec);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function model = yalmip2intlinprog(interfacedata)

options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
Q       = interfacedata.Q;
lb      = interfacedata.lb;
ub      = interfacedata.ub;

ops = options.intlinprog;
switch options.verbose
    case 0
        ops.Display = 'off';
    case 1
        ops.Display = 'final';
    otherwise
        ops.Display = 'iter';
end

if isempty(F_struc)
    Aeq = [];
    beq = [];
    A = [];
    b = [];
else
    Aeq = -F_struc(1:1:K.f,2:end);
    beq = F_struc(1:1:K.f,1);
    A =-F_struc(K.f+1:end,2:end);
    b = F_struc(K.f+1:end,1);
end

intcon = union(interfacedata.integer_variables,interfacedata.binary_variables);
if ~isempty(interfacedata.binary_variables)
    lb(interfacedata.binary_variables) = max(0,lb(interfacedata.binary_variables));
    ub(interfacedata.binary_variables) = min(1,ub(interfacedata.binary_variables));
end

model.c = c;
model.intcon = intcon;
model.A = A;
model.b = b;
model.Aeq = Aeq;
model.beq = beq;
model.lb = lb;
model.ub = ub;
model.ops = ops;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef AtmoWithNoDragModelValidator < AbstractLaunchVehicleDataValidator
    %AtmoWithNoDragModelValidator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        lvdData LvdData
    end

    methods
        function obj = AtmoWithNoDragModelValidator(lvdData)
            obj.lvdData = lvdData;
        end

        function [errors, warnings] = validate(obj)
            errors = LaunchVehicleDataValidationError.empty(0,1);
            warnings = LaunchVehicleDataValidationWarning.empty(0,1);

            warnEvtNums = [];
            evts = obj.lvdData.script.evts;
            for(i=1:length(evts))
                evt = evts(i);
                evtPropagator = evt.propagatorObj;
                propagatorEnum = evtPropagator.propagatorEnum;

                if(propagatorEnum == PropagatorEnum.ForceModel)
                    if(not(any(evtPropagator.forceModels == ForceModelsEnum.Drag)))
                        stateLogEntries = obj.lvdData.stateLog.getAllStateLogEntriesForEvent(evt);

                        for(j=1:length(stateLogEntries))
                            stateLogEntry = stateLogEntries(j);
                            bodyInfo = stateLogEntry.centralBody;
                            altitude = norm(stateLogEntry.position) - bodyInfo.radius;

                            if(altitude <= bodyInfo.atmohgt && altitude >= 0)
                                warnEvtNums(end+1) = evt.getEventNum(); %#ok<AGROW>
                                break;
                            end
                        end
                    end
                elseif(propagatorEnum == PropagatorEnum.TwoBody)
                    stateLogEntries = obj.lvdData.stateLog.getAllStateLogEntriesForEvent(evt);

                    for(j=1:length(stateLogEntries))
                        stateLogEntry = stateLogEntries(j);
                        bodyInfo = stateLogEntry.centralBody;
                        altitude = norm(stateLogEntry.position) - bodyInfo.radius;

                        if(altitude <= bodyInfo.atmohgt && altitude >= 0)
                            warnEvtNums(end+1) = evt.getEventNum(); %#ok<AGROW>
                            break;
                        end
                    end
                end
            end

            if(not(isempty(warnEvtNums)))
                eventStr = makeEventsStr(unique(warnEvtNums));
                str = sprintf('Drag model is disabled on events that exist in an atmosphere. (Events: %s)\n%s', eventStr);
                warnings(end+1) = LaunchVehicleDataValidationWarning(str);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function cost = costFun(u)
% cost = costFun(u)
%
% Cost is the integral of torque-squared.

cost = u.^2;

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function foundHandles = find(obj, nameClass, varargin)
% Finds all (handle) object properties whose class and property
% values match given values
%
%   foundHandles = obj.find(nameClass, 'PropertyName',
%                           'PropertyValue', ...)
%
%
% IN
%   nameClass       string with class name (default: MrCopyData)
%   PropertyName/   pairs of string containing name of property
%   PropertyValue   and value (or pattern) that has to be matched
%                   NOTE: if cells of values are given, all
%                   objects are returned that match any of the
%                   entries,
%                   e.g. 'name', {'mean', 'snr'} will return
%                   objects if they are named 'mean OR 'snr'
%
% OUT
%   foundHandles    cell(nHandles,1) of all object handles for
%                   objects that match the properties
%
%                   NOTE: function also returns handle to
%                   calling object, if its class and properties
%                   fulfill the given criteria
%
% EXAMPLE:
%   Y = MrCopyData();
%   Y.find('MrCopyData', 'name', 'coolCopy');
%
%   See also MrCopyData

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2016-04-19
% Copyright (C) 2016 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

foundHandles = {};
if nargin
    searchPropertyNames = varargin(1:2:end);
    searchPropertyValues = varargin(2:2:end);
    nSearchProperties = numel(searchPropertyNames);

else
    nSearchProperties = 0;
end

% check whether object itself fulfills criteria
if isa(obj, nameClass)

    doesMatchProperties = true;
    iSearchProperty = 1;

    % Check search properties as long as matching to values
    while doesMatchProperties && ...
            iSearchProperty <= nSearchProperties

        searchProperty = searchPropertyNames{iSearchProperty};
        searchValue = searchPropertyValues{iSearchProperty};
        if isa(obj.(searchProperty), 'MrCopyData')
            % recursive comparison for MrCopyData-properties
            doesMatchProperties = obj.(searchProperty).comp(...
                searchValue);
        else

            doesMatchProperties = isequal(obj.(searchProperty), ...
                searchValue);

            % allow pattern matching for strings or cell of
            % strings (matching any entry of cell)
            if ischar(obj.(searchProperty))
                cellSearchValue = cellstr(searchValue);
                nCellEntries = numel(cellSearchValue);
                iCellEntry = 1;

                % check for each entry in cell whether it
                % matches the string value of this object's
                % property
                while iCellEntry <= nCellEntries && ~doesMatchProperties
                    currentSearchValue = ...
                        cellSearchValue{iCellEntry};

                    % check whether pattern expression given, i.e.
                    % * in search value
                    isSearchPattern = ~isempty(strfind(currentSearchValue, ...
                        '*'));
                    if isSearchPattern
                        doesMatchProperties = ~isempty(regexp( ...
                            obj.(searchProperty), currentSearchValue, 'once'));
                    else
                        doesMatchProperties = isequal(obj.(searchProperty), ...
                            currentSearchValue);
                    end
                    iCellEntry = iCellEntry + 1;
                end
            end

        end

        iSearchProperty = iSearchProperty + 1;
    end

    if doesMatchProperties
        foundHandles = [foundHandles; {obj}];
    end
end

[sel, mobj] = get_properties_to_update(obj);
for k = sel(:)'
    pname = mobj.Properties{k}.Name;

    % Continue to check properties recursively for MrCopyData-properties
    if isa(obj.(pname), 'MrCopyData') % recursive comparison
        newFoundHandles = obj.(pname).find(nameClass, varargin{:});
        foundHandles = [foundHandles;newFoundHandles];
    else
        % cell of MrCopyData also treated
        if iscell(obj.(pname)) && length(obj.(pname)) ...
                && isa(obj.(pname){1}, 'MrCopyData')
            for c = 1:length(obj.(pname))
                newFoundHandles = obj.(pname){c}.find(nameClass, varargin{:});
                foundHandles = [foundHandles;newFoundHandles];
            end
        end
    end
end
% remove duplicate entries, if sub-object itself was returned and
% as a property of super-object
% foundHandles = unique(foundHandles);
end

author: translationalneuromodeling

repo: tapas


text:
function hermite_test01 ( )

%*****************************************************************************80
%
%% TEST01 uses f(x) = 1 + 2x + 3x^2 at x = 0, 1, 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2011
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  x =  [ 0.0; 1.0;  2.0 ];
  y =  [ 1.0; 6.0; 17.0 ];
  yp = [ 2.0; 8.0; 14.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  HERMITE computes the Hermite interpolant to data.\n' );
  fprintf ( 1, '  Here, f(x) = 1 + 2x + 3x^2.\n' );

  hermite_demo ( n, x, y, yp );

  return
end


text:
function grad = compute_NN_grad(W, X, resp, hout, nonlinearity, nonzero_grad)

% nonlinearity = 1 -- hyperbolic tangent
% nonlinearity = 2 -- logistic sigmoid (subtract 0.5 from last layer's outputs)
% nonlinearity = 3 -- logistic sigmoid + hyperbolic tangent (last layer)

ncases = size(resp{1},2);
if isa(X, 'gsingle')
  gputype = 'gsingle';
  onesncases = gones(ncases,1,gputype);
elseif isa(X, 'gdouble')
  gputype = 'gsingle';
  onesncases = gones(ncases,1,gputype);
else
  onesncases = ones(ncases,1);
end

nlayer = numel(W);

if (nonlinearity == 1)

  if (exist('nonzero_grad', 'var'))
    snonzero_grad = sum(nonzero_grad);
    if (snonzero_grad == 0)
      for i=nlayer:-1:1
	grad{i} = zeros(size(W{i}));
      end
      return;
    end

    for i=nlayer:-1:1
      if (i == nlayer)
	tmp = hout(:, nonzero_grad);
      else
	tmp = (W{i+1}' * backward{i+1});
      end

      if (size(tmp, 1) ~= size(resp{i}, 1))
	tmp = tmp(1:end-1, :);
      end

      backward{i} = tmp .* (1 - resp{i}(:,nonzero_grad).^2);

      if (i-1) >= 1
	if (size(W{i}, 2) == size(resp{i-1}, 1) + 1)
	  grad{i} = backward{i} * [resp{i-1}(:,nonzero_grad)' ones(snonzero_grad, 1)];
	else
	  grad{i} = backward{i} * resp{i-1}(:,nonzero_grad)';
	end
      else
	grad{i} = backward{i} * X(:, nonzero_grad)';
      end
    end
  else
    for i=nlayer:-1:1
      if (i == nlayer)
	tmp = hout;
      else
	tmp = (W{i+1}' * backward{i+1});
      end

      if (size(tmp, 1) ~= size(resp{i}, 1))
	tmp = tmp(1:end-1, :);
      end

      backward{i} = tmp .* (1 - resp{i}.^2);

      if (i-1) >= 1
	if (size(W{i}, 2) == size(resp{i-1}, 1) + 1)
	  grad{i} = backward{i} * [resp{i-1}' onesncases];
	else
	  grad{i} = backward{i} * resp{i-1}';
	end
      else
	grad{i} = backward{i} * X';
      end
    end
  end

elseif (nonlinearity == 2)

  if (exist('nonzero_grad', 'var'))
    snonzero_grad = sum(nonzero_grad);
    if (snonzero_grad == 0)
      for i=nlayer:-1:1
	grad{i} = zeros(size(W{i}));
      end
      return;
    end

    for i=nlayer:-1:1
      if (i == nlayer)
	tmp = hout(:, nonzero_grad);
      else
	tmp = (W{i+1}' * backward{i+1});
      end

      if (size(tmp, 1) ~= size(resp{i}, 1))
	tmp = tmp(1:end-1, :);
      end

      backward{i} = tmp .* (1 - resp{i}(:,nonzero_grad)) .* resp{i}(:,nonzero_grad);

      if (i-1) >= 1
	if (size(W{i}, 2) == size(resp{i-1}, 1) + 1)
	  grad{i} = backward{i} * [resp{i-1}(:,nonzero_grad)' gones(snonzero_grad,1,gputype)];
	else
	  grad{i} = backward{i} * resp{i-1}(:,nonzero_grad)';
	end
      else
	grad{i} = backward{i} * X(:, nonzero_grad)';
      end
    end
  else
    for i=nlayer:-1:1
      if (i == nlayer)
	tmp = hout;
      else
	tmp = (W{i+1}' * backward{i+1});
      end

      if (size(tmp, 1) ~= size(resp{i}, 1))
	tmp = tmp(1:end-1, :);
      end

      backward{i} = tmp .* (1 - resp{i}) .* resp{i};

      if (i-1) >= 1
	if (size(W{i}, 2) == size(resp{i-1}, 1) + 1)
	  grad{i} = backward{i} * [resp{i-1}' onesncases];
	else
	  grad{i} = backward{i} * resp{i-1}';
	end
      else
	grad{i} = backward{i} * X';
      end
    end
  end

elseif (nonlinearity == 3)

  if (exist('nonzero_grad', 'var'))
    snonzero_grad = sum(nonzero_grad);
    if (snonzero_grad == 0)
      for i=nlayer:-1:1
	grad{i} = zeros(size(W{i}));
      end
      return;
    end

    for i=nlayer:-1:1
      if (i == nlayer)
	tmp = hout(:, nonzero_grad);
      else
	tmp = (W{i+1}' * backward{i+1});
      end

      if (size(tmp, 1) ~= size(resp{i}, 1))
	tmp = tmp(1:end-1, :);
      end

      if (i == nlayer)
	backward{i} = tmp .* (1 - resp{i}(:,nonzero_grad)^2);
      else
	backward{i} = tmp .* (1 - resp{i}(:,nonzero_grad)) .* resp{i}(:,nonzero_grad);
      end

      if (i-1) >= 1
	if (size(W{i}, 2) == size(resp{i-1}, 1) + 1)
	  grad{i} = backward{i} * [resp{i-1}(:,nonzero_grad)' gones(snonzero_grad,1,gputype)];
	else
	  grad{i} = backward{i} * resp{i-1}(:,nonzero_grad)';
	end
      else
	grad{i} = backward{i} * X(:, nonzero_grad)';
      end
    end
  else
    for i=nlayer:-1:1
      if (i == nlayer)
	tmp = hout;
      else
	tmp = (W{i+1}' * backward{i+1});
      end

      if (size(tmp, 1) ~= size(resp{i}, 1))
	tmp = tmp(1:end-1, :);
      end

      if (i == nlayer)
	backward{i} = tmp .* (1 - resp{i}.^2);
      else
	backward{i} = tmp .* (1 - resp{i}) .* resp{i};
      end

      if (i-1) >= 1
	if (size(W{i}, 2) == size(resp{i-1}, 1) + 1)
	  grad{i} = backward{i} * [resp{i-1}' onesncases];
	else
	  grad{i} = backward{i} * resp{i-1}';
	end
      else
	grad{i} = backward{i} * X';
      end
    end
  end

end

author: norouzi

repo: hdml

save_path: github-repos/MATLAB/norouzi-hdml


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

function test1(testCase)
    N = 100;
    S = 100;
    gen = spx.data.noise.Basic(N, S);
    sigma = 1;
    X = gen.gaussian(sigma);
    variance = sum(sum(X.^2)) / (N * S);
    verifyEqual(testCase, sigma, variance, 'RelTol', .1);
end

author: indigits

repo: sparse-plex


text:
% PDELDEMO This program creates the left preconditionied PDE example
%          from Chapter 3
%
% The example has the form
%
% L(u) = f, 0 < x,y < 1
% with homogeneous Dirichlet BC
%
% In terms of the program below, the equation is
% Lu = - (u_xx + u_yy) + 20*u*(u_x + u_y)  = f
%
% Within the nonlinear map, pdefun.m, we apply fish2d.m as
% a preconditioner.
%
% The action of the laplacian, partial wrt x, and partial wrt y
% are computed in the routines lapmf, dxmf, and dymf. Here mf
% stands for MATRIX FREE.
%
global rhsf prhsf;
m=31; n=m*m;
h=1/(m+1);
h2=(m+1)*(m+1);
h1=(m+1);
tolh=1/h2;
mres=3;
%
% set up the equation
%
x=1:m;
e=x';
x=x'*h;
%
% Form the solution: sol(x,y)= 10 x y (1-x) (1-y) exp(x^4.5)
%
zsoly=x.*(1-x);
zsolx=zsoly.*exp(x.^4.5);
solt=10*zsolx*zsoly';
sol=solt(:);
clear rhs; clear axt; clear at; clear ayt; clear solt;

% fix the right side so that the solution is known
%
% unpreconditioned and right preconditioned
%
rhsf = lapmf(sol) + 20*sol.*(dxmf(sol)+dymf(sol));
%
% left preconditioned
%
prhsf = sol + 20*fish2d(sol.*(dxmf(sol)+dymf(sol)));
%
% Set the initial iterate to zero. Call the solver three times with
% each of the Krylov methods
%
u0=zeros(m*m,1);
%
tol=[1.d-1,1.d-1]*tolh;
%
% GMRES
%
x=zeros(n,1); parms = [40,40,.9,1];
[sol, it_histg, ierr] = nsoli(x,'pdeleft',tol,parms);
%
% BICGSTAB
%
x=zeros(n,1); parms = [40,40,.9,3];
[sol, it_histb, ierr] = nsoli(x,'pdeleft',tol,parms);
%
% TFQMR
%
x=zeros(n,1); parms = [40,40,.9,4];
[sol, it_histt, ierr] = nsoli(x,'pdeleft',tol,parms);
figure(1)
ng=length(it_histg(:,1));
nb=length(it_histb(:,1));
nt=length(it_histt(:,1));
semilogy(0:ng-1,it_histg(:,1)/it_histg(1,1),'-',...
0:nb-1,it_histb(:,1)/it_histg(1,1),'--',...
0:nt-1,it_histt(:,1)/it_histg(1,1),'-.');
xlabel('Nonlinear iterations');
ylabel('Relative residual');
legend('GMRES','BICGSTAB','TFQMR');
figure(2)
semilogy(it_histg(:,2),it_histg(:,1)/it_histg(1,1),'-',...
it_histb(:,2),it_histb(:,1)/it_histb(1,1),'--',...
it_histt(:,2),it_histt(:,1)/it_histt(1,1),'-.');
xlabel('Function evaluations');
ylabel('Relative residual');
legend('GMRES','BICGSTAB','TFQMR');


text:
function InitGUI(handles)
% Initialization function for GUI (runs once only)

% Copyright 2010 The MathWorks, Inc.
%
% Auth/Revision:
%   The MathWorks Consulting Group
%   $Author: rjackey $
%   $Revision: 241 $  $Date: 2010-01-26 15:17:00 -0500 (Tue, 26 Jan 2010) $
% -------------------------------------------------------------------------

% Get the handle to the main GUI
hFigure = handles.MainFig;

% Get the main list of models
ItemsList = getappdata(hFigure, 'ItemsList');

% Get initial selection (same as initial output)
InitSelectionIndex = getappdata(hFigure, 'OutputList');

% Get the title of the GUI
guiTitle = getappdata(handles.MainFig, 'guiTitle');

% Set the GUI to normal/modal
% set(hFigure, 'WindowStyle', 'normal');
set(hFigure, 'WindowStyle', 'modal');

% Center the GUI on the screen
movegui(hFigure, 'center')

% Set the title for the application
set(hFigure, 'Name', guiTitle);

%**********************************
%% Get some initialization variables
%**********************************

% Assign initial right side selection
ItemsIdx_R = InitSelectionIndex;

% Assign left side to all items, then empty the selected ones
ItemsIdx_L = 1:length(ItemsList);
ItemsIdx_L(ItemsIdx_R) = [];

%******************************
%% Initialize some appdatas
%******************************

% Set the index of selected models
setappdata(hFigure, 'ItemsIdx_L', ItemsIdx_L);
setappdata(hFigure, 'ItemsIdx_R', ItemsIdx_R);

%******************************************
%% Finish the rest of the GUI initializations
%******************************************

% Assign the tooltips to the GUI items
AssignTooltips(handles);

% Refresh the GUI
RefreshGUI(handles);

author: Sable

repo: mcbench-benchmarks


text:
function sphere_llt_grid_line_count_test ( )

%*****************************************************************************80
%
%% SPHERE_LLT_GRID_LINE_COUNT_TEST tests SPHERE_LLT_GRID_LINE_COUNT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2015
%
%  Author:
%
%    John Burkardt
%
  lat_num = 3;
  long_num = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_LLT_GRID_LINE_COUNT_TEST\n' );
  fprintf ( 1, '  SPHERE_LLT_GRID_LINE_COUNT counts the lines used for a\n' );
  fprintf ( 1, '  grid based on triangles defined by latitude and longitude\n' );
  fprintf ( 1, '  lines on a sphere in 3D.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     LAT_NUM    LONG_NUM   LINE_NUM\n' );
  for lat_num = 1 : 2 : 17
    fprintf ( 1, '\n' );
    long_num = 1;
    for long_log = 1 : 4
      long_num = long_num * 2;
      line_num = sphere_llt_grid_line_count ( lat_num, long_num );
      fprintf ( 1, '  %8d  %8d  %8d\n', lat_num, long_num, line_num );
    end
  end

  return
end

author: johannesgerer


text:
function fx2 = p03_fx2 ( x )

%*****************************************************************************80
%
%% P03_FX2 evaluates the second derivative of the function for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the abscissa.
%
%    Output, real FX2, the second derivative of the function at X.
%
  fx2 = exp ( - x ) * ( x - 2.0 );

  return
end


text:
function validJunctionMets = findMetabolicJunctions(model, nRxnsJnc)
% Finds metabolic branchpoints with different numbers of branches
%
% USAGE:
%
%    validJunctionMets = findMetabolicJunctions(model, nRxnsJnc)
%
% INPUTS:
%    model:                COBRA model structure
%    nRxnJnc:              Number of reactions to be considered a junction
%
% OUTPUT:
%    validJunctionMets:    List of junction metabolites
%
% .. Author: - Markus Herrgard 12/14/06

if (isfield(model,'c'))
    selRxnsC = (model.c == 0);
else
    selRxnsC = true(length(model.rxns),1);
end

%rxnGeneMat is a required field for this function, so if it does not exist,
%build it.
if ~isfield(model,'rxnGeneMat')
    model = buildRxnGeneMat(model);
end

[baseMetNames,compSymbols,uniqueMetNames,uniqueCompSymbols] = parseMetNames(model.mets);
uniqueMetNames = uniqueMetNames';
for i = 1:length(uniqueMetNames)
    sel = ismember(baseMetNames,uniqueMetNames{i});
    nRxnsMetUni(i) = sum(any(model.S(sel,selRxnsC) ~= 0,1));
end
nRxnsMetUni = full(nRxnsMetUni');
junctionMets = uniqueMetNames(nRxnsMetUni >= nRxnsJnc);

validJunctionMets = {};
for i = 1:length(junctionMets)
    sel = ismember(baseMetNames,junctionMets{i});
    if (length(unique(compSymbols(sel))) == 1)
        selRxns = any(model.S(sel,:) ~= 0,1) & selRxnsC';
        thisRxns = model.rxns(selRxns);
        geneMap = model.rxnGeneMat(findRxnIDs(model,thisRxns),:);
        selNonZero = any(geneMap,2);
        if (sum(selNonZero) == nRxnsJnc & ...
            size(unique(geneMap(selNonZero,:),'rows'),1) == nRxnsJnc)
            validJunctionMets{end+1} = junctionMets{i};
            if (verbFlag)
                fprintf('*** %s ***\n',junctionMets{i});
                for j = 1:length(thisRxns)
                    %fprintf('%s\t',thisRxns{j});
                    geneInd = find(model.rxnGeneMat(findRxnIDs(model,thisRxns{j}),:));
                    if (~isempty(geneInd))
                        thisGenes = model.genes(geneInd);
                        for k = 1:length(thisGenes)
                            fprintf('%s ',thisGenes{k});
                        end
                    end
                    fprintf('\n');
                end
            end
        end
    end
end

validJunctionMets = validJunctionMets';
for i = 1:length(validJunctionMets)
    validJunctionMets{i} = [validJunctionMets{i} '(c)'];
end

author: opencobra

repo: cobratoolbox


text:
function g = diskfun(f, varargin)
% DISKFUN returns a diskfun representing the BALLFUN F evaluated at a
% planar slice.
%
%   G = DISKFUN(F) is the slice of F in the XY plane.
%
%   G = DISKFUN(F, 'x', C) is the slice of F in the plane X = C,
%   scaled to the unit disk; G is a diskfun.
%
%   G = DISKFUN(F, 'y', C) is the slice of F in the plane Y = C,
%   scaled to the unit disk; G is a diskfun.
%
%   G = DISKFUN(F, 'z', C) is the slice of F in the plane Z = C,
%   scaled to the unit disk; G is a diskfun.
%
%   G = DISKFUN(F, PHI, THETA, PSI, C) rotates F using Euler angles phi, theta,
%   and psi with the ZXZ convention and then evaluates it in the plane Z = C
%   (C = 0 by default)
%
% % See also BALLFUN/SPHEREFUN.

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

if isempty( f )
   g = diskfun();
   return
end

% Parse user inputs to get the Euler angles
[phi, theta, psi, c] = parseInputs(varargin{:});

% Throw an error if c>1 or c<1
if abs(c) > 1
    error('CHEBFUN:BALLFUN:diskfun:sliceNotInBall',...
        ['The specified slice does not lie in unit ball.']);
end

% Rotate f using Euler angles phi, theta and psi
f = rotate(f, phi, theta, psi);

% Get the size
[m,n,~] = size(f);

% If n is odd, make it even
m = m + 1-mod(m,2);
n = n + mod(n,2);

% Evaluation points in [c, 1]
rho = chebpts(m)*sqrt(1-c^2);
rho = rho(ceil(m/2):end);

% Evaluation points in [-pi,pi)
lambda = pi*trigpts(n);

% Build the grid and evaluate at the plane Z = C
r  = sqrt(rho.^2+c^2);

theta = atan2(rho,c);
G = zeros(length(r),n);
for i = 1:length(r)
   G(i,:) = fevalm(f, r(i), lambda, theta(i));
end

% Return the diskfun
g = diskfun(real(G));
end

function [phi, theta, psi, c] = parseInputs(varargin)
% Parse user inputs to DISKFUN.
c = 0;
if nargin == 0
    phi = 0;
    theta = 0;
    psi = 0;
elseif ischar(varargin{1})
    phi = 0;
    if strcmp(varargin{1},'x')
        theta = -pi/2;
        psi = pi/2;
        % Evaluate at Y-Z plane
    elseif strcmp(varargin{1},'y')
        % Evaluate at X-Z plane
        theta = -pi/2;
        psi = 0;
    elseif strcmp(varargin{1},'z')
        % Evaluate at X-Y plane
        theta = 0;
        psi = 0;
    end
    if nargin >= 2
       c = varargin{2};
    end
else
    phi = varargin{1};
    theta = 0;
    psi = 0;
    if nargin >= 2
        theta = varargin{2};
    end
    if nargin >= 3
        psi = varargin{3};
    end
    if nargin >= 4
        c = varargin{4};
    end
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% LOGSPEC - plot mean log power spectra of submitted data on loglog scale
%             using PLOTDATA or PLOTTOPO formats
%
% Usage:
%        >> [spectra,freqs] = logspec(data,frames,srate);
%        >> [spectra,freqs] = logspec(data,frames,srate,'title',...
%                                    [loHz-hiHz],'chan_locs',rm_mean);
% Inputs:
%    data   = input data (chans,frames*epochs)
%    frames = data samples per epoch   {default length(data)}
%    srate  = data sampling rate in Hz {default 256 Hz};
%    'title' = plot title {default: none}
%    [loHz-hiHz] = [loHz hiHz] plotting limits
%       {default: [srate/fftlength srate/2]}
%    'chan_locs' = channel location file (ala TOPOPLOT)
%                   Else [rows cols] to plot data in a grid array
%    rm_mean = [0/1] 1 -> remove log mean spectrum from all
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 11-07-97
%
% See also: PLOTDATA, PLOTTOPO

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

% Changed PLOTDATA below to PLOTTOPO 11/12/99 -sm
% Mentioned new grid array option 12/22/99 -sm
% 01-25-02 reformated help & license, added links -ad

function [spectra,freqs] = logspec(data,frames,srate,titl,Hzlimits,chanlocs,rmmean)

if nargin < 1
  help logspec
  return
end
[rows,cols] = size(data);

icadefs  % read plotdata chan limit

if rows > MAXPLOTDATACHANS
   fprintf('logspec(): max plotdata() channels is %d.\n',MAXPLOTDATACHANS);
   return
end
if rows < 2
   fprintf('logspec(): min plotdata() channels is %d.\n',2);
   return
end
if nargin<7
  rmmean = 0; % default
end

if nargin < 5
 Hzlimits = 0;
end

if nargin < 4
  titl = ' ';
end

if nargin < 3
  srate = 256;
end

if nargin < 2,
  frames = cols;
end

epochs = fix(cols/frames);
if epochs*frames ~= cols
    fprintf('logspec() - frames does not divide data length.\n');
    return
end

fftlength = 2^floor(log(frames)/log(2));

spectra = zeros(rows,epochs*fftlength/2);
f2 = fftlength/2;
dB = 10/log(10);

if length(Hzlimits) < 2
 Hzlimits = [srate/fftlength srate/2];
end
if Hzlimits(2) <= Hzlimits(1)
   help logspec
   return
end

for e=1:epochs
  for r=1:rows
    [Pxx,freqs] = pwelch(data(r,(e-1)*frames+1:e*frames),fftlength,...
                          fftlength/4,fftlength,srate);
    spectra(r,(e-1)*f2+1:e*f2) = Pxx(2:f2+1)'; % omit DC bin
  end
end

clf
freqs = freqs(2:f2+1);
fsi = find(freqs >= Hzlimits(1) & freqs <= Hzlimits(2));
minf = freqs(fsi(1));
maxf = freqs(fsi(length(fsi)));
nfs = length(fsi);

showspec = zeros(rows,length(fsi)*epochs);
for e = 1:epochs
 showspec(:,(e-1)*nfs+1:e*nfs) = dB*log(spectra(:,(e-1)*f2+fsi));
end
% minspec = min(min(showspec));
% showspec = showspec-minspec; % make minimum 0 dB

showspec = blockave(showspec,nfs);

% meanspec = mean(showspec);
% showspec = showspec - ones(rows,1)*meanspec;

% >> plotdata(data,frames,limits,title,channames,colors,rtitle)
% diff = 0;
% MINUEND = 6;
% for r=1:rows
%  diff = diff - MINUEND;
%  showspec(r,:) = showspec(r<:)-diff;
% end
% semilogx(freqs(fsi),showspec');
% ax = axis;
% axis([minf maxf ax(3) ax(4)]);
% title(titl);

if nargin<6
% >> plotdata(data,frames,limits,title,channames,colors,rtitle,ydir)
     if rmmean
        showspec = showspec - ones(rows,1)*mean(showspec);
     end
     plotdata(showspec,nfs,[minf maxf 0 0],titl);
else
% >> plottopo(data,'chan_locs',frames,limits,title,channels,axsize,colors,ydir)
     if rmmean
        showspec = showspec - ones(rows,1)*mean(showspec);
     end
     plottopo(showspec,chanlocs,nfs,[minf maxf 0 0],titl);
end
ax = get(gcf,'children');
for a = ax
  set(a,'XScale','log')
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function GEplot(filename,Lat,Lon,style,opt11,opt12,opt21,opt22)
%
% function GEplot(filename,Lat,Lon,style,opt11,opt12,opt21,opt22)
%
% Description: creates a file in kmz format that can be opened into Google Earth.
%    GEplot uses the same syntax as the traditional plot function but
%    requires Latitude and Longitudd (WGS84) instead of x and y (note that Lat is
%    the first argument).
%    If you need to convert from UTM to Lat/Lon you may use utm2deg.m also
%    available at Matlab Central
%
% Arguments:
%    filename Example 'rafael', will become 'rafael.kmz'.  The same name
%             will appear inside Temporary Places in Google Earth as a layer.
%    dot_size Approximate size of the mark, in meters
%    Lat, Lon Vectors containing Latitudes and Longitudes.  The number of marks
%             created by this function is equal to the length of Lat/Lon vectors
%(opt)style   allows for specifying symbols and colors (line styles are not
%             supported by Google Earth currently. (see bellow)
%(opt)opt...   allows for specifying symbol size and line width (see bellow)
%
% Example:
%    GEplot('my_track',Lat,Lon,'o-r','MarkerSize',10,'LineWidth',3)
%    GEplot('my_track',Lat,Lon)
%
% Plot style parameters implemented in GEplot
%            color               symbol                line style
%            -----------------------------------------------------------
%            b     blue          .     point              -    solid
%            g     green         o     circle          (none)  no line
%            r     red           x     x-mark
%            c     cyan          +     plus
%            m     magenta       *     star
%            y     yellow        s     square
%            k     black         d     diamond
%            w     white (new)   S     filled square (new)
%                                D     filled diamond (new)
%                                O     filled circle=big dot (new)
%
% Plot properties: 'MarkerSize' 'LineWidth'
%
% Additional Notes:
% 'Hold on' and 'Hold off' were not implemented since one can generate a
%    .kmz file for each plot and open all simultaneously within GE.
% Unless you have a lot of data point it is recomended to show the symbols
%    since they are created in a separate folder so within Google Earth it
%    is very easy to show/hide the line or the symbols.
% Current kml/kmz format does not support different linestyles, just solid.
%    Nevertheless it is possible to define the opacity of the color (the
%    first FF in the color definition means no transparency).
% Within Matlab, it is possible to generate a second plot with the
%    same name, then you just need to select File/Revert within GE to update.
%
%
% Author: Rafael Palacios, Universidad Pontificia Comillas
% http://www.iit.upcomillas.es/palacios
% November 2006
% Version 1.1: Fixed an error while plotting graphs without symbols.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Argument checking
%
error(nargchk(3, 8, nargin));  %3 arguments required, 8 maximum
n1=length(Lat);
n2=length(Lon);
if (n1~=n2)
   error('Lat and Lon vectors should have the same length');
end
if (nargin==5 || nargin==7)
    error('size arguments must be "MarkerSize" or "LineWidth" strings followed by a number');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% symbol size and line width
%
markersize=7; %matlab default
linewidth=2;  %matlab default is 0.5, too thin for map overlay
if (nargin==6)
    if (strcmpi(opt11,'markersize')==1)
        markersize=opt12;
    elseif (strcmpi(opt11,'linewidth')==1)
        linewidth=opt12;
    else
        error('size arguments must be "MarkerSize" or "LineWidth" strings followed by a number');
    end
end
if (nargin==8)
    if (strcmpi(opt21,'markersize')==1)
        markersize=opt22;
    elseif (strcmpi(opt21,'linewidth')==1)
        linewidth=opt22;
    else
        error('size arguments must be "MarkerSize" or "LineWidth" strings followed by a number');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% symbol, line style and color
%
symbol='none';
iconfilename='none';
linestyle='-';
color='b';
colorstring='ffff0000';

if (nargin>=4)
    %linestyle
    if (strfind(style,'-'))
        linestyle='-';
    else
        linestyle='none';
    end

    %symbol
    if (strfind(style,'.')), symbol='.'; iconfilename='dot'; end
    if (strfind(style,'o')), symbol='o'; iconfilename='circle'; end
    if (strfind(style,'x')), symbol='x'; iconfilename='x'; end
    if (strfind(style,'+')), symbol='+'; iconfilename='plus'; end
    if (strfind(style,'*')), symbol='*'; iconfilename='star'; end
    if (strfind(style,'s')), symbol='s'; iconfilename='square'; end
    if (strfind(style,'d')), symbol='d'; iconfilename='diamond'; end
    if (strfind(style,'S')), symbol='S'; iconfilename='Ssquare'; end
    if (strfind(style,'D')), symbol='D'; iconfilename='Sdiamon'; end
    if (strfind(style,'O')), symbol='O'; iconfilename='dot'; end
    if (strfind(style,'0')), symbol='O'; iconfilename='dot'; end

    %color
    if (strfind(style,'b')), color='b'; colorstring='ffff0000'; end
    if (strfind(style,'g')), color='g'; colorstring='ff00ff00'; end
    if (strfind(style,'r')), color='r'; colorstring='ff0000ff'; end
    if (strfind(style,'c')), color='c'; colorstring='ffffff00'; end
    if (strfind(style,'m')), color='m'; colorstring='ffff00ff'; end
    if (strfind(style,'y')), color='y'; colorstring='ff00ffff'; end
    if (strfind(style,'k')), color='k'; colorstring='ff000000'; end
    if (strfind(style,'w')), color='w'; colorstring='ffffffff'; end
end

iconfilename=strcat('GEimages/',iconfilename,'_',color,'.png');
if (symbol=='.')
    markersize=markersize/5;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creating kml file
%
fp=fopen(strcat(filename,'.kml'),'w');
if (fp==-1)
    message=sprint('Unable to open file %s.kml',filename);
    error(message);
end
fprintf(fp,'<?xml version="1.0" encoding="UTF-8"?>\n');
fprintf(fp,'<kml xmlns="http://earth.google.com/kml/2.1">\n');
fprintf(fp,'<Document>\n');
fprintf(fp,'<name>%s</name>\n',strcat(filename,'.kml'));
fprintf(fp,'<description>Graph generated in Matlab (using GEplot.m by Rafael Palacios)</description>\n');
%
%Symbol styles definition
fprintf(fp,'<Style id="mystyle">\n');
fprintf(fp,'   <IconStyle>\n');
fprintf(fp,'      <scale>%.2f</scale>\n',markersize/14); %scale adjusted for .png image sizes
fprintf(fp,'      <Icon><href>%s</href></Icon>\n',iconfilename);
fprintf(fp,'   </IconStyle>\n');
fprintf(fp,'   <LineStyle>\n');
fprintf(fp,'      <color>%s</color>\n',colorstring);
fprintf(fp,'      <width>%d</width>\n',linewidth);
fprintf(fp,'   </LineStyle>\n');
fprintf(fp,'</Style>\n');
fprintf(fp,'\n');

if (linestyle=='-')
    fprintf(fp,'    <Placemark>\n');
    fprintf(fp,'      <description><![CDATA[Line created with Matlab GEplot.m]]></description>\n');
    fprintf(fp,'      <name>Line</name>\n');
    fprintf(fp,'      <visibility>1</visibility>\n');
    fprintf(fp,'      <open>1</open>\n');
    fprintf(fp,'      <styleUrl>mystyle</styleUrl>\n');
    fprintf(fp,'      <LineString>\n');
    fprintf(fp,'        <extrude>0</extrude>\n');
    fprintf(fp,'        <tessellate>0</tessellate>\n');
    fprintf(fp,'        <altitudeMode>clampToGround</altitudeMode>\n');
    fprintf(fp,'        <coordinates>\n');
    for k=1:n1
      fprintf(fp,'%11.6f, %11.6f, 0\n',Lon(k),Lat(k));
    end
    fprintf(fp,'        </coordinates>\n');
    fprintf(fp,'      </LineString>\n');
    fprintf(fp,'    </Placemark>\n');
end

if (strcmp(symbol,'none')==0)
fprintf(fp,'    <Folder>\n');
fprintf(fp,'      <name>Data points</name>\n');
for k=1:n1
    fprintf(fp,'      <Placemark>\n');
    fprintf(fp,'         <description><![CDATA[Point created with Matlab GEplot.m]]></description>\n');
%    fprintf(fp,'         <name>Point %d</name>\n',k);  %you may add point labels here
    fprintf(fp,'         <visibility>1</visibility>\n');
    fprintf(fp,'         <open>1</open>\n');
    fprintf(fp,'         <styleUrl>#mystyle</styleUrl>\n');
    fprintf(fp,'         <Point>\n');
    fprintf(fp,'           <coordinates>\n');
    fprintf(fp,'%11.6f, %11.6f, 0\n',Lon(k),Lat(k));
    fprintf(fp,'           </coordinates>\n');
    fprintf(fp,'         </Point>\n');
    fprintf(fp,'      </Placemark>\n');
end
fprintf(fp,'    </Folder>\n');
end

fprintf(fp,'</Document>\n');
fprintf(fp,'</kml>\n');

fclose(fp);

if (strcmp(symbol,'none')==1)
   zip(filename,{strcat(filename,'.kml')});
else
   zip(filename,{strcat(filename,'.kml'), iconfilename});
end
movefile(strcat(filename,'.zip'),strcat(filename,'.kmz'));
delete(strcat(filename,'.kml'));


text:
% APM Load Application
%
% app = apm_app(server,name)
%
% Function apm_app loads a model file (apm) and optionally
%   a data file (csv) to the APM server with the
%   following arguments:
%
%      app = application name
%   server = address of server
%     name = loads name.apm and name.csv
%
function [app] = apm_app(server,name)

% model name with .apm extension
app_model = [name '.apm'];

% data file with .csv extension (optional)
app_data = [name '.csv'];

% application name, use random number to avoid ip conflicts
app = [name '_' int2str(rand()*10000)];
app = lower(deblank(app));

% clear previous application
apm(server,app,'clear all');

% check that model file exists (required)
if (~exist(app_model,'file')),
    disp(['Error: file ' app_model ' does not exist']);
    app = [];
    return
else
    % load model file
    apm_load(server,app,app_model);
end

% check if data file exists (optional)
if (exist(app_data,'file')),
    % load data file
    csv_load(server,app,app_data);
end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function [c, ceq] = unitdisk(x)
% Unit disk. The vector x should be a row vector. This is to test the
% non-linear constraint functionality of pso.
%
% Non-linear constraints such that
% c(x) <= 0
% ceq(x) = 0
% Within a tolerance of options.TolCon

c = x*x' - 1 ;
ceq = [] ;


text:
function checkDir(path2outputs)

if ~exist(path2outputs, 'dir')
    mkdir(path2outputs);
end

author: shuochsu

repo: DeepVideoDeblurring


text:
function focusWindowSizePopup_Callback(hObject, eventdata, handles)
dropdown_items = cellstr(get(hObject,'String'));
focus_seconds = regexp(dropdown_items{get(hObject,'Value')},'([\d*.])*','match');
focus_seconds = str2double(focus_seconds{1});
handles.data.settings.focus_window_size = focus_seconds;
handles.data.saveSettings();

if ~isempty(handles.data.audiodata)
    update_fig(hObject, eventdata, handles);
else
    guidata(hObject, handles);
end

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
% This script generates paths of a subordinated Brownian motion with CIR-induced subordinator

% see A. Meucci (2009)
% "Review of Discrete and Continuous Processes in Finance - Theory and Applications"
% available at ssrn.com

% Code by A. Meucci, April 2009
% Most recent version available at www.symmys.com > Teaching > MATLAB

clc; clear; close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m=.1;
s=.4;

Kappa=.6;  % 2*Kappa*T_dot>Lambda^2;
T_dot=1;
Lambda=1;
T=252*10;
dt=1/252;
J=2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dB=sqrt(dt)*randn(J,T);
yt=1;
y=[];
d_Xs=[];
d_taus=[];
for t=1:T
    dy=Kappa*(T_dot-yt)*dt+ Lambda*sqrt(yt).*dB(:,t);
    yt=max(yt+dy,10^(-10));

    d_tau=yt*dt;
    dX=normrnd(m*d_tau,s*sqrt(d_tau));

    y=[y yt];
    d_taus=[d_taus d_tau];
    d_Xs=[d_Xs dX];

end
tau=cumsum(d_taus,2);
X=cumsum(d_Xs,2);

figure
subplot(2,1,1)
h3=plot(dt*[1:T],X(1,:),'k');
title('CIR-subordinated process')
grid on

subplot(2,1,2)
h1=plot(dt*[1:T],y(1,:));
hold on
h2=plot(dt*[1:T],tau(1,:),'r');
grid on
legend('CIR','stoch. time','location','northwest')


text:
function [trl, event] = ft_trialfun_trial(cfg)

% FT_TRIALFUN_TRIAL creates a trial definition that corresponds to the
% events that are returned by FT_READ_EVENT with type='trial'
%
% You can use this function as follows
%   cfg           = [];
%   cfg.dataset   = string, containing filename or directory
%   cfg.trialfun  = 'ft_trialfun_trial';
%   cfg           = definetrial(cfg);
%   data          = preprocessing(cfg);
%
% See also FT_DEFINETRIAL, FT_PREPROCESSING

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

if isfield(cfg, 'event')
  % for BCI applications events should be specified in the cfg
  % to prevent reading the same events many times
  event = cfg.event;
else
  event = ft_read_event(cfg.dataset);
end

sel = find(strcmp({event.type}, 'trial'));
trl = zeros(length(sel),3);

smp = [event.sample];

for i=1:length(sel)
  % determine the begin, end and offset for each trial and add it to the Nx3 matrix
  begsample = event(sel(i)).sample;
  endsample = begsample + event(sel(i)).duration - 1;
  offset    = event(sel(i)).offset;

  tmpsel    = find(smp==smp(sel(i))-offset);
  tmpval    = zeros(1,0);
  for k=1:length(tmpsel)
    if ~isempty(event(tmpsel(k)).value)
      % this assumes per trial that the triggers occur in a standardised
      % order, otherwise the entries per column will have different
      % meanings
      tmpval = [tmpval event(tmpsel(k)).value];
    end
  end

  trl(i,1:(length(tmpval)+3))  = [begsample endsample offset tmpval];
end

author: spm


text:
function sF = dtheta(sF)
% first derivative in direction theta

s = size(sF);
sF = reshape(sF, []);

sF.fhat(1) = 0; % exclude some special cases
fhat = zeros((sF.bandwidth+2)^2, length(sF));
for m = 0:sF.bandwidth+1
  if 0 <= m-1
    fhat(m^2+1:(m+1)^2, :) = (m-1)*sqrt((m^2-(-m:m)'.^2)/((2*m-1)*(2*m+1))).*[zeros(1, length(sF)); sF.fhat((m-1)^2+1:m^2, :); zeros(1, length(sF))];
  end
  if m+1 <= sF.bandwidth
    fhat(m^2+1:(m+1)^2, :) = fhat(m^2+1:(m+1)^2, :)-(m+2)*sqrt(((m+1)^2-(-m:m)'.^2)/((2*m+1)*(2*m+3))).*sF.fhat((m+1)^2+2:(m+2)^2-1, :);
  end
end

sF = S2FunHarmonic(fhat);
f = @(v) sF.eval(v)./max(sin(v.theta), eps);
sF = S2FunHarmonic.quadrature(f, 'bandwidth', sF.bandwidth);
sF = reshape(sF, s);

end

author: mtex-toolbox

repo: mtex