text:
function y = sum_square_abs( x, dim )

%SUM_SQUARE_ABS   Sum of the squares of absolute values.
%   For real arrays, SUM_SQUARE_ABS(X) computes the same result as
%   SUM_SQUARE(X). For complex arrays, SUM_SQUARE(X) first computes the
%   magnitudes of the elements of X, so it compute SUM_SQUARE_ABS(X).
%
%   Similarly, SUM_SQUARE_ABS(X,DIM) implements SUM_SQUARE(ABS(X),DIM).
%
%   Disciplined convex programming information:
%       SUM_SQUARE_ABS(X,...) is convex and nonmonotonic in X. Thus, when
%       used in CVX expressions, X must be affine. DIM must be constant.

narginchk(1,2);
y = conj( x ) .* x;
if nargin == 2,
    y = sum( y, dim );
else
    y = sum( y );
end

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

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function rs = surrogate2(s)

%tstoolbox/@signal/surrogate2
%   Syntax:
%     * rs = surrogate2(s)
%
%   create surrogate data for a scalar time series
%   see : James Theiler et al.'Using Surrogate Data to Detect Nonlinearity
%   in Time Series', APPENDIX : ALGORITHM II
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(1,1);

c = surrogate2(s.core); 		% call real working routine for parent core object
rs = signal(c, s);				% special constructor calling syntax for working routines

rs = addhistory(rs,  ['Surrogated with Theiler Algorithm II'] );
rs = addcommandlines(rs, 's = surrogate2(s');

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function varargout=savebnm(fname,varargin);
% a compact binary format for saving a matrix. (using float)
%
%
% savebnm(fname,M1[,M2,M3,...]);
%
% Can save multidimensional matrices containing: Nans,Imaginary numbers and Inf
%
% see loadbnm
%
% Aslak Grinsted 2004 april
%

% -------------------------------------------------------------------------
%   Copyright (C) 2002-2004, Aslak Grinsted
%   This software may be used, copied, or redistributed as long as it is not
%   sold and this copyright notice is reproduced on each copy made.  This
%   routine is provided as is without any express or implied warranties
%   whatsoever.

% TODO: save inputnames also

if (length(strfind(fname,'.'))==0)
    fname=[fname '.bnm'];
end

for ii=1:length(varargin) %validate arguments
    if ~isnumeric(varargin{ii})
        error('BNM format only supports numeric data!');
    end
end

[fid,message]=fopen(fname,'w','ieee-le');
if (fid<0)
    message=strrep(message,'Sorry. No help in figuring out the problem . . .','');
    error(['Could not open the specified file for writing. ' message ' file=' fname])
end

fwrite(fid,'BNM Matrix File');
version=2;

fwrite(fid,[version length(varargin)],'int16');

for ii=1:length(varargin)
    %[ii ftell(fid)]
    %M=varargin{ii};
    sz=size(varargin{ii});
    fwrite(fid,[length(sz) isreal(varargin{ii})],'int16');
    fwrite(fid,sz,'uint32');
    fwrite(fid,varargin{ii},'float');
    if ~isreal(varargin{ii})
        fwrite(fid,imag(varargin{ii}),'float');
    end
end
fclose(fid);

function savestring(fid,s);
if ~ischar(s)|(size(s,1)~=1)
    error('not a string')
end
fwrite(fid,length(s),'uint32');
fwrite(fid,s);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% The COBRAToolbox: testmultiProductionEnvelope.m
%
% Purpose:
%     - test the multiProductionEnvelope function
%
% Authors:
%     - Jacek Wachowiak
global CBTDIR
% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testMultiProductionEnvelope'));
cd(fileDir);

% test variables
model = getDistributedModel('ecoli_core_model.mat');
model.lb(36) = 0; % setting the model to anaerobic conditions
model.ub(36) = 0; % setting the model to anaerobic conditions
biomassRxn = model.rxns(13);
target = model.rxns(13); % objective = biomass reaction
deletions = model.rxns(21);
deletions2 = model.genes(21);
targetreactionsBasic = {'EX_ac(e)','EX_for(e)','EX_ethoh(e)'};

%reference data
refData_x=(0:(0.2117/19):0.2117)';
refData_lowerbound = zeros(20, 1); % EX_ac(e)
refData_lowerbound(18) = 0.0736;
refData_lowerbound(19) = 3.9862;
refData_lowerbound(20) = 8.5036;
 refData_upperBound = (10:(-1.4964/19):8.5036)';

% function outputs
% each column of targetValues corresponds to one reaction results, while biomassValues is the x axis
% for this test: 1	EX_ac(e), 2	EX_acald(e), 3	EX_akg(e), 4	EX_etoh(e), 5	EX_for(e), after that the results do not correspond to the graph,
% additionally there are always 2 curves, the data corresponds to the lower one
[biomassValues, targetLowerBounds, targetUpperBounds, targetReactions] = multiProductionEnvelope(model);
[biomassValues2, targetLowerBounds2, targetUpperBounds2, targetReactions2] = multiProductionEnvelope(model, deletions, biomassRxn);
%gene not reaction removal
[biomassValues3, targetLowerBounds3, targetUpperBounds3,  targetReactions3] = multiProductionEnvelope(model, deletions2, biomassRxn, 1, 20);

% tests - not all results are possible to be tested, the suitable were chosen
% x axis comparison
assert(isequal((abs(refData_x-biomassValues) < 1e-4), true(20, 1)));
assert(isequal((abs(refData_lowerbound-targetLowerBounds(:,ismember(targetReactions,'EX_ac(e)'))) < 1e-4), true(20, 1)));
assert(isequal((abs(refData_upperBound-targetUpperBounds(:, ismember(targetReactions,'EX_ac(e)'))) < 1e-4), true(20, 1)));

close all hidden force

% change to old directory
cd(currentDir);

author: opencobra

repo: cobratoolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright 2014 National Renewable Energy Laboratory and National
% Technology & Engineering Solutions of Sandia, LLC (NTESS).
% Under the terms of Contract DE-NA0003525 with NTESS,
% the U.S. Government retains certain rights in this software.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
% http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

classdef mooringClass<handle
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % The ``mooringClass`` creates a ``mooring`` object saved to the MATLAB
    % workspace. The ``mooringClass`` includes properties and methods used
    % to define cable connections relative to other bodies.
    % It is suggested that the ``mooringClass`` be used for connections between
    % bodies and the global reference frame.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % This class contains mooring parameters and settings
    properties (SetAccess = 'public', GetAccess = 'public')%input file
        initial                 = struct(...                               % (`structure`) Defines the initial displacement of the mooring.
            'displacement',         [0 0 0], ...                           %
            'axis',                 [0 1 0], ...                           %
            'angle',                0)                                     % (`structure`) Defines the initial displacement of the mooring. ``displacement`` (`3x1 float vector`) is defined as the initial displacement of the pto [m] in the following format [x y z], Default = [``0 0 0``].
        location                = [0 0 0]                                  % (`float 1 x 3`) Mooring Reference location. Default = ``[0 0 0]``
        matrix                  = struct(...                               % (`structure`) Defines the mooring parameters.
            'damping',              zeros(6,6), ...                        %
            'stiffness',            zeros(6,6), ...                        %
            'preTension',           [0 0 0 0 0 0])                         % (`structure`) Defines the mooring parameters. ``damping`` (`6x6 float matrix`) Matrix of damping coefficients, Default = ``zeros(6)``. ``stiffness`` (`6x6 float matrix`) Matrix of stiffness coefficients, Default = ``zeros(6)``. ``preTension`` (`6x6 float matrix`) Array of pretension force in each dof, Default = ``[0 0 0 0 0 0]``.
        moorDyn                 = 0                                        % (`integer`) Flag to indicate a MoorDyn block, 0 or 1. Default = ``0``
        moorDynLines            = 0                                        % (`integer`) Number of lines in MoorDyn. Default = ``0``
        moorDynNodes            = []                                       % (`integer`) number of nodes for each line. Default = ``'NOT DEFINED'``
        name                    = 'NOT DEFINED'                            % (`string`) Name of the mooring. Default = ``'NOT DEFINED'``
    end

    properties (SetAccess = 'private', GetAccess = 'public') %internal
        orientation             = []                                       % (`float 1 x 6`) Initial 6DOF location. Default = ``[0 0 0 0 0 0]``
        number                  = []                                       % (`integer`) Mooring number. Default = ``'NOT DEFINED'``
    end

    methods (Access = 'public')
        function obj = mooringClass(name)
            % This method initializes the mooringClass object
            if exist('name','var')
                obj.name = name;
            else
                error('The mooring class number(s) in the wecSimInputFile must be specified in ascending order starting from 1. The mooringClass() function should be called first to initialize each mooring line with a name.')
            end
        end

        function setInitDisp(obj, relCoord, axisAngleList, addLinDisp)
            % Function to set a mooring's initial displacement
            %
            % This function assumes that all rotations are about the same relative coordinate.
            % If not, the user should input a relative coordinate of 0,0,0 and
            % use the additional linear displacement parameter to set the cg or orientation
            % correctly.
            %
            % Parameters
            % ------------
            %    relCoord : [1 3] float vector
            %        Distance from x_rot to the body center of gravity or the constraint
            %        or pto location as defined by: relCoord = cg - x_rot. [m]
            %
            %    axisAngleList : [nAngle 4] float vector
            %        List of axes and angles of the rotations with the
            %        format: [n_x n_y n_z angle] (angle in rad)
            %        Rotations applied consecutively in order of dimension 1
            %
            %    addLinDisp : [1 3] float vector
            %        Initial linear displacement (in addition to the
            %        displacement caused by rotation) [m]
            %

            % initialize quantities before for loop
            axisList = axisAngleList(:,1:3);
            angleList = axisAngleList(:,4);
            nAngle = size(axisList,1);
            rotMat = eye(3);
            % Loop through all axes and angles.
            for i=1:nAngle
                rotMat = axisAngle2RotMat(axisList(i,:),angleList(i))*rotMat;
            end
            % calculate net axis-angle rotation
            [netAxis, netAngle] = rotMat2AxisAngle(rotMat);
            % calculate net displacement due to rotation
            rotatedRelCoord = relCoord*(rotMat');
            linDisp = rotatedRelCoord - relCoord;
            % apply rotation and displacement to object
            obj.initial.displacement = linDisp + addLinDisp;
            obj.initial.axis = netAxis;
            obj.initial.angle = netAngle;
        end

        function listInfo(obj)
            % Method to list mooring info
            fprintf('\n\t***** Mooring Name: %s *****\n',obj.name)
        end

        function obj = setLoc(obj)
            % This method sets mooring location
            obj.orientation = [obj.location + obj.initial.displacement 0 0 0];
        end

        function setNumber(obj,number)
            % Method to set the private number property
            obj.number = number;
        end
    end
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function K = nddisimXndsimKernCompute(disimKern, simKern, t1, t2)

% NDDISIMXNDSIMKERNCOMPUTE Compute a cross kernel between DISIM and SIM kernels with no decay in the SIM part.
% FORMAT
% DESC computes cross kernel terms between DISIM and SIM kernels for
% the multiple output kernel.
% ARG disimKern : the kernel structure associated with the DISIM
% kernel.
% ARG simKern : the kernel structure associated with the SIM
% kernel.
% ARG t : inputs for which kernel is to be computed.
% RETURN k : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between DISIM and SIM kernels for
% the multiple output kernel.
% ARG disimKern : the kernel structure associated with the DISIM
% kernel.
% ARG simKern : the kernel structure associated with the SIM
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% RETURN k : block of values from kernel matrix.
%
% SEEALSO : multiKernParamInit, multiKernCompute, disimKernParamInit, simKernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% COPYRIGHT : Antti Honkela, 2007-2009
%
% COPYRIGHT : Jaakko Peltonen, 2011

% KERN

%fprintf(1,'nddisimXndsimKernCompute step1\n');
%t1
%t2

if nargin < 4
  t2 = t1;
end
if size(t1, 2) > 1 | size(t2, 2) > 1
  error('Input can only have one column');
end

if (isempty(t1)) || (isempty(t2)),
  K = zeros(size(t1,1),size(t2,1));
  return;
end;

if isfield(disimKern,'delay'),
  t1=t1-disimKern.delay;

  % crude way to handle times below zero, just truncate to zero
  % since the cross-kernel value at t1=0 is zero, which is the same
  % as cross-kernel values for negative t.
  I=find(t1<0);
  t1(I)=0;
end;

if disimKern.inverseWidth ~= simKern.inverseWidth
  error('Kernels cannot be cross combined if they have different inverse widths.')
end
%if disimKern.di_decay ~= simKern.decay
%  error('Kernels cannot be cross combined if they have different driving input decays.');
%end
if disimKern.di_variance ~= simKern.variance
  error('Kernels cannot be cross combined if they have different driving input variances.');
end

dim1 = size(t1, 1);
dim2 = size(t2, 1);
%size(t1)
%t1
%size(t2)
%t2
t1Mat = t1(:, ones(1, dim2));
t2Mat = t2(:, ones(1, dim1))';
diffT = (t1Mat - t2Mat);

l = sqrt(2/disimKern.inverseWidth);
D_i = disimKern.decay;
%delta = disimKern.di_decay;

sigma = sqrt(2/simKern.inverseWidth);
if isfield(simKern, 'isNegativeS') && (simKern.isNegativeS == true)
    variancemultiplier = (simKern.sensitivity*simKern.sensitivity);
else
    variancemultiplier = simKern.variance;
end
variancemultiplier=variancemultiplier*sqrt(disimKern.variance);

k1a=-sqrt(pi)*variancemultiplier*sigma/(2*D_i*D_i)*...
    (exp((D_i*sigma/2)^2-D_i*t1Mat+D_i*t2Mat+lnDiffErfs(D_i*sigma/2+t2Mat/sigma,D_i*sigma/2+t2Mat/sigma-t1Mat/sigma))...
     +exp((D_i*sigma/2)^2-D_i*t1Mat+lnDiffErfs(D_i*sigma/2-t1Mat/sigma,D_i*sigma/2)));
k2a=-variancemultiplier*sigma*sigma/(2*D_i)*...
    (exp(-(diffT/sigma).^2)-exp(-(t1Mat/sigma).^2)+1-exp(-(t2Mat/sigma).^2));
k3a=-sqrt(pi)*variancemultiplier*sigma/(2*D_i)*...
    ((diffT-1/D_i).*erf(diffT/sigma)-(t1Mat-1/D_i).*erf(t1Mat/sigma)-(t2Mat+exp(-D_i*t1Mat)/D_i).*erf(t2Mat/sigma));

%
%k1b
%k2b
%pause
K=real(k1a+k2a+k3a);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
clear;
t0=0.01;
tf=0.05;
b0=[1.2 2.2 1.8];
[t,b]=ode45('dfun2',[t0, tf],b0);
plot(b(:,1),b(:,2));
hold on

xlabel('x1');
ylabel('x2');

author: sfvsfv

repo: Mathematical-modeling


text:
function DBListFigures(query)

%DBListFigures - List figures matching certain criteria.
%
%  USAGE
%
%    DBListFigures(query)
%
%    query          optional database query (WHERE clause)
%
%  EXAMPLE
%
%    DBListFigures('name="Spectrogram"');
%
%  SEE
%
%    See also DBGetVariables, DBGetFigures.
%

% Copyright (C) 2013 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Make sure MyM is installed and functional
CheckMyM;

if nargin < 1,
	query = '';
end

DBDisplay(DBGetFigures(query,'output','keys'));


text:
function cvx_optval = prod_inv( x, dim, p )

%DET_INV   Internal cvx version.

narginchk(1,3);
if ~isreal( x ),
    error( 'First argument must be real.' );
end
sx = size( x );
if nargin < 2 || isempty( dim ),
    dim = cvx_default_dimension( sx );
elseif ~cvx_check_dimension( dim ),
    error( 'Second argument must be a positive integer.' );
end
sx( end + 1 : dim ) = 1;
if nargin < 2,
    p = 1;
elseif ~isnumeric( p ) || ~isreal( p ) || numel( p ) ~=  1 || p <= 0,
    error( 'Third argument must be a positive scalar.' );
end

if cvx_isconstant( x ),

    cvx_optval = cvx( prod_inv( cvx_constant( x ), dim, p ) );

elseif sx( dim ) == 1,

    cvx_optval = inv_pos( x );

else

    sy = sx;
    sy( dim ) = 1;
    y = [];
    cvx_begin
        epigraph variable y( sy )
        geo_mean( cat( dim, x, y ), dim, [ ones(n,1) ; p ] ) >= 1; %#ok
    cvx_end

end

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


text:
function xDup=dupEls(x,numTimes)
%%DUPELS Return a vector with the elements of the vector x duplicated
%        numTimes.
%
%INPUTS: x A vectors whose elements are to be duplicated.
% numTimes The number of times the elements in x are to be duplicated.
%
%OUTPUTS: xDup A vector the same orientation (row or column vector) as x
%              where the elements of x have been duplicated numTimes. For
%              example, if x=[1;2;3], then dupEls(x,3) returns
%              [1;1;1;2;2;2;3;3;3].
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=length(x);
xDup=zeros(xDim*numTimes,1);

%If the input was a row vector.
if(size(x,1)==1)
    xDup=xDup';
end

sel=(0:(xDim-1))*numTimes;
for curDup=1:numTimes
    xDup(sel+curDup)=x;
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function determ = bab_determinant ( n, alpha, beta )

%*****************************************************************************80
%
%% BAB_DETERMINANT computes the determinant of the BAB matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real ALPHA, BETA, parameters that define the matrix.
%
%    Output, real DETERM, the determinant.
%
  determ_nm1 = alpha;

  if ( n == 1 )
    determ = determ_nm1;
    return
  end

  determ_nm2 = determ_nm1;
  determ_nm1 = alpha * alpha - beta * beta;

  if ( n == 2 )
    determ = determ_nm1;
    return
  end

  for i = n - 2 : -1 : 1

    determ = alpha * determ_nm1 - beta * beta * determ_nm2;

    determ_nm2 = determ_nm1;
    determ_nm1 = determ;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ index ] = gsp_good_graph_index( G, X, param )
%GSP_GOOD_GRAPH_INDEX Index testing how well a given graph G, matches some data X
%   Usage: gsp_good_graph(G, X);
%
%   Input parameters:
%       G          : the graph
%       X          : a data matrix
%       param      : structure of optional parameters
%   Output parameters:
%       index      : the computed index
%
%   A wrapper function with which one may test how well a given graph G,
%   matches some data X.
%
%   Example:::
%       G = gsp_2dgrid(16);
%       X = pinv(full(G.L)) * randn(G.N, G.N);
%       param.verbose = 1;
%       param.index = 'tcer';
%       index = gsp_good_graph_index(G, X, param)
%       param.index = 'stationarity';
%       index = gsp_good_graph_index(G, X, param)
%
%   Optional paramaters
%   -------------------
%
%   * *param.index*: 'tcer' or 'stationarity' (default 'tcer').
%
%   See also: gsp_learn_tcer, gsp_stationarity_ratio

% Author  : Andreas Loukas
% Date    : 15 Nov 2016

% Handle input
if nargin < 3, param = struct(); end
if not(isfield(param, 'index'));   param.index = 'tcer'; end;

switch param.index,

    case 'tcer',
        index = gsp_learn_tcer(G, X, param);

    case 'stationarity'
        index = gsp_stationarity_ratio(G, X*X', param);

    otherwise,
        error('uknown index.');

end
end

repo: gspbox


text:
% A Thin Plate Subjected to Uniform Traction
% T3 Implementation
% 2 elements
% clear memory
clear all;
clc;
close all;
% materials
E  = 30e6;     poisson = 0.30;  thickness = 1;

% matrix D
D=E/(1-poisson^2)*[1 poisson 0;poisson 1 0;0 0 (1-poisson)/2];

% trivial preprocessing
% numberElements: number of elements
numberElements=2;
% numberNodes: number of nodes
numberNodes=4;
% coordinates and connectivities
elementNodes=[1 3 2; 1 4 3];
nodeCoordinates=[0, 0; 0, 10; 20, 10; 20, 0];
drawingMesh(nodeCoordinates,elementNodes,'T3','k-o');

% GDof: global number of degrees of freedom
GDof=2*numberNodes;

% boundary conditions
prescribedDof=[1 2 3 4]';
% force vector
force=zeros(GDof,1);
force(5)=5000; force(7) =5000;

% calculation of the system stiffness matrix
stiffness=formStiffness2D(GDof,numberElements,...
    elementNodes,numberNodes,nodeCoordinates,D,thickness);

% solution
displacements=solution(GDof,prescribedDof,stiffness,force);

% output displacements
outputDisplacements(displacements, numberNodes, GDof);

outputStress(displacements,numberElements,...
    elementNodes,nodeCoordinates,D)

repo: Aero-matlab


text:
function x = complex_lorentz( sx, dim )

%COMPLEX_LORENTZ   Complex second-order cone.
%   COMPLEX_LORENTZ(N), where N is a positive integer, creates a column
%   variable of length N and a scalar variable, and constrains them
%   to lie in a second-order cone. That is, given the declaration
%       variable x(n) complex
%       variable y
%   the constraint
%       {x,y} == complex_lorentz(n)
%   is equivalent to
%       norm(x,2) <= y
%   The inequality form is more natural, and preferred in most cases. But
%   in fact, the COMPLEX_LORENTZ set form is used by CVX itself to convert
%   complex NORM()-based constraints to solvable form.
%
%   COMPLEX_LORENTZ(SX,DIM), where SX is a valid size vector and DIM is a
%   positive integer, creates an array variable of size SX and an array
%   variable of size SY (see below) and applies the second-order cone
%   constraint along dimension DIM. That is, given the declarations
%       sy = sx; sy(min(dim,length(sx)+1))=1;
%       variable x(sx) complex
%       variable y(sy)
%   the constraint
%       {x,y} == complex_lorentz(sx,dim)
%   is equivalent to
%       norms(x,2,dim) <= y
%   Again, the inequality form is preferred, but CVX uses the set form
%   internally. DIM is optional; if it is omitted, the first non-singleton
%   dimension is used.
%
%   LORENTZ(SX,DIM,CPLX) creates real second-order cones if CPLX is FALSE,
%   and complex second-order cones if CPLX is TRUE. The latter case is
%   equivalent to COMPLEX_LORENTZ(SX,DIM).
%
%   Disciplined convex programming information:
%       LORENTZ is a cvx set specification. See the user guide for
%       details on how to use sets.

if nargin == 1,
    x = lorentz( sx, [], true );
else
    x = lorentz( sx, dim, true );
end

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

author: yu-jiang

repo: radpbook

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


text:
function Y = uminus(Y)
%UMINUS (overloaded)

if isa(Y,'blkvar')
    Y = sdpvar(Y);
end

switch Y.typeflag
    case {0,1,2,3,4}
        Y.basis = -Y.basis;
    case 5
        error('Cone object cannot be negated');
    otherwise
end
% Reset info about conic terms
Y.conicinfo = [0 0];
Y.extra.opname='';

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [taskReport, essentialGenes, taskStructure]=checkTasksGenes(model,inputFile,printOutput,printOnlyFailed,getEssential,taskStructure)
% checkTasksGenes
%   Performs a set of simulations as defined in a task file. This function
%   is identical to "checkTasks", except it allows for the determination of
%   essential genes, rather than essential reactions.
%
%   model           a model structure
%   inputFile       a task list in Excel format. See the function
%                   parseTaskList for details (opt if taskStructure is
%                   supplied)
%   printOutput     true if the results of the test should be displayed
%                   (opt, default true)
%   printOnlyFailed true if only tasks that failed should be displayed
%                   (opt, default false)
%   getEssential    true if the essential genes should be calculated for
%                   all the tasks. (opt, default false)
%   taskStructure   structure with the tasks, as from parseTaskList. If
%                   this is supplied then inputFile is ignored (opt)
%
%   taskReport          structure with the results
%       id              cell array with the id of the task
%       description     cell array with the description of the task
%       ok              boolean array with true if the task was successful
%   essentialGenes      MxN matrix with the essential genes (M) for each
%                       task (N). An element is true if the corresponding
%                       gene is essential in the corresponding task.
%                       Failed tasks and SHOULD FAIL tasks are ignored.
%                       If getEssential is false, then essentialRxns will
%                       be returned as: false(nGenes,nTasks).
%   taskStructure       structure with the tasks, as from parseTaskList
%
%   This function is used for defining a set of tasks for a model to
%   perform. The tasks are defined by defining constraints on the model,
%   and if the problem is feasible, then the task is considered successful.
%   In general, each row can contain one constraint on uptakes, one
%   constraint on outputs, one new equation, and one change of reaction
%   bounds. If more bounds are needed to define the task, then several rows
%   can be used for each task.
%
%   Usage: [taskReport, essentialGenes, taskStructure]=checkTasksGenes(model,inputFile,...
%           printOutput,printOnlyFailed,getEssential,taskStructure)
%

if nargin<3 || isempty(printOutput)
    printOutput=true;
end
if nargin<4 || isempty(printOnlyFailed)
    printOnlyFailed=false;
end
if nargin<5 || isempty(getEssential)
    getEssential=false;
end

%Prepare the input model
model.b=zeros(numel(model.mets),2);

modelMets=upper(strcat(model.metNames,'[',model.comps(model.metComps),']'));
if ~isfield(model,'unconstrained')
    EM='Exchange metabolites should normally not be removed from the model when using checkTasks. Inputs and outputs are defined in the task file instead. Use importModel(file,false) to import a model with exchange metabolites remaining';
    dispEM(EM,false);
end

%Parse the task file
if nargin<6
    taskStructure=parseTaskList(inputFile);
end

essentialGenes=false(numel(model.genes),numel(taskStructure));

tModel=model;
taskReport=[];
for i=1:numel(taskStructure)
    taskReport.id{i,1}=taskStructure(i).id;
    taskReport.description{i,1}=taskStructure(i).description;
    %Set the inputs
    if ~isempty(taskStructure(i).inputs)
        [I, J]=ismember(upper(taskStructure(i).inputs),modelMets);
        J=J(I); %Only keep the ones with matches
        K=ismember(upper(taskStructure(i).inputs),'ALLMETS');
        L=~cellfun('isempty',strfind(upper(taskStructure(i).inputs),'ALLMETSIN'));
        %Check that all metabolites are either real metabolites or
        %ALLMETS/ALLMETSIN
        if ~all(I|K|L)
            fprintf(['ERROR: Could not find all inputs in "[' taskStructure(i).id '] ' taskStructure(i).description '"\n']);
            taskReport.ok(i,1)=false;
            tModel=model;
            continue;
        end
        if numel(J)~=numel(unique(J))
            EM=['The constraints on some input(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time'];
            dispEM(EM);
        end
        %If all metabolites should be added
        if any(K)
            %Check if ALLMETS is the first metabolite. Otherwise print a
            %warning since it will write over any other constraints that
            %are set
            if K(1)==0
                EM=['ALLMETS is used as an input in "[' taskStructure(i).id '] ' taskStructure(i).description '" but it it not the first metabolite in the list. Constraints defined for the metabolites before it will be over-written'];
                dispEM(EM,false);
            end
            %Use the first match of ALLMETS. There should only be one, but
            %still..
            tModel.b(:,1)=taskStructure(i).UBin(find(K,1))*-1;
        end
        %If metabolites in a specific compartment should be used
        if any(L)
            L=find(L);
            for j=1:numel(L)
                %The compartment defined
                compartment=upper(taskStructure(i).inputs{L(j)}(11:end-1));
                %Check if it exists in the model
                C=find(ismember(upper(model.comps),compartment));
                if any(C)
                    %Match to metabolites
                    tModel.b(model.metComps==C,1)=taskStructure(i).UBin(L(j))*-1;
                else
                    EM=['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist'];
                    dispEM(EM);
                end
            end
        end
        %Then add the normal constraints
        if any(J)
            tModel.b(J,1)=taskStructure(i).UBin(I)*-1;
            tModel.b(J,2)=taskStructure(i).LBin(I)*-1;
        end
    end
    %Set the outputs
    if ~isempty(taskStructure(i).outputs)
        [I, J]=ismember(upper(taskStructure(i).outputs),modelMets);
        J=J(I); %Only keep the ones with matches
        K=ismember(upper(taskStructure(i).outputs),'ALLMETS');
        L=~cellfun('isempty',strfind(upper(taskStructure(i).outputs),'ALLMETSIN'));
        %Check that all metabolites are either real metabolites or
        %ALLMETS/ALLMETSIN
        if ~all(I|K|L)
            fprintf(['ERROR: Could not find all outputs in "[' taskStructure(i).id '] ' taskStructure(i).description '"\n']);
            taskReport.ok(i,1)=false;
            tModel=model;
            continue;
        end
        if numel(J)~=numel(unique(J))
            EM=['The constraints on some output(s) in "[' taskStructure(i).id '] ' taskStructure(i).description '" are defined more than one time'];
            dispEM(EM);
        end
        %If all metabolites should be added
        if any(K)
            %Check if ALLMETS is the first metabolite. Otherwise print a
            %warning since it will write over any other constraints that
            %are set
            if K(1)==0
                EM=['ALLMETS is used as an output in "[' taskStructure(i).id '] ' taskStructure(i).description '" but it it not the first metabolite in the list. Constraints defined for the metabolites before it will be over-written'];
                dispEM(EM,false);
            end
            %Use the first match of ALLMETS. There should only be one, but
            %still..
            tModel.b(:,2)=taskStructure(i).UBout(find(K,1));
        end
        %If metabolites in a specific compartment should be used
        if any(L)
            L=find(L);
            for j=1:numel(L)
                %The compartment defined
                compartment=upper(taskStructure(i).outputs{L(j)}(11:end-1));
                %Check if it exists in the model
                C=find(ismember(upper(model.comps),compartment));
                if any(C)
                    %Match to metabolites
                    tModel.b(model.metComps==C,2)=taskStructure(i).UBout(L(j));
                else
                    EM=['The compartment defined for ALLMETSIN in "[' taskStructure(i).id '] ' taskStructure(i).description '" does not exist'];
                    dispEM(EM);
                end
            end
        end
        %Then add the normal constraints
        if any(J)
            tModel.b(J,1)=taskStructure(i).LBout(I);
            tModel.b(J,2)=taskStructure(i).UBout(I);
        end
    end
    %Add new rxns
    if ~isempty(taskStructure(i).equations)
        rxn.equations=taskStructure(i).equations;
        rxn.lb=taskStructure(i).LBequ;
        rxn.ub=taskStructure(i).UBequ;
        rxn.rxns=strcat({'TEMPORARY_'},num2str((1:numel(taskStructure(i).equations))'));
        %Allow for new metabolites to be added. This is because it should
        %be possible to add, say, a whole new pathway
        tModel=addRxns(tModel,rxn,3,[],true);
    end
    %Add changed bounds
    if ~isempty(taskStructure(i).changed)
        tModel=setParam(tModel,'lb',taskStructure(i).changed,taskStructure(i).LBrxn);
        tModel=setParam(tModel,'ub',taskStructure(i).changed,taskStructure(i).UBrxn);
    end

    %Solve and print
    sol=solveLP(tModel);
    if ~isempty(sol.x)
        if ~taskStructure(i).shouldFail
            taskReport.ok(i,1)=true;
            if printOnlyFailed==false && printOutput==true
                fprintf(['PASS: [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
            end
            %Calculate the essential reactions
            if getEssential==true
                if ~isfield(tModel,'rules')
                    % the getEssentialGenes function requires the model to
                    % have a "rules" field
                    tModel = generateRules(tModel);
                end
                [~, taskEssential]=getEssentialGenes(tModel);
                essentialGenes(taskEssential,i)=true;
            end
        else
            taskReport.ok(i,1)=false;
            if printOutput==true
                fprintf(['PASS (should fail): [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
            end
        end
    else
        if ~taskStructure(i).shouldFail
            taskReport.ok(i,1)=false;
            if printOutput==true
                fprintf(['FAIL: [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
            end
        else
            taskReport.ok(i,1)=true;
            if printOnlyFailed==false && printOutput==true
                fprintf(['FAIL (should fail): [' taskStructure(i).id '] ' taskStructure(i).description '\n']);
            end
        end
    end
    if taskStructure(i).printFluxes && ~isempty(sol.x)
        sol=solveLP(tModel,1);
        if ~isempty(sol.x)
            printFluxes(tModel,sol.x,false,10^-6,[],'%rxnID (%eqn):%flux\n');
            fprintf('\n');
        end
    end
    tModel=model;
end
end

author: SysBioChalmers

repo: Human-GEM


text:
function [params, names] = gibbsperiodicKernExtractParam(kern)

% GIBBSPERIODICKERNEXTRACTPARAM Extract parameters from the GIBBSPERIODIC kernel structure.
% FORMAT
% DESC extracts parameters from the Gibbs-kernel derived periodic
% kernel structure into a vector of parameters for optimisation.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
%
% FORMAT
% DESC extracts parameters and parameter names from the Gibbs-kernel derived periodic
% kernel structure.
% ARG kern : the kernel structure containing the parameters to be
% extracted.
% RETURN param : vector of parameters extracted from the kernel. If
% the field 'transforms' is not empty in the kernel structure, the
% parameters will be transformed before optimisation (for example
% positive only parameters could be logged before being returned).
% RETURN names : cell array of strings containing names for each
% parameter.
%
% SEEALSO gibbsperiodicKernParamInit, gibbsperiodicKernExpandParam, kernExtractParam, scg, conjgrad
%
% COPYRIGHT : Neil D. Lawrence, 2007

% KERN

params = zeros(1, kern.nParams);
if nargout < 2
  params(1:end-1) = modelExtractParam(kern.lengthScaleFunc);
  params(end) = kern.variance;
else
  names = cell(1, kern.nParams);
  [params(1:end-1), names(1:end-1)] = modelExtractParam(kern.lengthScaleFunc);
  params(end) = kern.variance;
  names{end} = 'variance';
end


text:
function h1 = imageplot(M,str, a,b,c)

% imageplot - diplay an image and a title
%
% Example of usages:
%   imageplot(M);
%   imageplot(M,title);
%   imageplot(M,title,1,2,1);   % to make subplot(1,2,1);
%
%   imageplot(M,options);
%
%   If you want to display several images:
%       imageplot({M1 M2}, {'title1', 'title2'});
%
%   Copyright (c) 2007 Gabriel Peyre

if nargin<2
    str = [];
end

options.null = 0;
if isstruct(str)
    options = str;
    str = '';
end

nbdims = nb_dims(M);

if iscell(M)
    q = length(M);
    if nargin<5
        c = 1;
    end
    if nargin<4
        a = ceil(q/4);
    end
    if nargin<3
        b = ceil(q/a);
    end
    if (c-1+q)>(a*b)
        warning('a and c parameters not large enough');
        a = ceil((c-1+q)/4);
        b = ceil((c-1+q)/a);
    end
    for i=1:q
        if iscell(str)
            str1 = str{i};
        else
            str1 = str;
        end
        h{i} = imageplot(M{i},str1, a,b,c-1+i);
    end
    global axlist;
    if not(isempty(axlist))
        linkaxes(axlist, 'xy');
    end

    if nargout>0
        if exist('h')
            h1 = h;
        else
            h1 = [];
        end
    end
    return;
end

if nargin==5
    global axlist;
    global imageplot_size;
    if c==1 || isempty(imageplot_size) || imageplot_size~=size(M,1)
        clear axlist;
        global axlist;
        axlist = [];
        imageplot_size = size(M,1);
    end
    axlist(end+1) = subplot(a,b,c);
end

if nbdims==1
    h = plot(M); axis tight;
elseif size(M,3)<=3
    % gray-scale or color image
    if size(M,3)==2
        M = cat(3,M, zeros(size(M,1),size(M,2)));
    end
    if not(isreal(M))
        if size(M,3)==1
            % turn into color matrix
            M = cat(3, real(M), imag(M), zeros(size(M,1),size(M,2)));
        else
            warning('Complex data');
            M = real(M);
        end
    end
    if size(M,3)==1
        colormap gray(256);
    else
        colormap jet(256);
    end
    h = imagesc(rescale(M)); axis image; axis off;
else

    if not(isfield(options, 'center') )
        options.center = .5; % here a value in [0,1]
    end
    if not(isfield(options, 'sigma'))
        options.sigma = .08; % control the width of the non-transparent region
    end
    a = compute_alpha_map('gaussian', options); % you can plot(a) to see the alphamap

    % volumetric image
    h = vol3d('cdata',rescale(M),'texture','2D');
    view(3);
    axis tight; % daspect([1 1 .4])
    colormap bone(256);
%    alphamap('rampup');
    % alphamap(.06 .* alphamap);
    vol3d(h);
end
if not(isempty(str))
    title(str);
end

if nargout>0
    if exist('h')
        h1 = h;
    else
        h1 = [];
    end
end

if nargin==5 && c==a*b
    linkaxes(axlist, 'xy');
end

function d = nb_dims(x)

% nb_dims - debugged version of ndims.
%
%   d = nb_dims(x);
%
%   Copyright (c) 2004 Gabriel Peyre

if isempty(x)
    d = 0;
    return;
end

d = ndims(x);

if d==2 && (size(x,1)==1 || size(x,2)==1)
    d = 1;
end

author: optimaltransport


text:
function[S] = adaptMedian_colfilt(I, w_max)
%
%

if nargin < 2
  w_max = 9;
end

[Nrows Ncols depth] = size(I);
S = uint8(zeros(Nrows, Ncols, depth));

for d = 1:depth,

% Get plane d
Ib = I(:,:,d);

% Initialization
w = 1;
p = 3;
cond = uint8(1);    % zero --> noiseless pixel

for l = w+2:2:w_max

  Imed = medfilt2(Ib, [p p]);
  Imin = colfilt(Ib, [p p], 'sliding', @(x) min(x));
  Imax = colfilt(Ib, [p p], 'sliding', @(x) max(x));

  m1 = uint8( ((Imin < Imed) .* (Imed < Imax)) );
  m2 = uint8( ((Imin < Ib) .* (Ib < Imax)) );

  if(l == w)
    S(:,:,d) =  m1.*(1-m2).*l;
  else
    S(:,:,d) = cond.*( m1.*(1-m2).*l ) + (1-cond).*S(:,:,d);
  end

  cond = cond.*(1-m1); % pixel that need further evaluation

  p = p+2;

end % _END_ FOR(l)

S(:,:,d) = cond.*(w_max+2) + (1-cond).*S(:,:,d);

end % _END_ FOR(d)

S = double(S);

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [KKTConstraints, details] = kkt(F,h,parametricVariables,ops);
%KKT Create KKT system
%
% [KKTConstraints, details] = kkt(Constraints,Objective,parameters)

if nargin == 2
    [KKTConstraints, details] = kkt(lmi(F),h);
elseif nargin==3
    [KKTConstraints, details] = kkt(lmi(F),h,parametricVariables);
else
    [KKTConstraints, details] = kkt(lmi(F),h,parametricVariables,ops);
end


text:
function  [yout, nfe] = mmid1(y, dy, xs, htot, nstep, odefun, la, kb, kk0)
% mmid.m is the algoritm for a 'modified midpoint' step in the Bulirsch-Stoer method
% for solving ode's, see W Press ea.in "Numerical Recipes" (1994) p723. This version
% is for parallel processing, see odeintp.m and mmid.m
%
% inputs : at xs(1..np) the dependent variable vector y(1..nv,1..np) and its derivative dy(1..nv,1..np)
% are input, also the total step htot(1..np) to be made  and the number of substeps used (nstep),
% a client function for calculating the derivatives (odefun) with some additional
% parameters la, kb(1..np), kk0 that are used only by odefun.
%
% output : the estimated vector for the next step yout(1..nv, 1..np) and number of function evaluations (nfe)

% D Vangheluwe 8 mrt 2005, revised 2 may 2005
% remark 1 : we did not take nv = number of variables as input argument as this can also
%   be obtained from the length of vector y or dy.

nv = size(y,1);
nfe = 0;
h = htot/nstep;
ym = y;
hm = repmat(h, nv, 1);
yn = y + hm .* dy;
x = xs + h;
yout = feval(odefun, x, y, la, kb, kk0);
nfe = nfe + 1;
h2 = 2 * hm;
for i = 2: nstep
   swap = ym + h2 .* yout;
   ym = yn;
   yn = swap;
   x = x + h;
   yout = feval(odefun, x, yn, la, kb, kk0);
   nfe = nfe + 1;
end
yout = 0.5 * (ym + yn + hm .* yout);
return

author: Sable

repo: mcbench-benchmarks


text:
%mrLoadAlignParams
%
%Loads in the file 'AlignParams.mat' from the current directory.
%If the file is not found, 'VolParams.mat' is loaded, and a new
%file 'AlignParams.mat' is created.  This way, mrLoadAlignParams
%can be used to seamlessly update 'VolParams' to 'AlignParams' to
%incorporate the switch from mrLoadVol to mrAlign.

%7/28/97 gmb   Wrote it.

ipSkip = 0;
curInplane = 0;
aTheta = 0;
cTheta = 0;

qt = '''';
if check4File('AlignParams')
  load AlignParams
else
  disp(['File ',qt,'AlignParams.mat',qt,' not found.']);
  if check4File('VolParams')
    disp(['Loading file ',qt,'VolParams.mat',qt]);
    disp(['and creating file ',qt,'AlignParams.mat',qt,'.']);
    load VolParams		    %holds curSag obXM obYM subject inplane_pix_size
    ipThickness = 1/inplane_pix_size(3);
    mrSaveAlignParams(obXM,obYM,subject,inplane_pix_size,...
	ipThickness,ipSkip,curSag,curInplane,aTheta,cTheta);
  else
    disp(['File ',qt,'VolParams.mat',qt,' not found.']);
    disp(['No files were loaded.']);
    return
  end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%kfunction ' '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros function.pane file
%
% Parameters:
% String: expr 'f(w,h,d,t,e) =', default: 'w + h': 'function expression'
% Integer: wsize ' ', default: 350: 'Width size of value data'
% Integer: hsize ' ', default: 350: 'Height size of value data'
% Integer: dsize ' ', default: 1: 'Depth size of value data'
% Integer: tsize ' ', default: 1: 'Time size of value data'
% Integer: esize ' ', default: 1: 'Elements size of value data'
%
% Example: kfunction( {'expr','w + h';'wsize',350;'hsize',350;'dsize',1;'tsize',1;'esize',1})
%
% Khoros helpfile follows below:

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [gIX,rankscore] = RankByStimLock(hfig,gIX)
fishset = getappdata(hfig,'fishset');

isStimAvr = getappdata(hfig,'isStimAvr');
isRawtime = getappdata(hfig,'isRawtime');
if isStimAvr == 1 || isRawtime == 1
    setappdata(hfig,'isStimAvr',0);
    setappdata(hfig,'isRawtime',0);
    global h_isStimAvr h_israwtime; %#ok<TLEV>
    h_isStimAvr.Value = 0;
    h_israwtime.Value = 0;
    UpdateTimeIndex(hfig);
end

[gIX, numU] = SqueezeGroupIX(gIX);
M = getappdata(hfig,'M');
C = FindClustermeans(gIX,M);

[~,~,H] = GetTrialAvrLongTrace(hfig,C);

% if fishset == 1
%     [~,~,H] = GetTrialAvrLongTrace(hfig,C);
% else
%     [~,~,H_multistim] = GetTrialAvrLongTrace(hfig,C);
%     if isempty(H_multistim)
%         errordlg('chosen stimulus range not suitable for stim-lock analysis');
%         rankscore = 1:numU;
%         return;
%     end
%     H = zeros(size(H_multistim,1),1);
%     for i = 1:size(H_multistim,1)
%         H(i) = min(H_multistim(i,:));
%     end
%
%     %     % instead of algebraic average along 2nd dimension, use
%     %     % inverse of geometric average... large value~low variation. geometric
%     %     % mean biases towards large values, i.e. good stim-lock of any stimulus
%     %     % is emphasized.
%     %     H = zeros(size(H_raw,1),1);
%     %     for i = 1:size(H_raw,1),
%     %         temp = sum((1./H_raw(i,:)).^2);
%     %         H(i) = 1./sqrt(temp);
%     %     end
%
% end

[gIX,rankscore] = SortGroupIXbyScore(H,gIX,numU);
end

author: xiuyechen

repo: FishExplorer


text:
function done = alreadydone(variable,method,goal_vexity)
global ALREADY_MODELLED ALREADY_MODELLED_INDEX
done = 0;
index = find(ALREADY_MODELLED_INDEX == variable);
if ~isempty(index)
    if ~isempty(ALREADY_MODELLED{index})
        if isequal(ALREADY_MODELLED{index}.method,method)
            if isequal(ALREADY_MODELLED{index}.goal_vexity,goal_vexity)
                done = 1;
            end
        elseif isequal(ALREADY_MODELLED{index}.method,'milp') || isequal(ALREADY_MODELLED{index}.method,'integer') || isequal(ALREADY_MODELLED{index}.method,'exact')
            done = 1;
        end
    end
end


text:
%% triSurf2ImSpec
% Below is a demonstration of the features of the |triSurf2ImSpec| function

%%
clear; close all; clc;

%% Syntax
% |[M,G]=triSurf2ImSpec(F,V,voxelSize);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%kktranspose 'Transpose Data on Dimensions'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ktranspose.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'input data object'
% OutputFile: o 'Output', required: 'resulting data object'
% Toggle: w ' ', default: 0: 'swap width and one other dimension'
% Toggle: h ' ', default: 0: 'swap height with one other dimension'
% Toggle: e ' ', default: 0: 'swap elements with one other dimension'
% Toggle: d ' ', default: 0: 'swap depth with one other dimension'
% Toggle: t ' ', default: 0: 'swap time with one other dimension'
%
% Example: o = kktranspose(i, {'i','';'o','';'w',0;'h',0;'e',0;'d',0;'t',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% ktranspose - Transpose Data across Dimensions
%
%  DESCRIPTION
% ktranspose will swap data across two dimensions.  Any of the polymorphic
% dimensions of WIDTH (W), HEIGHT (H), ELEMENTS (E), DEPTH (D),
% and TIME (T), can be switched with another.
%
% \f(CW
%
%                  'depth' planes of value data
%                ____/   ____       ____
%               +---+|  +---+|     +---+|  The value segment is a time
%     'height'--|   ||  |   || ... |   ||  series of volumes.
%               +---+   +-|-+      +---+
%                 |       |
%              'width'   each element in the value value
%                        segment is a vector of size
%                        'elements'.
%
%
% The dimensions which should be transposed are specified.
% On the graphical user interface accessible via
% cantata, these flags are represented as a series of buttons, with one
% for each dimension.  On the command line, these flags are specified
% as a letter : w,h,d,t,e.   Only two flags should be specified. The
% dimensions indicated by these two flags will be transposed. For example,
%
% 	ktranspose -i infile -o outfile -w -h
%
% will cause the height dimension and the width dimension to be
% switched. If the input file has width = 4, height = 2,
% the output dimension will have width = 2, height = 4.
%
%
%  "Map Data"
% If the source object has map data, the map data cannot
% have depth, time, or elements dimensions of greater than
% one.  The map width and height dimensions can be any size.
%
%  "Location & Time Data"
% If the source object has rectilinear or curvilinear
% location or time data, the orientation operator will fail.  If the
% source object has uniform location data, orientation will
% continue, but the location data will not be modified
% to reflect the reorientation.
%
%  "Map Data"
% If the object has mask data, it will be transposed along with the value data.
%
%  "Failure Modes"
% The transpose operator will fail if the source object does not contain
% value data, or if it contains rectilinear or curvilinear location or time
% data.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% kaxis, korient
%
%  RESTRICTIONS
% This program has not yet been updated to completely support the
% polymorphic data model.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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


text:
function [w,run] = train_cg2(x,w)
% Conjugate gradient method
% Data is columns of x, each column already scaled by the output y (+1 or -1).
% w is the starting guess for the parameters (a column).

% Written by Thomas P Minka

flops(0);
old_g = zeros(size(w));
old_u = [];
for iter = 1:1000
  old_w = w;
  % s1 = 1-sigma
  s1 = 1./(1+exp(w'*x));
  g = x*s1';
  a = s1.*(1-s1);
  if iter == 1
    u = g;
  else
    h = x*diag(a)*x';
    u = orthA([old_u g],h);
    u = u(:,cols(u));
  end

  % line search along u
  ug = u'*g;
  ux = u'*x;
  uhu = (ux.^2)*a';
  w = w + (ug/uhu)*u;
  old_g = g;

  % update memory
  old_u = [old_u u];
  if(cols(old_u) > 1)
    old_u(:,1) = [];
  end

  fl = flops;
  run.w(:,iter) = w;
  run.flops(iter) = fl;
  run.e(iter) = logProb(x,w);
  flops(fl);

  if max(abs(w - old_w)) < 1e-5
    break
  end
end
figure(2)
plot(run.e)
if iter == 1000
  warning('not enough iters')
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function [V,F] = cookie_cutter(P,E)
  % COOKIE_CUTTER Generate a cookie cutter model from a curve
  %
  % [V,F,W,G] = cookie_cutter(P,E)
  %
  % Inputs:
  %   P  #P by 2 list of positions (bbd should be ~40)
  %   E  #E by 3 list of edge indices into P
  % Outputs:
  %   V  #V by 3 list of positions of cutter part
  %   F  #F by 3 list of triangle indices on cutter part
  %
  % Example:
  %
  %     clf;
  %     axis([-40 40 -40 40]);axis equal;
  %     P = get_pencil_curve();
  %     E = [1:size(P,1);2:size(P,1) 1]';
  %     [V,F] = cookie_cutter(P,E);
  %     tsurf(F,V);axis equal;
  %
  function [V,O] = robust_overlap(S,E)
    sqr_len = mean(sum((S(E(:,1),:)-S(E(:,2),:)).^2,2));
    max_area = sqr_len/9;
    [V,F] = cdt(S,E,'TriangleFlags',sprintf('-qYY -a%g',max_area),'Quiet',true);
    w = winding_number(S,E,barycenter(V,F));
    F = F(w>0.5,:);
    [V,I] = remove_unreferenced(V,F);
    F = I(F);
    O = outline(F);
    [V,I] = remove_unreferenced(V,O);
    O = I(O);
    [V,~,I] = remove_duplicate_vertices(V,1e-3);
    O = I(O);
    % remove degnerate edges
    O = O(O(:,1) ~= O(:,2),:);
    % remove duplice edges
    [~,I] = unique(sort(O,2),'rows');
    O = O(I,:);
    assert(size(unique(sort(O,2),'rows'),1) == size(O,1));
    [V,F] = triangle(V,O,[],'Flags','-YY','Quiet');
    w = abs(winding_number(V,O,barycenter(V,F)));
    F = F(w>0.5,:);
    O = outline(F);
    [V,I] = remove_unreferenced(V,O);
    O = I(O);
    assert(size(unique(sort(O,2),'rows'),1) == size(O,1));

  end

  % thin wall
  offset = 0.5;
  [SI,SO,EI,EO] = offset_curve(P,offset);
  S = [SI;SO];
  E = [EI;size(SI,1)+EO];
  [S,E] = robust_overlap(S,E);

  % thick wall
  k_offset = 2;
  [KI,KO,KEI,KEO] = offset_curve(P,k_offset);
  K = [KI;KO];
  KE = [KEI;size(KI,1)+KEO];
  [K,KE] = robust_overlap(K,KE);

  % height of cutting wall
  h = 15;
  % height of lip
  th = 3;
  % height of stamper
  sth = 5;
  [CV,CF] = wedding_cake({K,S},{KE,E},[th,h]);

  tol = 0.8;
  [sP,~,sE,~] = offset_curve(P,offset+tol);
  sE = fliplr(sE);
  [sP,sE] = robust_overlap(sP,sE);
  % center
  c = mean(sP);
  % center should be inside
  assert(abs(winding_number(sP,sE,c))>0.5,'Uhoh, center for post fell outside curve. Too non-convex');
  % radius determined by closest point on curve
  r = min(pdist2(K,c));
  r = 0.9*r;
  r = min(r,4);
  theta = linspace(0,2*pi,100)';
  theta = theta(1:end-1);
  hP = bsxfun(@plus,c,r*[cos(theta) sin(theta)]);
  hE = [1:size(hP,1);2:size(hP,1) 1]';
  [SV,SF] = wedding_cake({sP,hP},{sE,hE},[sth,2*h]);

  % concatenate
  F = [CF;size(CV,1)+SF];
  %V = [CV;bsxfun(@plus,[max(CV(:,1))*2 0 0],SV)];
  % Cutter needs to be flipped!
  V = [bsxfun(@times,[-1 1 1],CV);bsxfun(@plus,[abs(max(CV(:,1))-min(CV(:,1)))*1.1 0 0],SV)];
end

author: alecjacobson

repo: gptoolbox


text:
function CP = copula(gamma)

% copula - Compute copular associated to a joint distribution gamma.
%
% CP = copula(gamma)
%
%   Copyright (c) 2017 Gabriel Peyre

N = size(gamma,1);

mu = sum(gamma,2);
nu = sum(gamma,1)';

% cumulative
cmu = [0;cumsum(mu)];
cnu = [0;cumsum(nu)];
% inverse cumulatives
t = (0:N)'/N;
icmu = interp1(cmu, t, t, 'spline');
icnu = interp1(cnu, t, t, 'spline');
%
Cgamma = cumsum(cumsum(gamma,1),2);
s = linspace(0,1,N);

icmu = max(min(icmu,1),0);
icnu = max(min(icnu,1),0);

[S2,S1] = meshgrid(s,s);
[T2,T1] = meshgrid(icnu,icmu);
CP = interp2(S2,S1,Cgamma,T2,T1);

end


text:
function [data] = readmarkerfile(folder);
% function reads a CTF data file (.ds), Read 'markerfile.mrk'
% marker is the read-only 'markerfile.mrk'
% creates structure 'data' which contains number_markers, number_samples,
%       marker_names, trial_times, marker.

%
%
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> %
%      <                                                      > %
%      <                    DISCLAIMER:                       > %
%      <                                                      > %
%      <  THIS PROGRAM IS INTENDED FOR RESEARCH PURPOSES ONLY. > %
%      <  THIS PROGRAM IS IN NO WAY INTENDED FOR CLINICAL OR   > %
%      <                    OFFICIAL USE.                     > %
%      <                                                      > %
%      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<> %
%
%
%

a = exist(fullfile(folder,'MarkerFile.mrk'));
if a == 0
  disp(' ');
  disp('Error: Your data could not be found. MarkerFile.mrk does not exist in this directory.');
  disp(' ');
  data = [];
  return
end   % error message if incorrect directory is displayed

markr =  textread(fullfile(folder,'MarkerFile.mrk'),'%s','delimiter','\n');
% creates markerfile.mrk in FOLDER directory. Type 'marker' to view.

% define markers
number_markers_id = strmatch('NUMBER OF MARKERS:',markr,'exact');
markers = str2num(markr{number_markers_id+1});

% define samples
number_samples_id = strmatch('NUMBER OF SAMPLES:',markr,'exact');
samples = str2num(char(markr(number_samples_id+1)));

for i = 1:length(samples)
  if samples(i) == 0
    disp('Error: one of the markers is bad.');
    data = [];
    return
  end
end

name_id = strmatch('NAME:',markr,'exact');
names = markr(name_id+1);
% defines names

trial = strmatch('LIST OF SAMPLES',markr)+2;

for i = 1:markers
  trials{i} = str2num(char(markr(trial(i):trial(i)+samples(i))));
  trials{i}(:,1) = trials{i}(:,1) + 1;
end

data = struct('number_markers',{markers},'number_samples',{samples},'marker_names',{names},'trial_times',{trials});
% defines data

return

author: PatternRecognition

repo: OpenBMI


text:
% test for heat diffusion on meshes

% load mesh
name = 'mushroom';
name = 'fandisk';
name = 'lion-head';
name = 'bunny';
name = 'skull';
name = 'horse';
name = 'nefertiti';
name = 'elephant-50kv';
options.name = name;
[vertex,face] = read_off([name '.off']);

rep = 'results/mesh-smoothing/';
if not(exist(rep))
    mkdir(rep);
end

nvert = size(vertex,2);
nface = size(face,2);

laplacian_type = 'distance';
laplacian_type = 'combinatorial';
laplacian_type = 'conformal';

%% symmetric laplacian
if not(strcmp(laplacian_type,'conformal'))
    options.symmetrize = 1;
    options.normalize = 0;
    L0 = compute_mesh_laplacian(vertex,face,laplacian_type,options);
    G0 = compute_mesh_gradient(vertex,face,laplacian_type,options);
    disp(['Error (should be 0): ' num2str(norm(L0-G0'*G0, 'fro')) '.']);
    options.normalize = 1;
    L1 = compute_mesh_laplacian(vertex,face,laplacian_type,options);
    G1 = compute_mesh_gradient(vertex,face,laplacian_type,options);
    disp(['Error (should be 0): ' num2str(norm(L1-G1'*G1, 'fro')) '.']);
end

%% un-symmetric laplacian
options.symmetrize = 0;
options.normalize = 1;
disp('--> Computing laplacian');
L = compute_mesh_laplacian(vertex,face,laplacian_type,options);

%% heat diffusion flow
Tlist = [0 10 40 200];
options.dt = 0.3;
clf;
for i=1:length(Tlist)
    options.Tmax = Tlist(i);
    vertex1 = perform_mesh_heat_diffusion(vertex,face,L,options);
    % display
    subplot(1,length(Tlist),i);
    plot_mesh(vertex1,face,options);
    shading interp; camlight; axis tight;
end
saveas(gcf, [rep name '-smoothing-' laplacian_type '.png'], 'png');

%% quadratic regularization
clf;
for i=1:length(Tlist)
    At = speye(nvert) + Tlist(i)*L;
    vertex1 = (At\vertex');
    % display
    subplot(1,length(Tlist),i);
    plot_mesh(vertex1,face,options);
    shading interp; camlight; axis tight;
end
saveas(gcf, [rep name '-quadreg-' laplacian_type '.png'], 'png');

%% non linear L1 minimization
options.symmetrize = 1;
options.normalize = 0;
disp('--> Computing gradient');
G = compute_mesh_gradient(vertex,face,laplacian_type,options);
Tlist = [0 0.01 0.02 0.05 0.1 0.2];
clf;
for i=1:length(Tlist)
    t = Tlist(i);
    vertex1 = vertex;
    options.niter = 500;
    options.lambda_max = t*2;
    options.lambda_min = t;
    if t>0
    for k=1:3
        [a,tmp,E] = perform_analysis_regularization(vertex(k,:)', G, options);
        vertex1(k,:) = a(:)';
    end
    end
    % display
    subplot(2,ceil(length(Tlist)/2),i);
    plot_mesh(vertex1,face,options);
    shading interp; camlight; axis tight;
end
saveas(gcf, [rep name '-l1reg-' laplacian_type '.png'], 'png');

author: gpeyre

repo: matlab-toolboxes


text:
function K=kernel_noise(M,dimker)

% Copyright (C) <2007>  <Francesc Moreno-Noguer, Vincent Lepetit, Pascal Fua>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the version 3 of the GNU General Public License
% as published by the Free Software Foundation.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Francesc Moreno-Noguer, CVLab-EPFL, September 2007.
% fmorenoguer@gmail.com, http://cvlab.epfl.ch/~fmoreno/

MtM=M'*M;
[V,S]=eig(MtM);

K=V(:,dimker:-1:1);

author: cvlab-epfl

repo: EPnP

save_path: github-repos/MATLAB/cvlab-epfl-EPnP


text:
function fismat = genfis4(Xin, Xout, fistype, sig_par, varargin)

%GENFIS4 Generates a FIS using CART
%
%   Given separate sets of input and output data, GENFIS4 generates a fuzzy
%   inference system (FIS) using CART algorithm. GENFIS4 accomplishes this
%   by extracting a set of rules that models the data behavior. The rule
%   extraction method first uses TREEINFO to determine the number
%   of rules and membership functions for the antecedents and consequents.
%
%   FIS = GENFIS4(XIN, XOUT) returns a Sugeno-type FIS given input data XIN
%   and output data XOUT. The matrices XIN and XOUT have one column per FIS
%   input and output, respectively.
%
%   FIS = GENFIS4(XIN, XOUT, TYPE) returns FIS of type specified by the
%   argument TYPE. It can take one of two values: 'mamdani' or 'sugeno'.
%
%   FIS = GENFIS4(XIN, XOUT, TYPE, SIGM_PAR) allows you to specify the
%   sigmoid parameters in the argument SIGM_PAR. The parameter is used to
%   form input membership function. The larger SIGM_PAR the closer
%   behaviour of fuzzy CART to behaviour of "crisp" CART. It should be
%   positive scalar or vector. The latter should have size equal to number
%   of FIS inputs. This argument also takes the value 'auto' in which case
%   GENFIS4 uses additional information containing in data set to calculate
%   SIGM_PAR.
%
%   FIS = GENFIS4(XIN, XOUT, TYPE, SIGM_PAR, VARARGIN) allows you to
%   specify options for the CART algorithm. Type HELP CLASSREGTREE and
%   HELP TEST for a list of options that can be specified for the CART
%   algorithm.
%
%   Examples:
%
%       Xin1 = 7 * rand(50, 1);
%       Xin2 = 20 * rand(50, 1) - 10;
%       Xin = [Xin1 Xin2];
%       Xout = 5 * rand(50, 1);
%       fis = genfis4(Xin, Xout);
%
%       fis = genfis4(Xin, Xout, 'mamdani', [10, 1]);
%       specifies the type of FIS and the sigmoid parameters desired.
%
%       fis = genfis4(Xin, Xout, 'mamdani', 'auto', {'minparent', 15});
%       specifies the type of FIS, the sigmoid parameters desired and CART
%       options.
%
%   See also CLASSREGTREE, TEST, TREEINFO, GENFIS3, ANFISX

%   Per Konstantin A. Sidelnikov, 2009.

%%%%%%%%%%%%%%%%%%%%
% Some constants
%%%%%%%%%%%%%%%%%%%%

% Number of standard deviations (can be changed)
NSTD = 1;
% Area under the gauss curve in the
% [mu - NSTD * sigma, mu + NSTD * sigma]
AREA = erf(NSTD / sqrt(2));
% Scaling factor
SCALE = log((1 + AREA) / (1 - AREA)) / NSTD;

%%%%%%%%%%%%%%%%%%%%
% Number of input arguments checking
%%%%%%%%%%%%%%%%%%%%

if nargin < 2
    error('FuzzyLogic:missingparams', ...
        'genfis4 requires input and output data to build a FIS.');
end

if nargin < 3
    fistype = 'sugeno';
end

if nargin < 4
    sig_par = 'auto';
end

% hardcoded for now
in_mftype = 'sigmf';
out_mftype = 'gaussmf';

%%%%%%%%%%%%%%%%%%%%%%%%%
% IO checking
%%%%%%%%%%%%%%%%%%%%%%%%%

% Check Xin, Xout
[numData, numInp] = size(Xin);
[numData2, numOutp] = size(Xout);

if numData ~= numData2
    % There's a mismatch in the input and output data matrix dimensions
    if numData == numOutp
        % The output data matrix should have been transposed, we'll fix it
        Xout = Xout';
        numOutp = numData2;
    else
        error('FuzzyLogic:dimensionmismatch', ...
            'Mismatched input and output data matrices.');
    end
end

if numOutp > 1
    error('FuzzyLogic:outputmismatch', ...
        'Output data must be a vector.');
end

%%%%%%%%%%%%%%%%%%%%
% Checking sig_par
%%%%%%%%%%%%%%%%%%%%

% Convert sig_par (if necessary) to numeric array
arrayNaN = NaN(1, numInp);
if ~isnumeric(sig_par)
    if ~isequal(sig_par, 'auto')
        error('FuzzyLogic:sigmoid', ...
            'Set sigmoid parameter to ''auto'' or a value greater than 0');
    end
    sig_par = arrayNaN;
elseif ~isscalar(sig_par)
    tmp = arrayNaN;
    tmp(1 : length(sig_par)) = sig_par;
    sig_par = tmp;
else
    sig_par = repmat(sig_par, 1, numInp);
end
% Avoid negative sig_par
sig_par = abs(sig_par);
% Length of sig_par must be the same as number of inputs
if length(sig_par) ~= numInp
    error('FuzzyLogic:sigmoid', ...
        ['Set number of sigmoid parameters to ', ...
        '1 or %d (number of inputs)'], numInp);
end

%%%%%%%%%%%%%%%%%%%%
% Creating tree
%%%%%%%%%%%%%%%%%%%%
ti = treeinfo(Xin, Xout, varargin{:});

n_node = numel(ti.node);
n_leaf = numel(ti.leaf);

in_mf = zeros(1, n_node);

%%%%%%%%%%%%%%%%%%%%%%
% Building FIS
%%%%%%%%%%%%%%%%%%%%%

% Initialize a FIS
str = sprintf('%s%g%g', fistype, numInp, numOutp);
fismat = newfis(str, fistype);

% Loop through and add inputs
for ind = 1 : numInp
    fismat = addvar(fismat, ...
        'input', ['in' num2str(ind)], minmax(Xin(:, ind)'));
end

% Loop through and add mf's
for ind = 1 : n_node
    var = ti.node(ind).variable;
    cut = ti.node(ind).cutpoint;
    s = ti.node(ind).sample;
    mfparams = computemfparams(in_mftype, Xin(s, var), ...
        cut, sig_par(var), SCALE);

    fismat = addmf(fismat, 'input', var, ...
        ['larger', num2str(ind)], in_mftype, mfparams);

    in_mf(ind) = length(fismat.input(var).mf);
end

% Add output
fismat = addvar(fismat, 'output', 'out', minmax(Xout'));

switch fistype
    case 'sugeno'
        % Loop through and add mf's
        for ind = 1 : n_leaf
            s = ti.leaf(ind).sample;
            mfparams = computemfparams('linear', Xin(s, :), Xout(s));

            fismat = addmf(fismat, 'output', 1, ...
                ['class', num2str(ind)], 'linear', mfparams);
        end
    case 'mamdani'
        % Loop through and add mf's
        for ind = 1 : n_leaf
            s = ti.leaf(ind).sample;
            mfparams = computemfparams(out_mftype, [], Xout(s));

            fismat = addmf(fismat, 'output', 1, ...
                ['class', num2str(ind)], out_mftype, mfparams);
        end
    otherwise
        error('FuzzyLogic:unknownfistype', ...
            'Unknown fistype specified');
end

% Create rules
rulelist = cell(n_leaf, 4);
for ind = 1 : n_leaf
    n = ti.branch(ind).nodes;
    e = ti.branch(ind).ineqs;

    var = [ti.node(n).variable];

    rulelist{ind, 1} = [var; in_mf(n) .* e];
    rulelist{ind, 2} = ind;
    rulelist{ind, 3} = 1;
    rulelist{ind, 4} = 1;
end

fismat = addrulex(fismat, rulelist);

function mfparams = computemfparams(mf, X, y, sig_par, scale)
%   This subfunction computes parameters of input and
%   output membership functions dependeding on mf's value.

switch mf
    case 'sigmf'
        % NaN's value of sig_par means its automatic calculation
        if isnan(sig_par)
            sigma = std(X);
            sig_par = scale / sigma;
        end
        mfparams = [sig_par, y];
    case 'gaussmf'
        sigma = std(y);
        % Check if y is a scalar or consists of identical elements
        if sigma == 0
            sigma = sqrt(eps);
        end
        c = mean(y);
        mfparams = [sigma, c];
    case 'linear'
        numData = size(X, 1);
        A = [X, ones(numData, 1)];
        % Using pinv instead of ldivide avoids
        % warning if A is close to singular
        mfparams = (pinv(A) * y)';
    otherwise
        error('FuzzyLogic:invalidmftype', ...
            'Unknown type of membership function specified.');
end

function pr = minmax(p)
%MINMAX Ranges of matrix rows.
%
%  Syntax
%
%    pr = minmax(p)
%
%  Description
%
%    MINMAX(P) takes one argument,
%      P - RxQ matrix.
%    and returns the Rx2 matrix PR of minimum and maximum values
%    for each row of P.
%
%    Alternately, P can be an MxN cell array of matrices.  Each matrix
%    P{i,j} should have Ri rows and Q columns.  In this case, MINMAX returns
%    an Mx1 cell array where the mth matrix is an Rix2 matrix of the
%    minimum and maximum values of elements for the matrics on the
%    ith row of P.
%
%  Examples
%
%    p = [0 1 2; -1 -2 -0.5]
%    pr = minmax(p)
%
%    p = {[0 1; -1 -2] [2 3 -2; 8 0 2]; [1 -2] [9 7 3]};
%    pr = minmax(p)

if iscell(p)
    m = size(p, 1);
    pr = cell(m, 1);
    for i = 1 : m
        pr{i} = minmax([p{i, :}]);
    end
elseif isa(p, 'double')
    pr = [min(p, [], 2), max(p, [], 2)];
else
    error('Argument has illegal type.');
end


text:
function [cost_total,grad_total] = RAECost(theta, alpha_cat, cat_size, beta, dictionary_length, hiddenSize, ...
    lambda, We_orig, data_cell, labels, freq_orig, sent_freq, f, f_prime)

[~, ~, ~, ~, ~, ~, ~, Wcat, bcat, We] = getW(1, theta, hiddenSize, cat_size, dictionary_length);

szWe = length(We(:));
szbcat = length(bcat(:));
szWcat = length(Wcat(:));
theta1 = theta;
theta1(end-szWe+1:end) = We;
theta1(end-szWe-szbcat-szWcat+1:end-szWe) = [];

theta2 = theta;
theta2(end-szWe+1:end) = We_orig(:) + We(:);
lambda2 = lambda;
lambda(3) = lambda2(4);
lambda(4) = lambda2(3);

% disp('DEBUGGING: DELETE AFTERWARDS')
% data_cell=data_cell(1:2);
% update W using Greedy Unsupervised RAE
[costRAE, gradRAE, allKids] = computeCostAndGradRAE([], theta1, 0, alpha_cat, cat_size, beta, dictionary_length, hiddenSize, ...
    (alpha_cat)*lambda, We_orig , data_cell, labels, freq_orig, f, f_prime);

WegradRAE = gradRAE(end-szWe+1:end);
gradRAE(end-szWe+1:end) = 0;
gradRAE = [gradRAE; zeros(szbcat+szWcat,1)];
gradRAE(end-szWe+1:end) = WegradRAE;

cost_total =  costRAE;
grad_total =  gradRAE;

author: jacoxu


text:
%% patch2EdgeIm
% Below is a demonstration of the features of the |patch2EdgeIm| function

%%
clear; close all; clc;

%% Syntax
% |[varargout]=patch2EdgeIm(varargin);|

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


text:
function [T,TI,E] = polygons_to_triangles(PI,PC)
  % POLYGONS_TO_TRIANGLES Given a polygonal mesh, replace each polygon with a
  % fan of triangles.
  %
  % [T,TI,E] = polygons_to_triangles(PI,PC)
  %
  % Inputs:
  %   PI  #PI stream of polygon indices into rows of some V
  %   PC  #polys+1 list of cumulative sum of dual face valences
  % Outputs:
  %   T  #T by 3 list of triangle indices into rows of V
  %   TI  #TI list of indices into input faces (1:#polys)
  %   E  #E by 2 list of edge indices into rows of V
  %
  % See also: dual
  I1 = repelem(PI(PC(1:end-1)+1),diff(PC)-2);
  I2 = PI;
  % remove first and last entries
  I2([PC(1:end-1)+1;PC(2:end)]) = [];
  I3 = PI;
  % remove first two entries
  I3([PC(1:end-1)+1;PC(1:end-1)+2]) = [];
  T = [I1 I2 I3];
  TI = repelem(1:numel(PC)-1,diff(PC)-2);
  if nargout>2
    I1 = PI;
    % order index per polygon
    Pj = [(1:numel(PI))' - repelem(PC(1:end-1),diff(PC))];
    % increment by one mod size
    Pj = mod(Pj,repelem(diff(PC),diff(PC)))+1;
    I2 = PI(repelem(PC(1:end-1),diff(PC)) + Pj);
    E = [I1 I2];
  end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function im_rec = tomo_recon_fbp(projmat_bp,angles)
%TOMO_RECON_FBP   Tomographic reconstruction using Filtered Back Projection
%   (BP) algorithm.
%
%   Phymhan
%   05-Aug-2013 17:23:35

[n_proj,D] = size(projmat_bp);
if length(angles) ~= n_proj
    fprintf('Size of proj_mat is incorrect.\r')
    im_rec = [];
    return
end
im_rec = zeros(D);
for k = 1:n_proj
    proj_vec = hpf(projmat_bp(k,:));
    im_rot = repmat(proj_vec/D,D,1);
    im_rot = imrotate(im_rot, angles(k), 'bilinear', 'crop');
    %DEBUG
    % image(im_rot);
    % waitforbuttonpress
    % %
    im_rec = im_rec+im_rot/n_proj;
    %DEBUG
    % imagesc(im_rec)
    % waitforbuttonpress
    % %
end
end


text:
function [W] = Btuph2W(Btuph)
% Convert power from British thermal units per hour to watts.
% Chad A. Greene 2012
W = Btuph*0.29307107;


text:
im = iread('lena.pgm');
whos im
about im
figure(1); idisp(im)
im(100,120)
im(100:104,120:124)
figure(2); plot(im(200,:))
eye = iroi(im);

im = iread('lena.pgm', 'double');
about im
idisp(im)

[im,tags]=iread('/Users/corkep/proj/Oxford/2012/images/lumix/P1000545.JPG');
idisp(im)
tags
tags.DigitalCamera

author: petercorke

repo: machinevision-toolbox-matlab


text:
function value = month_length_islamic ( y, m )

%*****************************************************************************80
%
%% MONTH_LENGTH_ISLAMIC returns the number of days in an Islamic month.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year in which the month occurred.
%
%    Input, integer M, the number of the month.
%
%    Output, integer MONTH_LENGTH_ISLAMIC, the number of days
%    in the month.
%
  mdays = [ 30, 29, 30, 29, 30, 29, 30, 29, 30, 29, 30, 29 ];
%
%  Copy the input.
%
  m2 = m;
  y2 = y;
%
%  Check the input.
%
  [ y2, m2, ierror ] = ym_check_islamic ( y2, m2 );

  if ( ierror ~= 0 )
    month_length_islamic = 0;
    return
  end
%
%  Get the number of days in the month.
%
  value = mdays(m2);
%
%  If necessary, add 1 day for a leap year.
%
  if ( m2 == 12 && year_is_leap_islamic ( y2 ) )
    value = value + 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout=m2v(varargin)
%
% vol=m2v(node,face,Nxyz)
%  or
% vol=m2v(node,face,xi,yi,zi)
%
% shortcut for mesh2vol, rasterizing a teterahedral mesh to a volume using graphics
%
% author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
%
% input/output: please see details in the help for mesh2vol
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

[varargout{1:nargout}]=mesh2vol(varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function dat= proc_variance(dat, nSections, calcStd)
%PROC_VARIANCE - computes the variance in equally spaced intervals
%
%Synopsis
% dat= proc_variance(dat, <nSections=1, calcStd=0>)
%
% IN   dat       - data structure of continuous or epoched data
%      nSections - number of intervals in which var is to be calculated
%      calcStd   - standard deviation is calculated instead of variance
%
% OUT  dat       - updated data structure
%
%Description
% calculate the variance in 'nSections' equally spaced intervals.
% works for cnt and epo structures.

if nargin==0,
  dat=[];return
end

misc_checkType(dat, 'STRUCT(x)');
misc_checkTypeIfExists('nSections', 'INT');
misc_checkTypeIfExists('calcStd','BOOL');
dat = misc_history(dat);

if ~exist('nSections','var'), nSections=1; end
if ischar(nSections) && strcmpi(nSections,'std'),
  nSections=1;
  calcStd=1;
end
if ~exist('calcStd','var') || (ischar(calcStd) && strcmpi(calcStd,'var')),
  calcStd= 0;
end
if ischar(calcStd) && strcmpi(nSections,'std'),
  calcStd=1;
end

[T, nChans, nMotos]= size(dat.x);
inter= round(linspace(1, T+1, nSections+1));
dat.t = [] ;

xo= zeros(nSections, nChans, nMotos);
for s= 1:nSections,
  Ti= inter(s):inter(s+1)-1;
  if length(Ti)==1,
    warning('calculating variance of scalar');
  end
  if calcStd,
    xo(s,:,:)= reshape(std(dat.x(Ti,:),0,1), [1, nChans, nMotos]);
  else
    if length(Ti)==1,
      xo(s,:,:)= dat.x(Ti,:);
    else
      if nChans*nMotos*length(Ti)<=10^6;
        xo(s,:,:)= reshape(var(dat.x(Ti,:)), [1, nChans, nMotos]);
      else
        for i=1:nMotos
          xo(s,:,i) = var(dat.x(Ti,:,i));
        end
      end
    end
  end
  dat.t(s) = Ti(end);
end

dat.x= xo;

author: bbci


text:
classdef ECGtask_Delineation_corrector < ECGtask

% ECGtask for ECGwrapper (for Matlab)
% ---------------------------------
%
% Description:
%
% Abstract class for defining ECGtask interface
%
%
% Author: Mariano Llamedo Soria (llamedom at frba.utn.edu.ar)
% Version: 0.1 beta
% Birthdate  : 18/2/2013
% Last update: 18/2/2013

    properties(GetAccess = public, Constant)
        name = 'Delineation_corrector';
        % Require the parent Wrapper object to this task
        target_units = 'Wrapper';
        doPayload = true;
    end

    properties( GetAccess = public, SetAccess = private)
        % if user = memory;
        % memory_constant is the fraction respect to user.MaxPossibleArrayBytes
        % which determines the maximum input data size.
        % Size > 1 means that this task can handle big data structures, as
        % this case.
        memory_constant = realmax;
        started = false;
    end

    properties( Access = private, Constant)

        fig_hdl = 1;
        cAnnotationOutputFields = { 'Pon' 'P' 'Poff' 'QRSon' 'Q' 'R' 'S' 'QRSoff' 'Ton' 'T' 'Toff'};

    end

    properties( Access = private )

    end

    properties

        progress_handle
        caller_variable = 'payload';
        tmp_path
        payload

    end

    methods

        function obj = ECGtask_QRS_detection(obj)

            obj.wavedet_config.setup.wavedet.QRS_detection_only = true;

        end

        function Start(obj, ECG_header, ECG_annotations)

            obj.started = true;

        end

        function payload = Process(obj, ECG, ECG_start_offset, ECG_sample_start_end_idx, ECG_header, ECG_annotations, ECG_annotations_start_end_idx )

            payload = [];

            obj.progress_handle.hide()

%             if( ~obj.started )
%                 obj.Start(ECG_header);
%                 if( ~obj.started )
%                     cprintf('*[1,0.5,0]', 'Task %s unable to be started for %s.\n', obj.name, ECG_header.recname);
%                     return
%                 end
%             end

            this_start_end = (ECG_sample_start_end_idx + ECG_start_offset - 1);
            aux_val = this_start_end ./ ECG_header.freq;
%             disp_string_title(1, sprintf( 'Correcting from %s to %s', Seconds2HMS(aux_val(1)), Seconds2HMS(aux_val(2)) ) );

            if( isempty(obj.payload) )
                Ann_struct2 = ECG_annotations;
            else
                Ann_struct = obj.payload;
                AnnNames = {};
                kk = 1;

                % Algorithm name
                for alg_name = rowvec(fieldnames(Ann_struct))

                    for lead_name = rowvec(fieldnames(Ann_struct.(alg_name{1})))

                        for wave_chosen = rowvec(obj.cAnnotationOutputFields)

                            wave_name = wave_chosen{1};

                            if( isfield(Ann_struct.(alg_name{1}).(lead_name{1}), wave_name ) )
                                aux_val = Ann_struct.(alg_name{1}).(lead_name{1});
                                if( isfield(aux_val, wave_name) )
                                    aux_val2 = ( aux_val.(wave_name) - aux_val.qrs ) * 1/ECG_header.freq;
                                    bAux = aux_val.qrs >= this_start_end(1) & aux_val.qrs <= this_start_end(2);
                                    lead_str = strrep(lead_name{1}, ' ', '_' );
                                    aux_str = [ wave_name '_' lead_str '_' alg_name{1}  ];
                                    Ann_struct2.(aux_str).time = [ (colvec(aux_val.qrs(bAux)) - ECG_start_offset + 1)  colvec(aux_val2(bAux))];
                                    AnnNames(kk,:) = { aux_str 'time' };
                                    kk = kk + 1;
                                end
                            end

                        end

                    end
                end

                % info about the series, requiered by the corrector
                % software.
                Ann_struct2.series_quality.AnnNames = AnnNames;
                Ann_struct2.series_quality.ratios = zeros(size(AnnNames,1),1);
                Ann_struct2.series_quality.estimated_labs = [];

            end

            if( ~isempty(ECG_start_offset) || ECG_start_offset > 1)
                 [~, iHours, iMins, iSeconds, iMilli ] = Seconds2HMS(ECG_start_offset/ECG_header.freq);
                 ECG_header.btime = datestr([2014,1,1,iHours, iMins,iSeconds], 'HH:MM:SS');
            end

            QRScorrector('ECG', ECG, 'QRS_annotations', Ann_struct2, 'Figure', figure(obj.fig_hdl) );

            disp_string_framed('*Blue', 'User interaction required' );

            aux_str = ['figure ' num2str(obj.fig_hdl) ''];

            aux_str2 = 'F5 (Run)';

            fprintf(1, 'This ECGtask allow user interaction. Press [CTRL + G] in %s to save results and press %s to continue.\n', aux_str, aux_str2)
            keyboard

            % last chance to save results
            if( ishandle(obj.fig_hdl) && (isempty(payload) || ~isstruct(payload)) )
                disp_string_framed('*[1,0.5,0]', 'Payload variable not saved' );
                fprintf(1, 'Press [CTRL + G] in %s to save results and press %s to continue.\n', aux_str, aux_str2)
                keyboard

                if( ~isempty(payload) && isstruct(payload) )
                    % save work done, to further improve in following
                    % invocations.
                    disp_string_framed('*Magenta', 'Corrections saved' );

                else
                    disp_string_framed('*[1,0.5,0]', 'Data not saved' );
                end

            else
                disp_string_framed('*Magenta', 'Corrections saved' );
            end

            if(ishandle(obj.fig_hdl))
                delete(obj.fig_hdl)
            end

            if(~(isempty(payload) || ~isstruct(payload)))
                % move QRS corrections according to the input offset
                for fn = rowvec(fieldnames(payload))
                    if( isfield(payload.(fn{1}), 'time' ) )
                        payload.(fn{1}).time = payload.(fn{1}).time + ECG_start_offset - 1;
                    end
                end
            end

            obj.progress_handle.show()

        end

        function payload = Finish(obj, payload, ECG_header)

        end

        function payload = Concatenate(obj, plA, plB)

            if( isempty(plA) )

                payload = plB;

            else

                fields = rowvec(unique( [rowvec(fieldnames(plA)) rowvec(fieldnames(plB))] ));

                diff_fields = setdiff(fields, {'series_quality'});

                for fn = diff_fields

                    if( isfield(plA, fn{1}) && isfield(plB, fn{1}) )
                        payload.(fn{1}).time = [ colvec(plA.(fn{1}).time); colvec(plB.(fn{1}).time) ];

                    elseif( ~isfield(plA, fn{1}) && isfield(plB, fn{1}) )
                        payload.(fn{1}).time = colvec(plB.(fn{1}).time);
                    end

                end

                aux_idx = find(strcmpi(fields, 'series_quality'));

                if( ~isempty(aux_idx) )

                    if( isfield(plA, 'series_quality') && isfield(plB, 'series_quality') )
                        payload.series_quality.ratios = [plA.series_quality.ratios plB.series_quality.ratios];
                        payload.series_quality.estimated_labs = cellfun(@(a,b)( [colvec(a);colvec(b)] ) , plA.series_quality.estimated_labs, plB.series_quality.estimated_labs, 'UniformOutput', false);
                        payload.series_quality.AnnNames = plA.series_quality.AnnNames;
                    end

                end

            end

        end

        %% property restriction functions

    end

    methods ( Access = private )

    end

end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function pass = test_stringConstructor(~)
%TEST_STRINGCONSTRUCTOR    Test that we can pass strings to the CHEBOP
%constructor

%% Linear problem, no x dependence
% String syntax
L1 = chebop('u``+u',[0 5]);
L1.lbc = 1;
L1.rbc = 2;
u1 = L1\0;

% Anonymous function syntax
L2 = chebop(@(u) diff(u,2) + u, [0 5]);
L2.lbc = 1;
L2.rbc = 2;
u2 = L2\0;

% These should be identical!
pass(1) = norm(u1-u2) == 0;

%% Linear problem, x dependence
% String syntax
L1 = chebop('u``+x*u',[0 5]);
L1.lbc = 1;
L1.rbc = 2;
u1 = L1\0;

% Anonymous function syntax
L2 = chebop(@(x,u) diff(u,2) + x.*u, [0 5]);
L2.lbc = 1;
L2.rbc = 2;
u2 = L2\0;

% These should be identical!
pass(2) = norm(u1-u2) == 0;

%% Nonlinear problem, no x dependence
% String syntax
tic
N1 = chebop('u``+sin(u)',[0 5]);
N1.lbc = 1;
N1.rbc = 2;
u1 = N1\0;

% Anonymous function syntax
N2 = chebop(@(u) diff(u,2) + sin(u), [0 5]);
N2.lbc = 1;
N2.rbc = 2;
u2 = N2\0;

% These should be identical!
pass(3) = norm(u1-u2) == 0;

%% Linear problem, x dependence
% String syntax
N1 = chebop('u``+x*sin(u)',[0 5]);
N1.lbc = 1;
N1.rbc = 2;
u1 = N1\0;

% Anonymous function syntax
N2 = chebop(@(x,u) diff(u,2) + x.*sin(u), [0 5]);
N2.lbc = 1;
N2.rbc = 2;
u2 = N2\0;

% These should be identical!
pass(4) = norm(u1-u2) == 0;

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun