text:
classdef TreeClassifier < Algorithm

    properties (Access = public)
        shouldTrain = false;
        maxNumSplits = 100;
    end

    properties (Access = private)
        classifier;
    end

    methods (Access = public)

        function obj = TreeClassifier(maxNumSplits)
            if nargin > 0
                obj.maxNumSplits = maxNumSplits;
            end
            obj.name = 'Tree';
            obj.inputPort = DataType.kTable;
            obj.outputPort = DataType.kTable;
        end

        function dataOut = compute(obj,data)
            if obj.shouldTrain
                obj.train(data);
                dataOut = [];
            else
                dataOut = obj.test(data);
            end
        end

        function train(obj,table)
            obj.classifier = fitctree(...
                table.features, ...
                table.label, ...
                'SplitCriterion', 'gdi', ...
                'MaxNumSplits', obj.maxNumSplits, ...
                'Surrogate', 'off');
        end

        function metrics = computeMetrics(obj,table)
            flops = timeit(@()obj.test(table)) / Constants.kReferenceComputingTime;
            memory = Helper.ComputeObjectSize(obj.classifier);
            outputSize = table.height * Constants.kClassificationResultBytes;
            metrics = Metric(flops,memory,outputSize);
        end

        function labels = test(obj,table)
            labels = predict(obj.classifier,table.features);
        end

        function str = toString(obj)
            str = sprintf('%s_%d',obj.name,obj.maxNumSplits);
        end

        function editableProperties = getEditableProperties(obj)
            editableProperties = Property('maxNumSplits',obj.maxNumSplits,1,100,PropertyType.kNumber);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function  [fx,dfdx,dfdp] = f_Qlearning(x,P,u,in)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [fx,dfdx,dfdP] = f_Qlearning(x,P,u,in)
% Reinforcement-learning evolution function for a n-armed bandit task
%
% An RL agent learns by trial and error. A bandit task is such that, after
% each action, the agent receives a feedback (reward if positive,
% punishment if negative). The RL agent updates its action values as
% follows:
% V(chosen action) = V(chosen action) + alpha*(feedback-V(chosen action))
% V(unchosen action) = V(unchosen action)
% IN:
%	- x: action values (n x 1)
%	- P: learning rate (will be sigmoid transformed)
%	- u: (1) previous action
%        (2) feedback received for previous action
%	- in: [useless]
% OUT:
%   - fx: updated action values
%   - dfdx/dfdP: gradients for VBA inversion

% /////////////////////////////////////////////////////////////////////////

% Get parameter values
% =========================================================================
% Some of the model parameters can only take values within a certain range.
% However, the VBA toolbox can only perform estimation of unbounded values.
% To solve this, we transform the values passed by VBA to map them on the
% acceptable range for each parameter.

% learning rate
alpha = VBA_sigmoid(P); % [-Inf,Inf] -> [0 1]

% Apply delta-rule to update action values
% =========================================================================

% get experimental conditions
prevActionIdx = u(1)+1; % action 0 is first index
feedback = u(2);

% start with previous values
fx = x;

% udpdate previous action value
delta = feedback - x(prevActionIdx);
fx(prevActionIdx) = x(prevActionIdx) + alpha*delta; % update chosen value

% Compute evolution function's gradient
% =========================================================================
% This is not necessary, as the toolbox will approximate those gradients if
% needed. However, providing the analytical gradient can higly speed-up the
% inversion.

n = numel(x);

% derivative w.r.t hidden state
% -------------------------------------------------------------------------
dfdx = eye(n);
dfdx(prevActionIdx,prevActionIdx) = 1 - alpha;

% derivative w.r.t parameters
% -------------------------------------------------------------------------
dfdp = zeros(1,n);
dfdp(prevActionIdx) = alpha*(1-alpha)*delta;

author: MBB-team

repo: VBA-toolbox


text:
function FLAM = fflaminar(RE)
% FFLAMINAR skin friction factor
%  FFLAMINAR(RE) returns the skin friction loss coefficient
%  for laminar flows in smooth pipes (Reynnolds number < 2300)
%  CALLED FUNCTION: none
%   Required Inputs are:
%    RE  - Reynolds number (-)
% ---------------------------------------------------------------
% The MATLAB function was created by Tibor Balint, December 1998
% TBoreal Research Corporation, Toronto, Ont. Canada
% (tibor@netcom.ca) and also, University of Warwick, UK
% ---------------------------------------------------------------

format long g;                  % set the format of the calculations

if (RE>2300)
   error('The Reynolds number for laminar pipe flows must be less than 2300')
end

%calculate the skin friction factor for laminer flow in pipes
FLAM=64/RE;

return                          %end of the function
% -------------- end of the function ----------------------

author: Sable

repo: mcbench-benchmarks


text:
function v = vec( x )

% VEC   Vectorize.
%    VEC(X), where X is a vector, matrix, or N-D array, returns a column vector
%    containing all of the elements of X; i.e., VEC(X)=X(:).

v = reshape( x, numel( x ), 1 );

% 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 [A_K_I,A_K_F,A_F_I,A_Rfl_F,A_Rfr_F,A_down_up] = Rotation_Matrices(phi,chi,delta)
%% Rotation from rear axle fixed coordinate system into vehicle fixed coordinate system
% Input parameters:
% phi       [rad]    Vehicle Rotation Angles [Phi Theta Psi] [3x1]
% chi       [rad]    Road Cardan Angles [chi_x chi_y Psi] [3x1]
% delta     [rad]    Steering Angles [deltaFL deltaFR 0 0] [1x4]
% Output parameters:
% A_K_I     [---]    Transformation Matrix A_K_I [3x3]
% A_K_F     [---]    Transformation Matrix A_K_F [3x3]
% A_F_I     [---]    Transformation Matrix A_F_I [3x3]
% A_Rfl_F   [---]    Transformation Matrix A_Rfl_F [3x3]
% A_Rfr_F   [---]    Transformation Matrix A_Rfr_F [3x3]
% A_down_up [---]    Transformation Matrix W := A_down_up [3x3]

%% Pre-Computation of trigonometric values
cphi=cos(phi);
sphi=sin(phi);
cchi=cos(chi);
schi=sin(chi);

%% Definition of rotation matrices

% rotation matrix from z-up coordinate system to z-down coordinate system
A_down_up = reshape([1 0 0,0 -1 0,0 0 -1],[3,3]);

% rotation matrix from inertial coordinate system to vehicle fixed coordinate system
A_K_I     = reshape([cphi(3)*cphi(2),-cphi(1)*sphi(3)+cphi(3)*sphi(2)*sphi(1),sphi(3)*sphi(1)+cphi(3)*cphi(1)*sphi(2),cphi(2)*sphi(3),cphi(3)*cphi(1)+sphi(3)*sphi(2)*sphi(1),-cphi(3)*sphi(1)+cphi(1)*sphi(3)*sphi(2),-sphi(2),cphi(2)*sphi(1),cphi(2)*cphi(1)],[3,3]);

% rotation matrix from rear axle fixed coordinate system to vehicle fixed coordinate system without road inclination
% A_K_F     = reshape([cphi(2),sphi(2)*sphi(1),cphi(1)*sphi(2),0,cphi(1),-sphi(1),-sphi(2),cphi(2)*sphi(1),cphi(2)*cphi(1)],[3,3]);

% rotation matrix from rear axle fixed coordinate system to vehicle fixed coordinate system over inertial coordinate system with road inclination
A_K_F     = reshape([cos(phi(2)-chi(2)),sin(phi(2)-chi(2))*sphi(1),sin(phi(2)-chi(2))*cphi(1),-sin(-chi(2)+phi(2))*schi(1),cchi(1)*cphi(1)+cphi(2)*cchi(2)*schi(1)*sphi(1)+sphi(2)*schi(1)*schi(2)*sphi(1),-cchi(1)*sphi(1)+cphi(2)*cchi(2)*cphi(1)*schi(1)+cphi(1)*sphi(2)*schi(1)*schi(2),-sin(-chi(2)+phi(2))*cchi(1),-cphi(1)*schi(1)+cphi(2)*cchi(1)*cchi(2)*sphi(1)+cchi(1)*sphi(2)*schi(2)*sphi(1),schi(1)*sphi(1)+cphi(2)*cchi(1)*cchi(2)*cphi(1)+cchi(1)*cphi(1)*sphi(2)*schi(2)],[3,3]);

% rotation matrix from inertial coordinate system to rear axle fixed coordinate system
A_F_I     = reshape([cphi(3)*cchi(2),-cchi(1)*sphi(3)+cphi(3)*schi(1)*schi(2),sphi(3)*schi(1)+cphi(3)*cchi(1)*schi(2),cchi(2)*sphi(3),cphi(3)*cchi(1)+sphi(3)*schi(1)*schi(2),-cphi(3)*schi(1)+cchi(1)*sphi(3)*schi(2),-schi(2),cchi(2)*schi(1),cchi(1)*cchi(2)],[3,3]);

% rotation matrix from rear axle fixed coordinate system to wheel fixed coordinate systems of front axle
A_Rfl_F   = reshape([cos(delta(1)),-sin(delta(1)),0,sin(delta(1)),cos(delta(1)),0,0,0,1],[3,3]);
A_Rfr_F   = reshape([cos(delta(2)),-sin(delta(2)),0,sin(delta(2)),cos(delta(2)),0,0,0,1],[3,3]);

end

author: TUMFTM


text:
function display(panel)
% Display function for BstPanel object.

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

% If displaying a template : do not display anything
if ((size(panel,1)==0) || (size(panel,2)==0))
    disp('Template Brainstorm panel...');
    return
end

for iPanel = 1:length(panel)
    % Display title and simple properties
    disp('Brainstorm window panel:');
    disp(sprintf('           name: %s', panel(iPanel).name));
    disp(sprintf('        jHandle: %s', char(panel(iPanel).jHandle.getClass())));
    switch (panel(iPanel).container.type)
        case ''
            disp(sprintf('      container: N/A'));
        case 'JavaWindow'
            disp(sprintf('      container: JavaWindow (%s)', char(panel(iPanel).container.handle{1})));
        case 'MatlabFigure'
            disp(sprintf('      container: MatlabFigure (%d)', panel(iPanel).container.handle{1}));
        case 'BrainstormTab'
            disp(sprintf('      container: BrainstormTab (''%s'')', panel(iPanel).container.handle{1}));
        case 'BrainstormPanel'
            disp(sprintf('      container: BrainstormPanel (''%s'')', panel(iPanel).container.handle{1}));
    end
    % Display sControls structure
    sControlsStr = '      sControls: ';
    sControlsFields = fieldnames(panel(iPanel).sControls);
    if isempty(sControlsFields)
        % No fields
        sControlsStr = [sControlsStr 'No controls registered'];
    else
        % Display all the fields
        for i=1:length(sControlsFields)
            if isjava(panel(iPanel).sControls.(sControlsFields{i}))
                controlDesc = class(panel(iPanel).sControls.(sControlsFields{i}));
            else
                controlDesc = 'User value';
            end
            if (i == 1)
                sControlsStr = sprintf('%s%s (%s)', sControlsStr, sControlsFields{i}, controlDesc);
            else
                sControlsStr = sprintf('%s\n                 %s (%s)', sControlsStr, sControlsFields{i}, controlDesc);
            end
        end
    end
    disp(sControlsStr);
    disp(' ');
end

author: brainstorm-tools


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DRAWROBOT3D(ROBOT, Q)	3D drawing with DH reference systems of the robot at the
% current joint coordinates Q. The robot parameters are stored in the
% variable ROBOT. Q = [q1 q2 ... qN] is a vector of joint values.
%
% If the robot possesses 3D graphics, the function DRAW_LINK is called to represent
% the link in 3D. Otherwise a line connecting each consecutive reference
% system is plotted.
%
% See also DENAVIT, DIRECTKINEMATIC, DRAW_LINK.
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche.
%   email: arturo.gil@umh.es date:   05/02/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%ligth;

%if q is not specified, then assume zero pose
if ~exist('q', 'var')
    q=zeros(1,robot.DOF)
end

h=figure(configuration.figure.robot);

%get adjusted view
[az,el] = view;

%retrieve the current zoom
zz=get(gca,'CameraViewAngle');

%delete figure
if ~exist('noclear')
    clf(h), grid on, hold on
end

xlabel('X (m)'), ylabel('Y (m)'), zlabel('Z (m)')

%apply the stored view
view(az,el);
%apply the zoom
set(gca,'CameraViewAngle', zz);

origin = [0 0 0];
destination = [];
color = ['r' 'g' 'b' 'c' 'm' 'y' 'k'];

%in case of a parallel robot, call drawrobot3d for each serial link
if isfield(robot, 'parallel')
    drawrobot3d_par(robot, q);
    return;
end

if isfield(robot, 'skip_graphics')
    return;
end

V=[];

T = robot.T0;
for i=1:robot.DOF+1

    if robot.graphical.has_graphics
        draw_link(robot, i, T);
    else
        draw_link(robot, i, T);

        destination=T(1:3,4);
        col = randperm(6);
        line([origin(1) destination(1)],[origin(2) destination(2)],[origin(3) destination(3)], 'Color',  color(col(1)) ,'LineWidth', 3 );
        origin = destination;
    end
    %draw D-H axes
    if robot.graphical.draw_axes
        draw_axes(T, sprintf('X_%d',i-1), sprintf('Y_%d',i-1), sprintf('Z_%d',i-1), robot.graphical.axes_scale);
    end
    %obtain transformation to draw the next link
    if i<= robot.DOF
       T=T*dh(robot, q, i);
    end
end
%save for later use
T06robot=T;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DRAW A TOOL AT THE ROBOT'S END
% if there is a tool attached, draw it!
% repeat the prior commands
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isfield(robot, 'tool')

    % Activate tool. Tools must be defined in a binary way.
    % This means that they are open or closed depending
    % on an joint
    if robot.tool.tool_open
        q_tool=ones(1,robot.tool.DOF);
    else
        q_tool=zeros(1,robot.tool.DOF);
    end
    robot.tool.T0=T;
    for i=1:robot.tool.DOF+1
        if robot.tool.graphical.has_graphics
            draw_link(robot.tool, i, T);
        else
            draw_link(robot.tool, i, T);

            destination=T(1:3,4);
            col = randperm(6);
            line([origin(1) destination(1)],[origin(2) destination(2)],[origin(3) destination(3)], 'Color',  color(col(1)) ,'LineWidth', 3 );
            origin = destination;
        end
%         if robot.tool.graphical.draw_axes
%             draw_axes(T, sprintf('X_{tool%d}',i-1), sprintf('Y_{tool%d}',i-1), sprintf('Z_{tool%d}',i-1));
%         end
        %obtain transformation to draw the next link
        if (i)<= robot.tool.DOF
            T=T*dh(robot.tool, q_tool, i);
        end
    end
    %finally draw the Tool Center Point (TCP)
    if robot.tool.graphical.draw_axes
        draw_axes(T06robot*robot.tool.TCP, 'X_{toolTCP}', 'Y_{toolTCP}', 'Z_{toolTCP}');
    end

    %save for later use
    T07=T06robot*robot.tool.TCP;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  draw equipment such as tables, conveyor velts... etc
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isfield(robot, 'equipment')
    for i=1:length(robot.equipment),
        if robot.equipment{i}.graphical.has_graphics
            draw_link(robot.equipment{i}, 1, robot.equipment{i}.T0);
        end
        if robot.equipment{i}.graphical.draw_axes
            draw_axes(robot.equipment{i}.T0, sprintf('X_{env%d}',0), sprintf('Y_{env%d}',0), sprintf('Z_{env%d}',0));
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   DRAW A PIECE GRABBED
% this draws a piece grabbed at the robots end, or at the last known position
% the variable robot.tool_activated=1 asumes that the piece has been grabbed by the end tool
% robot.tool_activated==0 maintains the last known location for the pieced
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isfield(robot, 'piece')
    for i=1:size(robot.piece,2)
        pieza_dibujar = robot.piece{i};
        if isfield(pieza_dibujar, 'graphical')
            if pieza_dibujar.graphical.has_graphics
                if isfield(robot, 'tool')
                    if robot.tool.piece_gripped==1
                        %update the last known position and orientation of the robot
                        %Trel is computed when the simulation_grip_piece function is
                        %executed. Trel is the relative position and orientation of
                        %the tool and the piece that assures that the piece is
                        %picked at a constant and visually effective orientation.
                        draw_link(pieza_dibujar, 1, T07*(robot.tool.Trel));
                    end
                else
                    draw_link(pieza_dibujar, 1, pieza_dibujar.T0);
                end
            end
        end
        if pieza_dibujar.graphical.draw_axes
            draw_axes(pieza_dibujar.T0, sprintf('X_{piece%d}',0), sprintf('Y_{piece%d}',0), sprintf('Z_{piece%d}',0));
        end
    end
end

end

repo: ARTE


text:
% adjustcylinder() - Adjust 3d object coordinates to match a pair of points
%
% Usage:
%   >> [x y z] = adjustcylinder( x, y, z, pos1, pos2);
%
% Inputs:
%  x,y,z      - 3-D point coordinates
%  pos1       - position of first point [x y z]
%  pos2       - position of second point [x y z]
%
% Outputs:
%  x,y,z      - updated 3-D point coordinates
%
% Author: Arnaud Delorme, CNL / Salk Institute, 30 Mai 2003

% Copyright (C) 2003 Arnaud Delorme
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [x, y, z] = adjustcylinder2( h, pos1, pos2);

    % figure; plot3(x(2,:),y(2,:),z(2,:)); [ x(2,:)' y(2,:)' z(2,:)']

    % stretch z coordinates to match for vector length
    % ------------------------------------------------
    dist = sqrt(sum((pos1-pos2).^2));
    z = get(h, 'zdata');
    zrange = max(z(:)) - min(z(:));
    set(h, 'zdata', get(h, 'zdata') /zrange*dist);

    % rotate in 3-D to match vector angle [0 0 1] -> vector angle)
    % only have to rotate in the x-z and y-z plane
    % --------------------------------------------
    vectrot = [ pos2(1)-pos1(1) pos2(2)-pos1(2) pos2(3)-pos1(3)];
    [thvect phivect] = cart2sph( vectrot(1), vectrot(2), vectrot(3) );

    rotate(h, [0 0 1], thvect/pi*180, [0 0 0]);
    rotate(h, [thvect+pi/2 0]/pi*180, (pi/2-phivect)/pi*180, [0 0 0]);

    x = get(h, 'xdata') + pos1(1);
    y = get(h, 'ydata') + pos1(2);
    z = get(h, 'zdata') + pos1(3);

    set(h, 'xdata', x);
    set(h, 'ydata', y);
    set(h, 'zdata', z);
    return;

author: PatternRecognition

repo: OpenBMI


text:
function szmM = calcSZM(quantizedM, nL, szmType)
% function szmM = calcSZM(quantizedM, nL, szmType)
%
% This function calculates the Size-Zone matrix for the passed quantized
% image.
%
% INPUTS:
%       quantizedM: quantized 3d matrix obtained, for example, by
%       imquantize_cerr.m
%       nL: Number of gray levels.
%       szmType: flag, 1 or 2.
%                   1: 3D zones
%                   2: 2D zones
% OUTPUT:
%       szmM: size-zone matrix of size (nL x L)
%
% EXAMPLE:
%
% numRows = 10;
% numCols = 10;
% numSlcs = 1;
%
% % number of gray levels
% nL = 3;
%
% % create an image with random numbers
% imgM = randi(nL,numRows,numCols,numSlcs);
%
% % set option to add run lengths from all directions
% szmType = 1;
%
% % call the rlm calculator
% szmType = calcSZM(imgM, nL, szmType);
%
%
% APA, 03/30/2017

if szmType == 1
    numNeighbors = 26;
else
    numNeighbors = 8;
end

szmM = sparse(nL,numel(quantizedM));
maxSiz = 0;
for level = 1:nL
    connM = bwlabeln(quantizedM==level, numNeighbors);
    if any(connM(:))
        regiosSizV = accumarray(connM(connM > 0),1);
        if ~isempty(regiosSizV)
            maxSiz = max(maxSiz,max(regiosSizV));
        end
        szmM(level,:) = accumarray(regiosSizV,1,[size(szmM,2) 1])';
    end
end
szmM = szmM(:,1:maxSiz);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function cmap= cmap_whitered(m, varargin)
%CMAP_WHITERED - Colormap going from white to red
%
%Synopsis:
% MAP= cmap_whitered(M, <OPT>)
%
%Input:
% M: Size of the colormap (number of entries). Default value: Same size
%    as current colormap
% OPT: Struct or property/value list of optinal properties:
%  .MinSat - minimal saturation (HSV model) from which fading is started
%  .MinVal - minimal value (in HSV model) from which fading is started
%  .MaxVal - maximal value (in HSV model) to which fading is performed
%
%Output:
% MAP: A colormap matrix of size [M 3]
%
%Example:
% clf;
% colormap(cmap_whitered(15));
% imagesc(toeplitz(1:15)); colorbar;
%
%See also COLORMAP, HSV2RGB, CMAP_HSVFADE

props= {'MinSat',  0.4,  'DOUBLE[1]'
        'MinVal',  0,    'DOUBLE[1]'
        'MaxVal',  0.8,  'DOUBLE[1]'};

if nargin<1 | isempty(m),
  m= size(get(gcf,'colormap'),1);
end

opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

m1= floor(m/2.5);
m2= m-m1-1;
map1= cmap_hsvFade(m1+1, 0, 1, [opt.MinSat 1]);
map2= cmap_hsvFade(m2+1, 0, [opt.MaxVal opt.MinVal], 1);

cmap= flipud([map1; map2(2:end,:)]);

author: bbci


text:
function [ vX, vV, mX ] = ProjectedGdFista( vX, vV, hG, hP, numIterations, stepSize )
% ----------------------------------------------------------------------------------------------- %
% [ vX, vV, mX ] = ProjectedGdFista( vX, vV, hG, hP, numIterations, stepSize )
% Solves an objective function by the Projected Gradient Descent /
% Projected Gradient Method (PGM) with FISTA Acceleration.
% Input:
%   - vX                -   Input Vector.
%                           The starting point for the gradient iterations.
%                           Structure: Vector (numElements x 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vV                -   Acceleration Vector.
%                           The buffer to hold the accelerated direction
%                           vector.
%                           Structure: Vector (numElements x 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - hG                -   Gradient Function Handler.
%                           A function handler which accepts a vector of
%                           size (numElements x 1) and returns the gradient
%                           in the form of a vector (numElements x 1).
%                           Structure: Function Handler.
%                           Type: NA.
%                           Range: NA.
%   - hP                -   Gradient Function Handler.
%                           A function handler which accepts a vector of
%                           size (numElements x 1) and returns the gradient
%                           in the form of a vector (numElements x 1).
%                           Structure: Function Handler.
%                           Type: NA.
%                           Range: NA.
%   - numIterations     -   Number of Iterations.
%                           Sets the number of iterations of the gradient
%                           descent.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - stepSize          -   Step Size.
%                           Sets the step size of the Gradient step.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
% Output:
%   - vX                -   Output Vector.
%                           The end point for the gradient iterations.
%                           Structure: Vector (numElements x 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vV                -   Acceleration Vector.
%                           The buffer to hold the accelerated direction
%                           vector of last iteration. May used by a wrapper
%                           for warm start.
%                           Structure: Vector (numElements x 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - mX                -   The Path Matrix.
%                           The points of the gradient iterations.
%                           Structure: Matrix (numElements x numIterations).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  Geoff Gordon & Ryan Tibshirani - Accelerated First Order Methods (https://www.cs.cmu.edu/~ggordon/10725-F12/slides/09-acceleration.pdf).
% Remarks:
%   1.  B
% TODO:
%   1.  C
% Release Notes:
%   -   1.0.000     26/12/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

mX = zeros(size(vX, 1), numIterations);

vG          = zeros(size(vX, 1), 1);
mX(:, 1)    = vX;
mX(:, 2)    = vX - stepSize * hG(vX);

for ii = 3:numIterations
    vV(:) = mX(:, ii - 1) + (((ii - 2) / (ii + 1)) * (mX(:, ii - 1) - mX(:, ii - 2)));
    vG(:) = hG(vV); %<! The gradient
    vX(:) = hP(vV - (stepSize * vG)); %<! Projection step

    mX(:, ii) = vX;
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function [weight] = anticor_expert(data, weight_o, w)
% This program generates portfolio for a specified parameter setting.
% Anticor expert
% while the original version is based on loop, which is slow in matlab,
% I replace the code using matrix operations
% Most previous experiments run based on the commented version, at the end
% of the file.
%
% function [weight] = anticor_expert(data, weight_o, w)
%
% weight: experts portfolio, used for next rebalance/combination
%
% data: market sequence vectors
% weight_o: last portfolio, also can be last price relative adjusted.
% w: window size for the experts
%
% Example: [weight] = anticor_expert(data, weight_o, 30)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[T, N]=size(data);
weight = weight_o;

% 1. Return the current portfolio if t < 2w
if (T >= 2*w),
    % 2. Compute LX1, LX2
    LX1 = log(data((T-2*w+1):(T-w),:));
    LX2 = log(data((T-w+1):T,:));

    % 2. Compute averages of LX1 and LX2
    mu1 = mean(LX1);
    mu2 = mean(LX2);

    % M_cor = corr(LX1, LX2);
    M_cov = zeros(N, N);
    M_cor = zeros(N, N);

    % I replace the following loop using matrix operations
    n_LX1 = LX1(:, :) - repmat(mu1, [w, 1]);
    n_LX2 = LX2(:, :) - repmat(mu2, [w, 1]);

    Std1 = diag(n_LX1'*n_LX1)/(w-1);
    Std2 = diag(n_LX2'*n_LX2)/(w-1);
    Std12 = Std1*Std2';
    M_cov = n_LX1'*n_LX2/(w-1);

    M_cor(Std12==0) = 0;
    M_cor(Std12~=0) = M_cov(Std12~=0)./sqrt(Std12(Std12~=0));

    % claim varialbe, matrixization
    claim = zeros(N);
    w_mu2 = repmat(mu2', [1, N]);
    w_mu1 = repmat(mu2, [N, 1]);
    s_12 = (w_mu2 >= w_mu1) & (M_cor > 0);
    claim(s_12) = claim(s_12) + M_cor(s_12);

    diag_M_cor = diag(M_cor);
    cor1 = max(0, repmat(-diag_M_cor, [1, N]));
    cor2 = max(0, repmat(-diag_M_cor', [N, 1]));
    claim(s_12) = claim(s_12) + cor1(s_12) + cor2(s_12);

    transfer = zeros(N);
    sum_claim = repmat(sum(claim, 2), [1, N]);
    s_1 = abs(sum_claim) > 0;
    w_weight_o = repmat(weight_o, [1, N]);
    transfer(s_1) = w_weight_o(s_1).*claim(s_1)./sum_claim(s_1);

    % Use matrix operations, replacing the following loop
    transfer_ij = transfer'-transfer;
    weight = weight - sum(transfer_ij)';
end
end

%     for i=1:N,
%         for j = 1:N,
%             M_cov(i, j) = (LX1(:, i)-mu1(:, i))'*(LX2(:, j) - mu2(:, j))/(w-1);
%
%             % Std1 = Std(LX1(:, i)
%             Std1 = (LX1(:, i) - mu1(:, i))'*(LX1(:, i) - mu1(:, i))/(w-1);
%             Std2 = (LX2(:, j) - mu2(:, j))'*(LX2(:, j) - mu2(:, j))/(w-1);
%
%             if (Std12(i, j) == 0),
%                 M_cor(i, j) = 0;
%             else
%                 M_cor(i, j) = M_cov(i, j) / sqrt(Std1(i, 1)*Std2(j, 1));
%             end
%         end
%     end

%
%     for i = 1:N,
%         for j = 1:N,
%             if ((mu2(i) >= mu2(j)) && (M_cor(i, j) > 0)),
%                 claim(i, j) = claim(i, j) + M_cor(i, j);
% %                 claim(i, j) = claim(i, j) + max(0, -M_cor(i, i));
% %                 claim(i, j) = claim(i, j) + max(0, -M_cor(j, j));
%             end
%         end
%     end
%     for i = 1:N,
%         sum_claim_i = sum(claim(i, :));
%         if (abs(sum_claim_i) > 0)
%             for j = 1:N,
%                 transfer(i, j) = weight_o(i)*claim(i, j)/sum_claim_i;
%             end
%         end
%     end
%     for i = 1:N,
%         for j = 1:N,
%             weight(i) = weight(i) - transfer(i, j) + transfer(j, i);
%         end
%     end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function plotIterativeValdata(valdata, fignum)

cthresh = [valdata(1, :).cthresh];
err = zeros(size(cthresh));
nregions = zeros(size(cthresh));
for f = 1:size(valdata,1)
    for k = 1:numel(cthresh)
        err(k) = err(k) + (1-valdata(f,k).conservation) / size(valdata,1);
        nregions(k) = nregions(k) + valdata(f,k).nregions / size(valdata,1);
    end
end
figure(fignum), subplot(1,2,1), hold off, plot(cthresh, err), title('thresh vs. err')
figure(fignum), subplot(1,2,2), plot(err, nregions), title('conf: err vs. nregions')
drawnow;

author: Cloud-CV

repo: object-proposals


text:
function F=dwt_sr_fuse(I1,I2,zt,D,overlap,epsilon)
%    DWT-SR
%    Input:
%    I1 - input image A
%    I2 - input image B
%    zt - maximum decomposition level
%    D  - Dictionary for sparse representation
%    overlap - the overlapped pixels between two neighbor patches
%    epsilon - sparse reconstuction error
%    Output:
%    F  - fused image
%
%    The code is edited by Yu Liu, 01-09-2014.

%-------------------------------------------------------------------------%
%                                DWT
%-------------------------------------------------------------------------%
I1=double(I1);
I2=double(I2);
tempA=I1;
tempB=I2;

X=cell(zt,4);
Y=cell(zt,4);
Z=cell(zt,4);
for i=1:zt
    [X{i,1},X{i,2},X{i,3},X{i,4}]=dwt2(tempA,'db1','mode','per');
    tempA=X{i,1};
    [Y{i,1},Y{i,2},Y{i,3},Y{i,4}]=dwt2(tempB,'db1','mode','per');
    tempB=Y{i,1};
end

%-------------------------------------------------------------------------%
%                               low-pass fusion
%-------------------------------------------------------------------------%
Z{zt,1}=sparse_fusion(X{zt,1},Y{zt,1},D,overlap,epsilon);

%-------------------------------------------------------------------------%
%                               high-pass fusion
%-------------------------------------------------------------------------%
for i=zt:-1:1
    for j=2:4
        Z{i,j}=selc(X{i,j},Y{i,j},3);
    end
end

%-------------------------------------------------------------------------%
%                               IDWT
%-------------------------------------------------------------------------%
for i=zt:-1:1
    if i>1
        Z{i-1,1}=idwt2(Z{i,1},Z{i,2},Z{i,3},Z{i,4},'db1','mode','per');
    else
        F=idwt2(Z{i,1},Z{i,2},Z{i,3},Z{i,4},'db1','mode','per');
    end
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
%WAVELET  1D Wavelet transform with optional singificance testing
%
%   [WAVE,PERIOD,SCALE,COI] = wavelet(Y,DT,PAD,DJ,S0,J1,MOTHER,PARAM)
%
%   Computes the wavelet transform of the vector Y (length N),
%   with sampling rate DT.
%
%   By default, the Morlet wavelet (k0=6) is used.
%   The wavelet basis is normalized to have total energy=1 at all scales.
%
%
% INPUTS:
%
%    Y = the time series of length N.
%    DT = amount of time between each Y value, i.e. the sampling time.
%
% OUTPUTS:
%
%    WAVE is the WAVELET transform of Y. This is a complex array
%    of dimensions (N,J1+1). FLOAT(WAVE) gives the WAVELET amplitude,
%    ATAN(IMAGINARY(WAVE),FLOAT(WAVE) gives the WAVELET phase.
%    The WAVELET power spectrum is ABS(WAVE)^2.
%    Its units are sigma^2 (the time series variance).
%
%
% OPTIONAL INPUTS:
%
% *** Note *** setting any of the following to -1 will cause the default
%               value to be used.
%
%    PAD = if set to 1 (default is 0), pad time series with enough zeroes to get
%         N up to the next higher power of 2. This prevents wraparound
%         from the end of the time series to the beginning, and also
%         speeds up the FFT's used to do the wavelet transform.
%         This will not eliminate all edge effects (see COI below).
%
%    DJ = the spacing between discrete scales. Default is 0.25.
%         A smaller # will give better scale resolution, but be slower to plot.
%
%    S0 = the smallest scale of the wavelet.  Default is 2*DT.
%
%    J1 = the # of scales minus one. Scales range from S0 up to S0*2^(J1*DJ),
%        to give a total of (J1+1) scales. Default is J1 = (LOG2(N DT/S0))/DJ.
%
%    MOTHER = the mother wavelet function.
%             The choices are 'MORLET', 'PAUL', or 'DOG'
%
%    PARAM = the mother wavelet parameter.
%            For 'MORLET' this is k0 (wavenumber), default is 6.
%            For 'PAUL' this is m (order), default is 4.
%            For 'DOG' this is m (m-th derivative), default is 2.
%
%
% OPTIONAL OUTPUTS:
%
%    PERIOD = the vector of "Fourier" periods (in time units) that corresponds
%           to the SCALEs.
%
%    SCALE = the vector of scale indices, given by S0*2^(j*DJ), j=0...J1
%            where J1+1 is the total # of scales.
%
%    COI = if specified, then return the Cone-of-Influence, which is a vector
%        of N points that contains the maximum period of useful information
%        at that particular time.
%        Periods greater than this are subject to edge effects.
%        This can be used to plot COI lines on a contour plot by doing:
%
%              contour(time,log(period),log(power))
%              plot(time,log(coi),'k')
%
%----------------------------------------------------------------------------
%   Copyright (C) 1995-2004, Christopher Torrence and Gilbert P. Compo
%
%   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.
%
% Notice: Please acknowledge the use of the above software in any publications:
%    ``Wavelet software was provided by C. Torrence and G. Compo,
%      and is available at URL: http://paos.colorado.edu/research/wavelets/''.
%
% Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
%            Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, 61-78.
%
% Please send a copy of such publications to either C. Torrence or G. Compo:
%  Dr. Christopher Torrence               Dr. Gilbert P. Compo
%  Research Systems, Inc.                 Climate Diagnostics Center
%  4990 Pearl East Circle                 325 Broadway R/CDC1
%  Boulder, CO 80301, USA                 Boulder, CO 80305-3328, USA
%  E-mail: chris[AT]rsinc[DOT]com         E-mail: compo[AT]colorado[DOT]edu
%----------------------------------------------------------------------------
function [wave,period,scale,coi] = ...
	wavelet(Y,dt,pad,dj,s0,J1,mother,param);

if (nargin < 8), param = -1;, end
if (nargin < 7), mother = -1;, end
if (nargin < 6), J1 = -1;, end
if (nargin < 5), s0 = -1;, end
if (nargin < 4), dj = -1;, end
if (nargin < 3), pad = 0;, end
if (nargin < 2)
	error('Must input a vector Y and sampling time DT')
end

n1 = length(Y);

if (s0 == -1), s0=2*dt;, end
if (dj == -1), dj = 1./4.;, end
if (J1 == -1), J1=fix((log(n1*dt/s0)/log(2))/dj);, end
if (mother == -1), mother = 'MORLET';, end

%....construct time series to analyze, pad if necessary
x(1:n1) = Y - mean(Y);
if (pad == 1)
	base2 = fix(log(n1)/log(2) + 0.4999);   % power of 2 nearest to N
	x = [x,zeros(1,2^(base2+1)-n1)];
end
n = length(x);

%....construct wavenumber array used in transform [Eqn(5)]
k = [1:fix(n/2)];
k = k.*((2.*pi)/(n*dt));
k = [0., k, -k(fix((n-1)/2):-1:1)];

%....compute FFT of the (padded) time series
f = fft(x);    % [Eqn(3)]

%....construct SCALE array & empty PERIOD & WAVE arrays
scale = s0*2.^((0:J1)*dj);
period = scale;
wave = zeros(J1+1,n);  % define the wavelet array
wave = wave + i*wave;  % make it complex

% loop through all scales and compute transform
for a1 = 1:J1+1
	[daughter,fourier_factor,coi,dofmin]=wave_bases(mother,k,scale(a1),param);
	wave(a1,:) = ifft(f.*daughter);  % wavelet transform[Eqn(4)]
end

period = fourier_factor*scale;
coi = coi*dt*[1E-5,1:((n1+1)/2-1),fliplr((1:(n1/2-1))),1E-5];  % COI [Sec.3g]
wave = wave(:,1:n1);  % get rid of padding before returning

return

%WAVE_BASES  1D Wavelet functions Morlet, Paul, or DOG
%
%  [DAUGHTER,FOURIER_FACTOR,COI,DOFMIN] = ...
%      wave_bases(MOTHER,K,SCALE,PARAM);
%
%   Computes the wavelet function as a function of Fourier frequency,
%   used for the wavelet transform in Fourier space.
%   (This program is called automatically by WAVELET)
%
% INPUTS:
%
%    MOTHER = a string, equal to 'MORLET' or 'PAUL' or 'DOG'
%    K = a vector, the Fourier frequencies at which to calculate the wavelet
%    SCALE = a number, the wavelet scale
%    PARAM = the nondimensional parameter for the wavelet function
%
% OUTPUTS:
%
%    DAUGHTER = a vector, the wavelet function
%    FOURIER_FACTOR = the ratio of Fourier period to scale
%    COI = a number, the cone-of-influence size at the scale
%    DOFMIN = a number, degrees of freedom for each point in the wavelet power
%             (either 2 for Morlet and Paul, or 1 for the DOG)
%
%----------------------------------------------------------------------------
%   Copyright (C) 1995-1998, Christopher Torrence and Gilbert P. Compo
%   University of Colorado, Program in Atmospheric and Oceanic Sciences.
%   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.
%----------------------------------------------------------------------------
function [daughter,fourier_factor,coi,dofmin] = ...
	wave_bases(mother,k,scale,param);

mother = upper(mother);
n = length(k);

if (strcmp(mother,'MORLET'))  %-----------------------------------  Morlet
	if (param == -1), param = 6.; end
	k0 = param;
	expnt = -(scale.*k - k0).^2/2.*(k > 0.);
	norm = sqrt(scale*k(2))*(pi^(-0.25))*sqrt(n);    % total energy=N   [Eqn(7)]
	daughter = norm*exp(expnt);
	daughter = daughter.*(k > 0.);     % Heaviside step function
	fourier_factor = (4*pi)/(k0 + sqrt(2 + k0^2)); % Scale-->Fourier [Sec.3h]
	coi = fourier_factor/sqrt(2);                  % Cone-of-influence [Sec.3g]
	dofmin = 2;                                    % Degrees of freedom
elseif (strcmp(mother,'PAUL'))  %--------------------------------  Paul
	if (param == -1), param = 4.; end
	m = param;
	expnt = -(scale.*k).*(k > 0.);
	norm = sqrt(scale*k(2))*(2^m/sqrt(m*prod(2:(2*m-1))))*sqrt(n);
	daughter = norm*((scale.*k).^m).*exp(expnt);
	daughter = daughter.*(k > 0.);     % Heaviside step function
	fourier_factor = 4*pi/(2*m+1);
	coi = fourier_factor*sqrt(2);
	dofmin = 2;
elseif (strcmp(mother,'DOG'))  %--------------------------------  DOG
	if (param == -1), param = 2.;, end
	m = param;
	expnt = -(scale.*k).^2 ./ 2.0;
	norm = sqrt(scale*k(2)/gamma(m+0.5))*sqrt(n);
	daughter = -norm*(i^m)*((scale.*k).^m).*exp(expnt);
	fourier_factor = 2*pi*sqrt(2./(2*m+1));
	coi = fourier_factor/sqrt(2);
	dofmin = 1;
else
	error('Mother must be one of MORLET,PAUL,DOG')
end

return

% end of code

author: jramshur

repo: HRVAS

save_path: github-repos/MATLAB/jramshur-HRVAS


text:
%% Generate best fit function
target_vy = 0.0;%[-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8];
% target_vy = [0];
target_vx = [-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4];
% target_vx = [0];
T = 0.4;
subfolder_name = 'library4';
if ~exist(fullfile('local', subfolder_name, 'periodic_fitting'), 'dir')
    mkdir(fullfile('local', subfolder_name, 'periodic_fitting'));
end
joint_names = {    'BasePosX'
    'BasePosY'
    'BasePosZ'
    'BaseRotX'
    'BaseRotY'
    'BaseRotZ'
    'qHRight'
    'qARight'
    'qBRight'
    'fourBarARight'
    'fourBarBRight'
    'qHLeft'
    'qALeft'
    'qBLeft'
    'fourBarALeft'
    'fourBarBLeft'};

N = 16;
N_vx = length(target_vx);
N_vy = length(target_vy);
t_all = cell(N_vx, N_vy);
q_all = cell(N_vx, N_vy);
dq_all = cell(N_vx, N_vy);
px_all = cell(N_vx, N_vy);
py_all = cell(N_vx, N_vy);
vx_all = cell(N_vx, N_vy);
vy_all = cell(N_vx, N_vy);
counter = 1;
for i = 1:N_vx
    vx = target_vx(i);
    for j = 1:N_vy
        vy = target_vy(j);
        data_name = fullfile('local', subfolder_name, sprintf('gait_X%0.1f_Y%.1f.mat', vx, vy));
        param = load(data_name);

        t_all{i,j} = [param.gait(1).tspan]./T;%, param.gait(3).tspan];
        q_all{i,j} = [param.gait(1).states.x];%, param. gait(3).states.x];
        dq_all{i,j} = [param.gait(1).states.dx];%, param. gait(3).states.dx];
        px_all{i,j} = [param.gait(1).states.x(1,:)];%ones(size(t_all{i,j}))*vx;
        py_all{i,j} = [param.gait(1).states.x(2,:)];%ones(size(t_all{i,j}))*vy;
        vx_all{i,j} = [param.gait(1).states.dx(1,:)];%ones(size(t_all{i,j}))*vx;
        vy_all{i,j} = [param.gait(1).states.dx(2,:)];%ones(size(t_all{i,j}))*vy;

    end
end

inpt = [horzcat(t_all{:})
%     horzcat(px_all{:})
%     horzcat(py_all{:})
    horzcat(vx_all{:})
%     horzcat(vx_all{:})
    ];
q_outpt = horzcat(q_all{:});
dq_outpt = horzcat(dq_all{:});
%%
for i = 1:N
    q_outpt_specific = q_outpt(i, :);
    dq_outpt_specific = dq_outpt(i, :);
    opt.neuralFitting(q_outpt_specific, inpt, 5, fullfile('local', subfolder_name, 'periodic_fitting', sprintf('sagittal_library_q%i', i)));
    opt.neuralFitting(dq_outpt_specific, inpt, 5, fullfile('local', subfolder_name, 'periodic_fitting', sprintf('sagittal_library_dq%i', i)));
end
addpath(fullfile('local', subfolder_name, 'periodic_fitting'));
%% Plot features
% joint angle
for n = 1:N
    f = figure(1000+n);
    f.Name = joint_names{n};
    set(f, 'WindowStyle', 'docked');

    ax = axes(f);
    hold(ax);

    scatter3(ax, inpt(1,:),inpt(2,:),q_outpt(n,:), 'r');

    [S, V] = meshgrid(0:0.01:1, -2:0.01:2);
    L = zeros(size(S));
    for i = 1:size(S, 1)
        for j = 1:size(S, 2)
            L(i, j) = feval(sprintf('sagittal_library_q%i', n), ([S(i, j); V(i, j)]));
            %L(i, j) = l(n);
        end
    end

    surface(ax, S, V, L);
end
%% joint velocity
for n = 1:N
    f = figure(2000+n);
    f.Name = joint_names{n};
    set(f, 'WindowStyle', 'docked');

    ax = axes(f);
    hold(ax);

    scatter3(ax, inpt(1,:),inpt(2,:),dq_outpt(n,:), 'r');

    [S, V] = meshgrid(0:0.01:1, -2:0.01:2);
    L = zeros(size(S));
    for i = 1:size(S, 1)
        for j = 1:size(S, 2)
            L(i, j) = feval(sprintf('sagittal_library_dq%i', n), ([S(i, j); V(i, j)]));
        end
    end

    surface(ax, S, V, L);
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [indexes]=train_test_random_new(y,n,nall)
% function to ramdonly select training samples and testing samples from the
% whole set of ground truth.
% alltrain is the ground truth
% % clc
% alltrain = alltrain';
% y = alltrain(2,:); % Indiana
% y = alltrain(1,:); % Salinas
% yindex = [];
% train = [];
% lys = 0;
K = max(y);
% pK = 0.0001,

% generate the  training set
indexes = [];
for i = 1:K
    index1 = find(y == i);
    per_index1 = randperm(length(index1));
    if length(index1)>n
        indexes = [indexes ;index1(per_index1(1:n))'];
    else
        indexes = [indexes ;index1(per_index1(1:round(length(index1)/2)))'];
    end
end
indexes = indexes(:);
indexes_all = [1:length(y)];
indexes_all(indexes) = [];
n_new = nall - length(indexes);
per_indexall = randperm(length(indexes_all));
indexes_new = indexes_all(per_indexall(1:n_new));
indexes = [indexes;indexes_new'];
indexes = indexes(:);

% indexes = indexes';
% train = y(indexes);
% y(indexes) = [];
% test = y;

% for k_iter = 1:K
%     index_k = y == k_iter;
%     index_k = find(index_k);
%     if length(index_k) > n
%         index_k_random =  ceil(length(index_k).*rand(n,1));
%         index_k_random = sort(index_k_random);
%         index_k_random1 = [index_k_random(2:n);index_k_random(1)];
%         resid = index_k_random1 - index_k_random;
%         resid0 = resid == 0;
%         index_k_random(resid0) = [];
%         train_k = alltrain(:,index_k(index_k_random));
%         yindex = [yindex,index_k(index_k_random)];
%     else
%         n1 = ceil(length(index_k)/2);
%         index_k_random =  ceil(length(index_k).*rand(n1,1));
%         index_k_random = sort(index_k_random);
%         index_k_random1 = [index_k_random(2:n1);index_k_random(1)];
%         resid = index_k_random1 - index_k_random;
%         resid0 = resid == 0;
%         index_k_random(resid0) = [];
%         train_k = alltrain(:,index_k(index_k_random));
%         yindex = [yindex,index_k(index_k_random)];
%     end
%     train = [train,train_k];
% end
%
% trainold = alltrain;
% alltrain(:,yindex) = [];
% test = alltrain;
% for i = 1:K
%     ly = length(find(y==i));
%     nf = ceil(pK*ly);
%     f = lys + ceil(ly.*rand(nf,1));
%     f = sort(f);
%     findex = [];
%     for j = 1:(length(f)-1)
%         if f(j) == f(j+1)
%             findex = [findex,j];
%         end
%     end
%     f(findex) = [];
%     train10 = alltrain(:,f);
%     yindex = [yindex,f'];
%     train = [train,train10];
%     lys = lys + ly;
% end
%
% alltrain(:,yindex) = [];
% test = alltrain;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%% Simplicial Complex in Two Dimensions
%
% We dsecribe the data structure of the simplicial complex associated to a
% two dimensional trianglulation give by |node,elem| . The |node| records
% the coordinates of vertices and |elem| is the pointer from local to
% global incices of vertices. See <meshdoc.html meshdoc> for details.
%
% A brief summary of ordering and orientation
%
% * edge: asecond ordering, i.e. edge(:,1)<edge(:,2)
%
% * elem: positive ordering or ascend ordering. The default one is positive
% ordering and the asecond ordering is used for edge elements.
%
% * local edge:
%   consistent orientation [2 3; 3 1; 1 2];
%   asecond orientation    [2 3; 1 3; 1 2];
% The default one is the consistent orientation which will be
% counter-clockwise if elem is positive ordered. The asecond orientation
% will be used for edge elements.
%
% * elem2edgeSign: records the inconsistency of the edge
%  orientation and the induced orientation.
%
% Local edge is the induced orientation [2 3; 3 1; 1 2].
%  elem2edgeSign(1:NT,1:4) records the elementwise inconsistency of local
%  edge orientation and global edge orientation. It can be obtained from
%  |dofedge|. This is used to correct the consistency of the local basis
%  with the global basis.
%
% Both elem and local edges are ascend ordering. Then elem2edgeSign =
%  [1 -1 1] records the inconsistency of the edge orientation and the
%  induced orientation. This is used to construct differential operators.
%
% Function to call
%
%  [elem2edge,edge,elem2edgeSign,edgeSign] = dofedge(elem);
%
% Functions to read on the usage
%
%  PoissonRT0;

%%
% The basic data structure of a mesh consists of node and elem:
node = [0,0; 1,0; 1,1; 0,1];    % nodes
elem = [2,3,1; 4,1,3];          % elements
%%
% In iFEM, |N, NE, NT| represents the muber of vertice, edges, triangles
% respectively.
%
N = size(node,1); NT = size(elem,1); % NE = size(edge,1);

%%
% The corresponding simplicial complex consists of vertices, edges and
% triangles. We shall discuss the following three issues:
%
% * *Indexing* of simplexes
% * *Ordering* of vertices of simplexes
% * *Orientatoin* of simplexes
%
% The indexing and ordering are related and the ordering and orientation
% are mixed together. However the indexing has nothing to do with the
% orientation. The indexing and ordering are the combinarotry structure,
% i.e. only |elem| is needed, while the orientation also depends on |node|,
% the geometry emembdding of vertices.
%
% For indexing, ordering and orientation, there are always local and global
% version. The relation between the local and global version is the most
% complicated issue.

%% Indexing of Simplexes
%
% The indexing refers to the numbering of simplexes, e.g., which edge is
% numbered as the first one. There are two types of the indexing: local and
% global. Each simplex in the simplicial complex has a unique index which
% is called the global index. In one triangle, the three vertices and three
% edges have their local index from 1:3.
%
% In the assembling procedure of finite element methods, an element-wise
% matrix using local indexing is first computed and then assembled to get a
% big matrix using global indexing. Thus the pointer from the local
% indexing to the global indexing is indispensible. For bases independent of
% the ordering and orientation, e.g., |P1| and |P2| elements, the pointer
% is sufficient, otherwise, the inconsistency of local ordering/orientation
% and global ordering/orientation should be taken into account.
%
% *Indexing pointers of vertices*
%
% The |NT| by 3 matrix |elem| is indeed the pointer from the local to the
% global indices of vertices of triangles. For example |elem(t,1)=25| means
% the first vertex of the triangle t is the 25-th vertex.
%
% Similiary, the |NE| by 2 matrix |edge| records the pointer from the local
% to the global indices of vertices of edges.
%
% *Local indexing of edges*
%
% The triangle constists of three vertices indexed as [1,2,3]. Each
% triangle contains three edges. There are two indexing schemes for edges.
%
% * Opposite indexing     |locEdge = [2,3; 3,1; 1,2]|
%
% In |locEdge|, the |i-th| edge is opposite to the |i-th| verices and thus
% called _opposite indexing_.
%
% * Lexicogrphic indexing |locEdgel = [1,2; 1,3; 2,3]|
%
% In |locEdgel|, the indexing is induced from the lexicographic ordering of
% the three edges.
%
% For 2-D triangulations, *we shall always chose opposite indexing*. The
% lexicographic indexing is mainly used in the construction of |face2edge|
% of 3-D triangulations; see <sc3doc.html Simplicial Complex in Three
% Dimensions> for details. Note that the ordering of vertices of each edge
% will not change the indexing. For example, |locEdge = [2,3; 1,3; 1,2]|
% use the same opposite indexing but different ordering. Chosing |[1 3]| or
% |[3 1]| for the second edge will depend on the consideration of
% orientation and ordering.
%
% *Global indexing of edges*
%
% One can easily collect all edges elementwise. The issue is the
% duplication. For example, each interior edge will be counted twice. The
% |unique| funciton is applied such that each edge has a unique global
% index.
%
totalEdge = uint32([elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])]);
sortedTotalEdge = sort(totalEdge,2);
[edge,tempvar,je] = unique(sortedTotalEdge,'rows');
NE = size(edge,1);
%%
% *Edge pointer*
%
% |elem2edge(1:NT,1:3)| records the pointer from the local index to the
% global index of edges. For example, |elem2edge(t,1)| = 10 means the first
% edge of triangle |t| (which is formed by [2 3] vertices of |t|) is the
% |10-th| one in the |edge| array.
%
% Such information is stored in the third output of |unique| function.
%%
elem2edge = uint32(reshape(je,NT,3));
%%
% Note that the pointer |elem2edge| depends on the local indexing of edges
% used in the generation of |totalEdge|. Here the opposite indexing of
% three local edges is used.

%% Ordering of Vertices
%
% We discuss the ordering of vertices of simplexes. Again there are local
% ordering and global ordering. They may not be consistent and a sign array
% is used to record the inconsistency if any.
%
% The local ordering refers to the ordering of local veritces of a simplex.
% The local ordering could be used in the formulation of the local basis
% and thus the ordering does matter.
%
% The global ordering refers to the ordering of the global index of
% vertices of a simplex.
%
% *elem*. The local ordering is always [1,2,3]. Any permutation of three
% veritces of a triangle still represents the same triangle. Such
% freedom provide a room to record more information like:
%
% * global ordering of vertices
% * orientation of triangles
% * refinement rule
%
% Two types of ordering of |elem| is of particular importance
%
% * Positive ordering
% * Ascend ordering
%
% In the positive ordering, the three vertices are ordered such that the
% signed area, det(v12,v13), is positive. If |elem| is not positive
% ordered, |elem = fixorder(node,elem)| will compute the signed area by
% |simplexvolume(node,elem)| and switch the vertices for triangles with
% negative areas.
%
% For 2-D triangulations, three vertices of a triangle in 2-D is sorted
% counter-cloclwise and the first vertex is chosen as the newest vertex.
% Such ordering enables the efficient implementation of local refinement
% and coarsening in 2-D; see <bisectdoc.html Bisection in Two Dimensions>
% and <coarsendoc.html Coarsening in Two Dimensions>. *Such ordering scheme
% is the default choice and used in most places*.
%
%
% In ascend ordering, the vertices of |elem| is sorted such that
% |elem(t,1)<elem(t,2)<elem(t,3)|. This can be easily achieved by |elem =
% sort(elem,2)|. Howevery, one has to rotate the boundary flag
% accordingly using |sortelem|.
bdFlag = setboundary(node,elem,'Dirichlet');
display('Before rotation'); display(elem); display(bdFlag);

[elem,bdFlag] = sortelem(elem,bdFlag);
display('After rotation'); display(elem); display(bdFlag);

%%
% Ascend ordering will benefit the construction of local bases for high
% order basis or basis with orientation.
%
% We may switch the default positive ordering of |elem| to ascend ordering
% when generating data structure for finite element basis. However such
% sorting is always hidden in the subroutines when a finite element basis
% requiring ordering is generated; see |PoissonRT0| and |PoissonBDM1|.

%%
% *edge*. The global ordering of edges is always ascended, i.e.
edge(:,1) < edge(:,2);
%%
% Indeed in the generation of |edge|, the |totalEdge| is sorted to the
% ascend ordering such that |unique| can be applied.
%%
%
% Recall that we always use the opposite indexing of edges. We could use
% _ascend ordering_ , i.e.
locEdge = [2 3; 1 3; 1 2]; % Ascend ordering of local edges
%%
% or the _consistent (orientation) ordering_
locEdge = [2 3; 3 1; 1 2]; % Consistent ordering of local edges
%%
% It is consistent since the local ordering orientation of edges is
% consistent with the induced orientation of the three edges; see the
% discussion on the orientation.
%
% There might be an inconsistency between the local and global ordering
% (even for the consistent orientation ordering). That is
% |edge(elem2edge(t,1),1)| may not be smaller than
% |edge(elem2edge(t,1),2)|. It will be more clear from the discussion of
% the corresponding orientation.

%% Orientation of Simplexes
%
% The orientation of a triangle is either positive or negative, the
% orientation of an edge is given by a tangential or normal vector. A
% normal vector is obtained by rotate a given tangential vector by 90
% degree clockwise. For example, when edges are given counter clockwise
% orientation, the corresponding normal vector is the outwards normal
% vector.
%
% The orientation of a d-simplex will induce an orientation of its d-1
% subcomplex and is called _induced orientation_. For example, a positive
% orientated triangle will induce the counter clockwise orientation of its
% three edges.
%
% The ordering of vertices will naturally introduce an orientation and will
% be called _ordering orientation_. More specifically
%
% * The vector from |edge(:,1)| to |edge(:,2)| defines an orientation of
% edges.
% * The sign of det(v12,v13) defines an orientation of the triangle.
%
% The orientation of a simplex in the simplicial complex should be
% uniquely determined which will be called _global orientation_. It can be
% chosen as the global ordering orientation but not always the case.
%
% Inside one triangle, the local orientation of three edges is more
% involved. The local ordering of edges will introduce a local ordering
% orientation. The orientation of the triangle will also induce an induced
% orientation. The local ordering orientation is used in the computing of
% local bases and the induced orientation is used when computing the
% differential operator locally. They may or may not be consisitent with
% the global orientation of edges.
%
% In general, there will be an inconsistency of the following types of
% orientation and apporipate data structure should be constructured to
% record such inconsistency.
%
% * a global orientation
% * the global ordering orientation
% * the local ordering orientation
% * the local induced orientation
%
% *elem*. The orientation of a triangle is either positive or negative. For
% the global ordering orientation, it is the sign of the signed area
% (output of |simplexvolume|).
%
[Dlambda,area,elemSign] = gradbasis(node,elem);
%%
% In the output of |gradbasis|, |area| is always positive and an additional
% array |elemSign| is used to record the sign of the signed area.
%
% |Dlambda(t,:,k)| is the gradient of $\lambda_k$. Therefore the outward
% normal direction of the kth face can be obtained by |-Dlambda(t,:,k)|
% which is independent of the ordering and orientation.
%
% *edge*. For 2-D triangulations, *we shall always chose global ordering
% orientation*, i.e., from the lower index to the bigger index.
%
% The local ordering orientation is implicitly used when computing finite
% element basis in each element. For example, the edge element on edge |[i j]|
% in |locEdge| is defined as
%
% $$\phi_{i,j} = 2(\lambda_i \nabla \lambda_j - \lambda_j \nabla \lambda_i).$$
%
% Permutation of |[i j]| to |[j i]| will change the sign of the basis. Note
% that this is locally, i.e., element by element.
%
% The global basis associated to an edge, however, depends only on the
% global orientation of this edge. We introduce |elem2edgeSign(1:NT, 1:3)|
% to record the inconsistency of a local ordering orientation and a global
% orientation.
%
% For the consistent local ordering [2 3; 3 1; 1 2] and global ascend
% ordering orientation, the elem2edgeSign can be generated as follows:

elem2edgeSign = ones(NT,3,'int8');
totalEdge = uint32([elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])]);
idx = (totalEdge(:,1)>totalEdge(:,2));
elem2edgeSign(idx) = -1;

%%
% There is one more inconsistency between the induced orientation and the
% global orientation of edges. If a triangle is positive orientated, the
% induced edge orientation should be given by the outwards (relative to a
% triangle) normal vector. This induced orientation may not be consistent
% with the global orientation of edges.
%
% It depends on the ordering of |elem| and |locEdge|. When |elem| is
% positive ordered and |locEdge| is consistently ordered, it is.
%
% For ascend ordering of |elem| and |edge|, we denote the direction as +1
% if the direction of an edge is the same with the induced direction in a
% certain elem, and -1 otherwise. Then the consistency is given by
elem2edgeSign = [+1 -1 +1];

%%
% The second sign is -1 because the local edge in ascend ordering is [1
% 3] not [3 1].
%%
% The |elem2edgeSign| will be used when assembling differential operators.
% For example, when computing |div| operators on a positive orientated
% triangle, the edge should have outwards normal direction or equivalently
% the counter clockwise orientation.
%%
% We summarize the two popular ordering and orientation schemes below.

%% Positive Ordering and Orientation
%
% The vertice of the |eleme| is sorted such that the area is always
% positive. i.e. the three vertices of the elem are ordered
% counter-clockwisely. Furthermore the first vertex is always the newest
% vertex of the triangle for the easy of local mesh refinement and
% coarsening.
%
% The local edge is using opposite indexing and consistent ordering
%
% |locEdge = [2,3; 3,1; 1,2]|
%
% The ascend ordering is used for global edges, i.e., |edge| is sorted s.t.
%
% |edge(:,1) < edge(:,2)|
%
% The inconsistency of the orientation is recorded in|elem2edgeSign|.
%
% *Where to use*: This is the default ordering and orientation scheme.
%
%% Ascend Ordering and Orientation
%
% *Ascend ordering*. The vertices of |elem| is sorted such that
%
% |elem(t,1)<elem(t,2)<elem(t,3)|
%
% The local edge is also in the ascend ordering
%
% |locEdge = [2,3; 1,3; 1,2]|
%
% The ascend ordering is used for edges, i.e., |edge| is sorted s.t.
%
% |edge(:,1) < edge(:,2)|
%
% It is easy to see the benefit of ascend ordering: the ordering of local
% edges is consistent with the global ones and so is the corresponding
% orientation.
%
% *Orientation*. We choose the ordering orientation.
%
% |elem: sign(det(v12,v13))|
%
% |edge: from the node with smaller global index to bigger one|
%
% For |edge|, the orientation of the global ordering and the local
% ordering is consistent. But they are not consistent with the induced
% orientation inside one triangle. Such inconsistency is recorded in
%
% |elem2edgeSign = [+1 -1 +1]|
%
%
% *Where to use*: for H(curl) and H(div) elements and high order (cubic and
% above) H(grad) elements.

%% An Example
%
node = [0,0; 1,0; 1,1; 0,1];    % nodes
elem = [2,3,1; 4,1,3];          % elements
% [node,elem] = uniformrefine(node,elem);
%%
% Poistive ordering and orientation
[elem2edge,edge,elem2edgeSign] = dofedge(elem);
% plot
set(gcf,'Units','normal');
set(gcf,'Position',[0,0.25,0.25,0.25]);
showmesh(node,elem);
findnode(node);
findelem(node,elem);
findedge(node,edge);
%%
display(elem);
display(edge);
display(elem2edge);
display(elem2edgeSign);

%%
% Ascend Ordering and Orientation
bdFlag = setboundary(node,elem,'Dirichlet');
[elem,bdFlag] = sortelem(elem,bdFlag);
[elem2edge,edge,elem2edgeSign] = dofedge(elem);
display(elem);
display(edge);
display(elem2edge);
display(elem2edgeSign);
elem2edgeSign = [1,-1,1];

repo: Aero-matlab


text:
function ExportAedat2Imu6(aedat)

%{
This function exports data to a .aedat file.
The .aedat file format is documented here:

http://inilabs.com/support/software/fileformat/
%}

dbstop if error

if ~exist('aedat', 'var')
	error('Missing input')
end

% Create the file
if ~isfield(aedat.exportParams, 'filePath')
    error('Missing file path and name')
end

f = fopen(aedat.exportParams.filePath, 'w', 'b');

% Simple - events only - assume DAVIS

% CRLF \r\n is needed to not break header parsing in jAER
fprintf(f,'#!AER-DAT2.0\r\n');
fprintf(f,'# This is a raw AE data file created by an export function in the AedatTools library\r\n');
fprintf(f,'# Data format is int32 address, int32 timestamp (8 bytes total), repeated for each event\r\n');
fprintf(f,'# Timestamps tick is 1 us\r\n');
% Put the source in - use an override if it has been given
if isfield(aedat.exportParams, 'source')
    fprintf(f,['# AEChip: ' aedat.exportParams.source '\r\n']);
else
    fprintf(f,['# AEChip: ' aedat.info.source '\r\n']);
end
fprintf(f,'# End of ASCII Header\r\n');

% DAVIS

accelX = aedat.data.imu6.accelX; % conversion from g to full scale, and shift bits
accelX = int16(accelX * 8192); % conversion from g to full scale 16 range
accelX = [accelX zeros(aedat.data.imu6.numEvents, 1, 'int16')];
accelX = accelX';
accelX = accelX(:);
accelX = typecast(accelX, 'uint32');
accelX = bitshift(accelX, 12); % shift bits
accelX = accelX + 2 ^ 31 + 2 ^ 11; % imu flag bits

accelY = aedat.data.imu6.accelY; % conversion from g to full scale, and shift bits
accelY = int16(accelY * 8192); % conversion from g to full scale 16 range
accelY = [accelY zeros(aedat.data.imu6.numEvents, 1, 'int16')];
accelY = accelY';
accelY = accelY(:);
accelY = typecast(accelY, 'uint32');
accelY = bitshift(accelY, 12); % shift bits
accelY = accelY + 2 ^ 31 + 2 ^ 11; % imu flag bits

accelZ = aedat.data.imu6.accelZ; % conversion from g to full scale, and shift bits
accelZ = int16(accelZ * 8192); % conversion from g to full scale 16 range
accelZ = [accelZ zeros(aedat.data.imu6.numEvents, 1, 'int16')];
accelZ = accelZ';
accelZ = accelZ(:);
accelZ = typecast(accelZ, 'uint32');
accelZ = bitshift(accelZ, 12); % shift bits
accelZ = accelZ + 2 ^ 31 + 2 ^ 11; % imu flag bits

temp = aedat.data.imu6.temperature; % conversion from g to full scale, and shift bits
temp = int16((temp - 35) * 340); % conversion from K to full scale 16 range
temp = [temp zeros(aedat.data.imu6.numEvents, 1, 'int16')];
temp = temp';
temp = temp(:);
temp = typecast(temp, 'uint32');
temp = bitshift(temp, 12); % shift bits
temp = temp + 2 ^ 31 + 2 ^ 11; % imu flag bits

gyroX = aedat.data.imu6.gyroX; % conversion from g to full scale, and shift bits
gyroX = int16(gyroX * 65.5); % conversion from g to full scale 16 range
gyroX = [gyroX zeros(aedat.data.imu6.numEvents, 1, 'int16')];
gyroX = gyroX';
gyroX = gyroX(:);
gyroX = typecast(gyroX, 'uint32');
gyroX = bitshift(gyroX, 12); % shift bits
gyroX = gyroX + 2 ^ 31 + 2 ^ 11; % imu flag bits

gyroY = aedat.data.imu6.gyroY; % conversion from g to full scale, and shift bits
gyroY = int16(gyroY * 65.5); % conversion from g to full scale 16 range
gyroY = [gyroY zeros(aedat.data.imu6.numEvents, 1, 'int16')];
gyroY = gyroY';
gyroY = gyroY(:);
gyroY = typecast(gyroY, 'uint32');
gyroY = bitshift(gyroY, 12); % shift bits
gyroY = gyroY + 2 ^ 31 + 2 ^ 11; % imu flag bits

gyroZ = aedat.data.imu6.gyroZ; % conversion from g to full scale, and shift bits
gyroZ = int16(gyroZ * 65.5); % conversion from g to full scale 16 range
gyroZ = [gyroZ zeros(aedat.data.imu6.numEvents, 1, 'int16')];
gyroZ = gyroZ';
gyroZ = gyroZ(:);
gyroZ = typecast(gyroZ, 'uint32');
gyroZ = bitshift(gyroZ, 12); % shift bits
gyroZ = gyroZ + 2 ^ 31 + 2 ^ 11; % imu flag bits

allData = [accelX accelY accelZ temp gyroX gyroY gyroZ];
allData = allData';
allData = allData(:);

timeStamps = uint32(aedat.data.imu6.timeStamp(:));
timeStamps = repmat(timeStamps', 7 , 1);
timeStamps = timeStamps(:);

output = zeros(1, 2 * aedat.data.imu6.numEvents * 7, 'uint32'); % allocate horizontal vector to hold output data
output(1 : 2 : end) = allData;
output(2 : 2 : end) = timeStamps; % set even elements to timestamps

% write addresses and timestamps
count=fwrite(f, output, 'uint32') / 2; % write 4 byte data
fclose(f);
fprintf('wrote %d events to %s\n',count, aedat.exportParams.filePath);

author: panpanfei


text:
function [data] = fiff_read_evoked(fname,setno)
%
% [data] = fiff_read_evoked(fname,setno)
%
% Read one evoked data set
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.11  2009/03/30 11:37:37  msh
%   Added copying of measurement info blocks from the original like in mne_browse_raw
%
%   Revision 1.10  2008/03/28 13:42:01  msh
%   Accommodate evoked files from xplotter, which may have channel information
%   local to the evoked block
%
%   Revision 1.9  2007/11/13 10:55:32  msh
%   Specify the second argument to all calls to the exist function
%
%   Revision 1.8  2006/10/14 11:47:32  msh
%   Omitted listing of nepoch from the diagnostic output
%
%   Revision 1.7  2006/09/21 19:27:16  msh
%   fiff_read_evoked did not tolerate missing comment
%
%   Revision 1.6  2006/05/03 19:03:19  msh
%   Eliminated the use of cast function for Matlab 6.5 compatibility
%
%   Revision 1.5  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.4  2006/04/18 20:44:46  msh
%   Added reading of forward solution.
%   Use length instead of size when appropriate
%
%   Revision 1.3  2006/04/12 13:13:51  msh
%   Added fiff_find_evoked.m
%   Use aspect_kind field name instead of aspect_type
%
%   Revision 1.2  2006/04/12 10:29:02  msh
%   Made evoked data writing compatible with the structures returned in reading.
%
%   Revision 1.1  2006/04/10 23:26:54  msh
%   Added fiff reading routines
%
%

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

me='MNE:fiff_read_evoked';

if nargin == 1
    setno = 1;
elseif nargin ~= 2
    error(me,'Incorrect number of arguments');
end
if setno <= 0
    error(me,'Data set selector must be positive');
end
%
%   Open the file
%
fprintf(1,'Reading %s ...\n',fname);
[ fid, tree ] = fiff_open(fname);
%
%   Read the measurement info
%
[ info, meas ] = fiff_read_meas_info(fid,tree);
info.filename = fname;
%
%   Locate the data of interest
%
processed = fiff_dir_tree_find(meas,FIFF.FIFFB_PROCESSED_DATA);
if length(processed) == 0
    fclose(fid);
    error(me,'Could not find processed data');
end
%
evoked = fiff_dir_tree_find(meas,FIFF.FIFFB_EVOKED);
if length(evoked) == 0
    fclose(fid);
    error(me,'Could not find evoked data');
end
%
%   Identify the aspects
%
naspect = 0;
is_smsh = [];
for k = 1:length(evoked)
    sets(k).aspects = fiff_dir_tree_find(evoked(k),FIFF.FIFFB_ASPECT);
    sets(k).naspect = length(sets(k).aspects);
    if sets(k).naspect > 0
        is_smsh = [ is_smsh zeros(1,sets(k).naspect) ];
        naspect = naspect + sets(k).naspect;
    end
    saspects  = fiff_dir_tree_find(evoked(k), FIFF.FIFFB_SMSH_ASPECT);
    nsaspects = length(saspects);
    if nsaspects > 0
        sets(k).naspect = sets(k).naspect + nsaspects;
        sets(k).aspects = [ sets(k).aspects saspects ];
        is_smsh = [ is_smsh ones(1,sets(k).naspect) ];
        naspect = naspect + nsaspects;
    end
end
fprintf(1,'\t%d evoked data sets containing a total of %d data aspects in %s\n',length(evoked),naspect,fname);
if setno > naspect || setno < 1
    fclose(fid);
    error(me,'Data set selector out of range');
end
%
%   Next locate the evoked data set
%
p = 0;
goon = true;
for k = 1:length(evoked)
    for a = 1:sets(k).naspect
        p = p + 1;
        if p == setno
            my_evoked = evoked(k);
            my_aspect = sets(k).aspects(a);
            goon = false;
            break;
        end
    end
    if ~goon
        break;
    end
end
%
%   The desired data should have been found but better to check
%
if ~exist('my_evoked','var') || ~exist('my_aspect','var')
    fclose(fid);
    error(me,'Desired data set not found');
end
%
%   Now find the data in the evoked block
%
nchan = 0;
sfreq = -1;
q = 0;
for k = 1:my_evoked.nent
    kind = my_evoked.dir(k).kind;
    pos  = my_evoked.dir(k).pos;
    switch kind
        case FIFF.FIFF_COMMENT
            tag = fiff_read_tag(fid,pos);
            comment = tag.data;
        case FIFF.FIFF_FIRST_SAMPLE
            tag = fiff_read_tag(fid,pos);
            first = tag.data;
        case FIFF.FIFF_LAST_SAMPLE
            tag = fiff_read_tag(fid,pos);
            last = tag.data;
        case FIFF.FIFF_NCHAN
            tag = fiff_read_tag(fid,pos);
            nchan = tag.data;
        case FIFF.FIFF_SFREQ
            tag = fiff_read_tag(fid,pos);
            sfreq = tag.data;
        case FIFF.FIFF_CH_INFO
            q = q+1;
            tag = fiff_read_tag(fid,pos);
            chs(q) = tag.data;
    end
end
if ~exist('comment','var')
    comment = 'No comment';
end
%
%   Local channel information?
%
if nchan > 0
    if ~exist('chs','var')
        fclose(fid);
        error(me, ...
            'Local channel information was not found when it was expected.');
    end
    if length(chs) ~= nchan
        fclose(fid);
        error(me, ...
            'Number of channels and number of channel definitions are different');
    end
    info.chs   = chs;
    info.nchan = nchan;
    fprintf(1, ...
        '\tFound channel information in evoked data. nchan = %d\n',nchan);
    if sfreq > 0
        info.sfreq = sfreq;
    end
end
nsamp = last-first+1;
fprintf(1,'\tFound the data of interest:\n');
fprintf(1,'\t\tt = %10.2f ... %10.2f ms (%s)\n',...
    1000*double(first)/info.sfreq,1000*double(last)/info.sfreq,comment);
if ~isempty(info.comps)
    fprintf(1,'\t\t%d CTF compensation matrices available\n',length(info.comps));
end
%
% Read the data in the aspect block
%
nepoch = 0;
for k = 1:my_aspect.nent
    kind = my_aspect.dir(k).kind;
    pos  = my_aspect.dir(k).pos;
    switch kind
        case FIFF.FIFF_COMMENT
            tag = fiff_read_tag(fid,pos);
            comment = tag.data;
        case FIFF.FIFF_ASPECT_KIND
            tag = fiff_read_tag(fid,pos);
            aspect_kind = tag.data;
        case FIFF.FIFF_NAVE
            tag = fiff_read_tag(fid,pos);
            nave = tag.data;
        case FIFF.FIFF_EPOCH
            nepoch = nepoch + 1;
            tag = fiff_read_tag(fid,pos);
            epoch(nepoch) = tag;
    end
end
if ~exist('nave','var')
    nave = 1;
end
fprintf(1,'\t\tnave = %d aspect type = %d\n',...
    nave,aspect_kind);
if nepoch ~= 1 && nepoch ~= info.nchan
    fclose(fid);
    error(me,'Number of epoch tags is unreasonable (nepoch = %d nchan = %d)',nepoch,info.nchan);
end
%
if nepoch == 1
    %
    %   Only one epoch
    %
    all = epoch(1).data;
    %
    %   May need a transpose if the number of channels is one
    %
    if size(all,2) == 1 && info.nchan == 1
        all = all';
    end
else
    %
    %   Put the old style epochs together
    %
    all = epoch(1).data';
    for k = 2:nepoch
        all = [ all ; epoch(k).data' ];
    end
end
if size(all,2) ~= nsamp
    fclose(fid);
    error(me,'Incorrect number of samples (%d instead of %d)',...
        size(all,2),nsamp);
end
%
%   Calibrate
%
for k = 1:info.nchan
    cals(k) = info.chs(k).cal;
end
all = diag(cals)*all;
%
%   Put it all together
%
data.info = info;
data.evoked.aspect_kind = aspect_kind;
data.evoked.is_smsh     = is_smsh(setno);
if exist('nave','var')
    data.evoked.nave = nave;
else
    data.evoked.nave  = 1;
end
data.evoked.first = first;
data.evoked.last  = last;
if exist('comment','var')
    data.evoked.comment = comment;
end
%
%   Times for convenience and the actual epoch data
%
data.evoked.times = double(data.evoked.first:1:data.evoked.last)/data.info.sfreq;
data.evoked.epochs = all;

fclose(fid);

return;

    function [tag] = find_tag(node,findkind)

        for p = 1:node.nent
            kind = node.dir(p).kind;
            pos  = node.dir(p).pos;
            if kind == findkind
                tag = fiff_read_tag(fid,pos);
                return;
            end
        end
        tag = [];
        return
    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [A, b] = sllinrega(X, Y, varargin)
%SLLINREGA Performs Augmented Multivariate Linear Regression
%
% $ Syntax $
%   - [A, b] = sllinrega(X, Y, ...)
%
% $ Arguments $
%   - X:        The sample matrix of x
%   - Y:        The sample matrix of y
%   - A:        The solved transform matrix
%   - b:        The solved shift vector
%
% $ Description $
%   - [A, b] = sllinrega(X, Y, ...) solves the regression problem given
%     by the following formula:
%           y = A * x + b
%     in least square error sense. The samples are stored in X and Y
%     in column-wise manner.
%     You can specify properties for regression as in sllinreg.
%
% $ Remarks $
%   - The implementation is based on sllinreg with an augmented
%     formulation as follows:
%           y = [A, b] * [x; 1]
%
% $ History $
%   - Created by Dahua Lin, on Sep 15th, 2006
%

%% parse and verify input arguments

if nargin < 2
    raise_lackinput('sllinrega', 2);
end

if ~isnumeric(X) || ~isnumeric(Y) || ndims(X) ~= 2 || ndims(Y) ~= 2
    error('sltoolbox:invalidarg', ...
        'The X and Y should be both 2D numeric matrices');
end

%% main

% augment formulation
[dx, nx] = size(X);
Xa = [X; ones(1, nx)];

% solve
Aa = sllinreg(Xa, Y, varargin{:});
clear Xa;

% extract
A = Aa(:, 1:dx);
b = Aa(:, dx+1);

author: lmthang


text:
function r=rot4y(theta)
    % r = rot4y(theta)
    %
    % roty produces a 4x4 rotation matrix representing
    % a rotation by theta radians about the y axis.
    %
    %	Argument definitions:
    %
    %	theta = rotation angle in radians
    c = cos(theta);
    s = sin(theta);
    r = [c  0  s 0;
        0  1  0 0;
        -s  0  c 0;
        0  0  0 1];

author: CelsoReyes


text:
classdef AltitudeOptimizationVariable < AbstractOptimizationVariable
    %AltitudeOptimizationVariable Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) AltitudeTermCondition = AltitudeTermCondition(0);

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        useTf(1,1) = false;
    end

    methods
        function obj = AltitudeOptimizationVariable(varObj)
            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = [];

            if(obj.useTf)
                x = obj.varObj.altitude;
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            lb = obj.lb(obj.useTf);
            ub = obj.ub(obj.useTf);
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lwrBnd;
            ub = obj.uprBnd;
        end

        function setBndsForVariable(obj, lb, ub)
            obj.lb = lb;
            obj.ub = ub;
        end

        function useTf = getUseTfForVariable(obj)
            useTf = obj.useTf;
        end

        function setUseTfForVariable(obj, useTf)
            obj.useTf = useTf;
        end

        function updateObjWithVarValue(obj, x)
            obj.varObj.altitude = x;
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
            nameStrs = {sprintf('Event %i Altitude Termination Condition', evtNum)};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function GUI = sessionGUI_plotTimeCourse(rois, scans, viewer);
%
%  GUI = sessionGUI_plotTimeCourse(<rois>, <scans=scan group>, <viewer>);
%
% Open a Time Course UI for the selected ROIs.
%
% If a viewer is open, will open a time course UI for the currently-selected ROI
% in the selected viewer. If no viewers are open, will open it for the
% selected ROI in the ROI listbox in the session GUI. These are both
% over-ridden if an roi list and/or viewer are manually specified.
%
%
%
% ras, 07/06.
mrGlobals2;

if notDefined('viewer')
    if isempty(GUI.settings.viewer) | GUI.settings.viewer==0
        viewer = 0;
    else
         viewer = GUI.viewers(GUI.settings.viewer);
    end
end

if notDefined('rois')
    if viewer==0 % no viewer, get from saved list
        rois = guiGet('selectedRois');
    else
        rois = {mrViewGet(viewer, 'selectedRoi')};
    end
end

if ~checkfields(GUI, 'controls', 'tcui')
    GUI.controls.tcui = [];
end

% check if an event-related scan group is assigned.
% If so, use that; otherwise use the GUI selected scans.
if notDefined('scans')
    scans = GUI.settings.scan;
    dt = GUI.settings.dataType;
    if isfield(dataTYPES(dt).scanParams(scans(1)), 'scanGroup') & ...
            ~isempty(dataTYPES(dt).scanParams(scans(1)).scanGroup)
        [scans dt] = er_getScanGroup(INPLANE{1});
    end
end

if notDefined('dt'), dt = guiGet('DataType'); end

% special case: for GLMs data type, the time series will be
% in a different d.t. We don't need to ask the user to confirm
% in this case, it should be understood.
if isequal(dataTYPES(dt).name, 'GLMs')
	queryFlag = 1;
else
	queryFlag = 0;
end

for i = 1:length(rois)
    tc = timeCourseUI(INPLANE{1}, rois{i}, scans, dt, queryFlag);
    GUI.controls.tcui(end+1) = tc.ui.fig;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [x, y, model, pca_model] = train(total)
% X = features
% Y = label
x = [];
y = {};

for n = 1:length(total)-1
    display(['Reading img',num2str(n)])
    for i = 1:8
        for j = 1:8
            y_temp = total{n+1,3}{i,j};
            if strcmp(y_temp,'')
                continue;
            end
            fields = fieldnames(total{n+1,2}{i,j});
            x_temp = zeros(1,length(fields));
            for f = 1:length(fields)
                x_temp(f) = [getfield(total{n+1,2}{i,j},fields{f})];
            end

            x = [ x ; x_temp ];
            y = [ y ; y_temp ];
        end
    end
end

[x, pca_x,pca_coeff] = preprocess(x);
save('pca_coeff.mat','pca_coeff')

model = fitcknn(x,y,'NumNeighbors',15);
save('model.mat','model')

pca_model = fitcknn(pca_x,y,'NumNeighbors',15);
save('pca_model.mat','pca_model')


text:
function [X_scal] = test_pretreatment(X,param)

% pretreatment for test data
%
% INPUT
% X:        data matrix (samples x variables)
% param:    output data structure from data_pretreatment routine
%
% OUTPUT
% X_scal:       pretreated data matrix (samples x variables)
%
% version 1.0 - september 2009
% Davide Ballabio
% Milano Chemometrics and QSAR Research Group
% www.disat.unimib.it/chm

a = param.a;
s = param.s;
m = param.m;
M = param.M;
pret_type = param.pret_type;

if strcmp(pret_type,'cent')
    amat = repmat(a,size(X,1),1);
    X_scal = X - amat;
elseif strcmp(pret_type,'scal')
    f = find(s>0);
    smat = repmat(s,size(X,1),1);
    X_scal = zeros(size(X,1),size(X,2));
    X_scal = X(:,f)./smat(:,f);
elseif strcmp(pret_type,'auto')
    f = find(s>0);
    amat = repmat(a,size(X,1),1);
    smat = repmat(s,size(X,1),1);
    X_scal = zeros(size(X,1),size(X,2));
    X_scal(:,f) = (X(:,f) - amat(:,f))./smat(:,f);
elseif strcmp(pret_type,'rang')
    f = find(M - m > 0);
    mmat = repmat(m,size(X,1),1);
    Mmat = repmat(M,size(X,1),1);
    X_scal = zeros(size(X,1),size(X,2));
    X_scal(:,f) = (X(:,f) - mmat(:,f))./(Mmat(:,f) - mmat(:,f));
else
    X_scal = X;
end

author: kmansouri

repo: OPERA

save_path: github-repos/MATLAB/kmansouri-OPERA


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

function [A, x, b, k] = problem_1()
    m = 100;
    n = 1000;
    k = 4;
    A = spx.dict.simple.gaussian_dict(m, n);
    gen = spx.data.synthetic.SparseSignalGenerator(n, k);
    % create a sparse vector
    x =  gen.biGaussian();
    b = A*x;
end

function [dict, reps, signals, k] = problem_2()
    m = 200;
    n = 1000;
    k = 10;
    s = 500;
    dict = spx.dict.simple.gaussian_dict(m, n);
    gen = spx.data.synthetic.SparseSignalGenerator(n, k, s);
    % create a sparse vector
    reps =  gen.biGaussian();
    signals = dict*reps;
end

function test_naive_omp_1(testCase)
    [A, x, b, k] = problem_1();
    solver = spx.pursuit.single.OrthogonalMatchingPursuit(A, k);
    result = solver.solve(b);
    cmpare = spx.commons.SparseSignalsComparison(x, result.z, k);
    %cmpare.summarize();
    verifyTrue(testCase, cmpare.all_have_matching_supports(1.0));
end

function test_naive_omp_2(testCase)
    [dict, reps, signals, k] = problem_2();
    solver = spx.pursuit.single.OrthogonalMatchingPursuit(dict, k);
    ns = size(signals, 2);
    dd = size(dict, 2);
    recovered = zeros(dd, ns);
    for s=1:ns
        signal = signals(:, s);
        result = solver.solve(signal);
        recovered(:, s) = result.z;
    end
    cmpare = spx.commons.SparseSignalsComparison(reps, recovered, k);
    % cmpare.summarize();
    verifyTrue(testCase, cmpare.all_have_matching_supports(1.0));
end

function test_omp_qr_1(testCase)
    [A, x, b, k] = problem_1();
    solver = spx.pursuit.single.OrthogonalMatchingPursuit(A, k);
    result = solver.solve_qr(b);
    cmpare = spx.commons.SparseSignalsComparison(x, result.z, k);
    %cmpare.summarize();
    verifyTrue(testCase, cmpare.all_have_matching_supports(1.0));
end

function test_omp_qr_2(testCase)
    [dict, reps, signals, k] = problem_2();
    solver = spx.pursuit.single.OrthogonalMatchingPursuit(dict, k);
    ns = size(signals, 2);
    dd = size(dict, 2);
    recovered = zeros(dd, ns);
    for s=1:ns
        signal = signals(:, s);
        result = solver.solve_qr(signal);
        recovered(:, s) = result.z;
    end
    cmpare = spx.commons.SparseSignalsComparison(reps, recovered, k);
    % cmpare.summarize();
    verifyTrue(testCase, cmpare.all_have_matching_supports(1.0));
end

author: indigits

repo: sparse-plex


text:
%                           buildTestProblem.m
%
% This function creates and outputs random generated data and measurements
% according to user's choice. It is invoked in test*.m in
% order to build a test problem.
%
% Inputs:
%   m(integer): number of measurements.
%   n(integer): length of the unknown signal.
%   isComplex(boolean, default=true): whether the signal and measurement
%     matrix is complex. isNonNegativeOnly(boolean, default=false): whether
%     the signal is real and non-negative.
%   dataType(string, default='gaussian'): it currently supports
%     ['gaussian', 'fourier'].
%
% Outputs:
%   A: m x n measurement matrix/function handle.
%   xt: n x 1 vector, true signal.
%   b0: m x 1 vector, measurements.
%   At: A n x m matrix/function handle that is the transpose of A.
%
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

function [A, xt, b0, At] = buildTestProblem(m, n, isComplex, isNonNegativeOnly, dataType)
if ~exist('isComplex','var') | isempty(isComplex)
    isComplex = true;
end
if ~exist('isNonNegativeOnly','var') | isempty(isNonNegativeOnly)
    isNonNegativeOnly = false;
end
if ~exist('dataType','var')
    dataType = 'Gaussian';
end

switch lower(dataType)
    case 'gaussian'
        A = (mvnrnd(zeros(1, n), eye(n)/2, m) + isComplex * 1i * ...
            mvnrnd(zeros(1, n), eye(n)/2, m));
        At = A';
        xt = (mvnrnd(zeros(1, n), eye(n)/2) + isComplex * 1i * ...
            mvnrnd(zeros(1, n), eye(n)/2))';
        b0 = abs(A*xt);
    case 'fourier'
        %  Define the Fourier measurement operator.
        %  The operator 'A' maps an n-vector into an m-vector, then
        %  computes the fft on that m-vector to produce m measurements.

        % rips first 'length' entries from a vector
        rip = @(x,length) x(1:length);
        A = @(x) fft([x;zeros(m-n,1)]);
        At = @(x) rip(m*ifft(x),n);     % transpose of FM
        xt = (mvnrnd(zeros(1, n), eye(n)/2) + isComplex * 1i * ...
            mvnrnd(zeros(1, n), eye(n)/2))';
        b0 = abs(A(xt)); % Compute the phaseless measurements

    otherwise
        error('invalid dataType: %s',dataType);
end
end

author: tomgoldstein

repo: phasepack-matlab


text:
classdef nme_shunt_acp < mp.nme_shunt_ac & mp.form_acp

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

%     properties
%     end     %% properties

%     methods
%     end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%NS1
% Numerical solution of 2-D instationary Navier-Stokes equations.
% Rectangular domain; nonuniform Cartesian staggered grid.
% Time discretization with omega-scheme.
% Pressure-correction method.
% Central or upwind scheme for inertia terms.

% Copyright 2001 P. Wesseling
% This program and its subprograms may be freely used, modified and distributed
% under the GNU General Public License: http://www.gnu.org/copyleft/gpl.html

% Theory is given in Chapter 6 of
% 	Principles of Computational Fluid Dynamics, by P. Wesseling
% 	Springer-Verlag, Berlin etc., 2001. ISBN 3-540-67853-0
% See http://dutita0.twi.tudelft.nl/nw/users/wesseling/

% Theory is also given in Chapter 5 of
%	Computational Fluid Dynamics
%	Lecture Notes by P. Wesseling
%	Department of Applied Mathematical Analysis, ITS, TU Delft
% See http://ta.twi.tudelft.nl/nw/users/wesseling/

% This program generates Figures  and  in the Lecture Notes.

% Functions called: startup,problem_specification, grid_generation, draw_grid,
%		    viscous_matrix, grad_and_div, initial_condition,
%		    right_hand_side, inertia_matrix, errornorm, isobarplot,
%		    streamlineplot, velocity_profile, relchangeplot

startup, format compact, clear all
time0 = clock;
global geval n u0 v0 J K yu yv yseglen alpha
% ........................................Input..............................
geval = 3;	% Enter 1 for horizontal Poiseuille flow to the right
		%  or 2 for vertical Poiseuille flow
		%  or 3 for backward facing step
		%  or 4 for driven cavity
		%  or 5 for uniform flow under angle alpha without segmentation
		%  or 6 for uniform flow under angle alpha with segmentation
		%  or 7 for horizontal Poiseuille flow to the left

problem_specification	% Prepare file problem_specification.m as part of input
% User must prepare files ubd.m, vbd.m, pbd.m to specify boundary values

% .................................End of input......................

grid_generation
draw_grid
h = waitbar(0,'Computing...');
viscous_matrix		% Generate viscous matrices Bu, Bv and Masku, Maskv
grad_and_div		% Generate pressure gradient matrices Pu, Pv and
			%  divergence matrices Du, Dv
initial_condition	% Generate initial u0, v0, p0
tic
[Lp,Up] = lu([Du Dv]*[Pu;Pv]);	% LU decomposition of pressure correction matrix
tijd = toc; disp(['Pres. corr. matrix LU time = ',num2str(tijd)])
if max(abs(sum(([Du Dv]*[Pu;Pv])'))) < 10*eps
  % Do consistency check
end

nstep = floor(tend/dt);	% Number of time steps
%nstep = 2;
relchange = zeros(nstep,1);	% For plotting relative change per time step
t = 0; n = 0;
right_hand_side		% Generate right-hand sides ru, rv; has to be moved
			% inside time loop if boundary conditions time-dependent
u1 = u0; v1 = v0;
for n = 1:nstep		% Time-stepping with omega scheme
  t = n*dt;
  u0 = u1; v0 = v1;
%  right_hand_side	% Generate right-hand sides ru, rv; can be moved outside
			% time loop if boundary conditions not time-dependent

  % Navier-Stokes, predictor:
  inertia_matrix	% Generate inertia matrices and adapt right-hand sides
  rum = (ru + rum)'; rvm = (rv + rvm)'; rum = rum(:); rvm = rvm(:);
  if n == 1, tic, end
  Cu = dt*(Bu + Cu); Cv = dt*(Bv + Cv);
  u1 = (speye(length(u0)) + omega*Cu)\...
        (Masku*u0 + dt*(rum - Pu*p0) - (1-omega)*Cu*u0);
  v1 = (speye(length(v0)) + omega*Cv)\...
        (Maskv*v0 + dt*(rvm - Pv*p0) - (1-omega)*Cv*v0);

  % End of Navier-Stokes predictor

  % Stokes:
  %rum = ru'; rum = rum(:); rvm = rv'; rvm = rvm(:);
  %if n == 1, tic, end
  %u1 = Bu\(dt*(rum - Pu*p0) + Masku*u0);
  %v1 = Bv\(dt*(rvm - Pv*p0) + Maskv*v0);
  % End of Stokes

  if n == 1, tijd = toc; disp(['velocity solve time = ',num2str(tijd)]), tic,end
  dp = Up\(Lp\(Du*u1 + Dv*v1));	p0 = p0 + dp/dt;
  u1 = u1 - Pu*dp;	v1 = v1 - Pv*dp;
  if n == 1, tijd = toc; disp(['press. corr. time = ',num2str(tijd)]),end
  mm = max([norm(u1,inf),norm(v1,inf)]); relchange(n) = max([norm(u1-u0,inf)...
  /mm,norm(v1-v0,inf)/mm, norm(dp,inf)/(norm(p0,inf) + 0.5*mm^2)]);
  waitbar(n/nstep)
end

disp(['Relative change = ',num2str(relchange(end))])
divergence = norm(Du*u1 + Dv*v1,inf);
if (geval==1)|(geval==2)
  errornorm		% Maximum norm of error in outflow profile if exact
			% solution is available
end
close(h)
tic
isobarplot		% Plot of isobars and velocity vectors
streamlineplot 		% Plot of streamlines and velocity vectors
velocity_profile	% Plot of velocity profiles
relchangeplot		% Plot of relative change per time step
tijd = toc; disp(['Plotting time = ',num2str(tijd)])
disp(['Total time = ',num2str(etime(clock,time0))])


text:
function colors = RainbowColors_ss(numcolors)
% Generates RGB Triplets based on the idealized colors below:
% 0. red    [1 0 0]
% 1. orange [1 .5 0]
% 2. yellow [1 1 0]
% 2.5. green  [0 1 0]
% 3. cyan   [0 1 1]
% 4. blue   [0 0 1]
% 5. purple [.5 0 .5]
%
% INPUT
% -numcolors: number of colors you want out
% OUTPUT
% - colors a 3 column matrix signifying Red Green and Blue values with each
% column, number of rows = numcolors
%
% Brendon Watson 2015

cidxs = linspace(0,1,numcolors);
r = redfunc(cidxs);
g = greenfunc(cidxs);
b = bluefunc(cidxs);

colors = [r g b];

1;

function y = redfunc(x)
for a = 1:length(x)
    if x(a)<2/5
        y(a) = 1;
    elseif x(a)>=2/5 && x(a)<2.5/5
        y(a) = 1-5*(x(a)-2/5);
    elseif x(a)>=2.5/5 && x(a)<4/5
        y(a) = 0;
    elseif x(a)>=4/5
        y(a) = 0+.5*5*(x(a)-4/5);
    end
end
y = y';

function y = greenfunc(x)
for a = 1:length(x)
    if x(a)<2/5
        y(a) = .5*5*(x(a));
    elseif x(a)>=2/5 && x(a)<3/5
        y(a) = 1;
    elseif x(a)>=3/5 && x(a)<4/5
        y(a) = 1-5*(x(a)-3/5);
    elseif x(a)>=4/5
        y(a) = 0;
    end
end
y = y';

function y = bluefunc(x)
for a = 1:length(x)
    if x(a)<3/5
        y(a) = 0;
    elseif x(a)>=2.5/5 && x(a)<3/5
        y(a) = 5*(x(a)-2.5/5);
    elseif x(a)>=3/5 && x(a)<4/5
        y(a) = 1;
    elseif x(a)>=4/5
        y(a) = 1-.5*5*(x(a)-4/5);
    end
end
y = y';

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode