text:
function db_rename_condition( oldPath, newPath, isMove, isUpdateStudyPath )
% DB_RENAME_CONDITION: Rename a condition in the Brainstorm database.
%
% USAGE:  db_rename_condition( oldPath, newPath, isMove=1, isUpdateStudyPath=1)
%         db_rename_condition( oldPath, newPath )
%
% INPUT:
%    - oldPath : Condition path to modify ("SubjectName/ConditionName", or "*/ConditionName")
%    - newPath : New condition path
%    - isMove  : If 1, actually rename the folders
%                Set to 0 when renaming a subject (moving entire subject folder at once)
%    - isUpdateStudyPath: If 1, modify the study path
%                         If 0, only update the file links

% @=============================================================================
% 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, 2009-2014

% Parse inputs
if (nargin < 4) || isempty(isUpdateStudyPath)
    isUpdateStudyPath = 1;
end
if (nargin < 3) || isempty(isMove)
    isMove = 1;
end
% No modification
if strcmpi(oldPath, newPath)
    return;
end
oldPath = strrep(oldPath, '\', '/');
newPath = strrep(newPath, '\', '/');

%% ===== MULTIPLE RENAMING =====
if any(oldPath == '*')
    % Get all the studies concerned by the modification
    [sStudies, iStudies] = bst_get('StudyWithCondition', oldPath);
    if isempty(sStudies)
        disp(['RENAME> Condition "' oldPath '" does not exist.']);
        return;
    end
    % Recursive calls
    for i = 1:length(sStudies)
        % Get subject/condition names
        [SubjectName, oldCond] = bst_fileparts(bst_fileparts(sStudies(i).FileName), 1);
        [tmp__, newCond] = bst_fileparts(newPath, 1);
        % Rename
        db_rename_condition(bst_fullfile(SubjectName,oldCond), bst_fullfile(SubjectName,newCond), isMove);
    end
    return
end
% Progress bar
bst_progress('start', 'Rename condition', ['Rename: ' oldPath ' => ' newPath]);

%% ===== PREPARING =====
% Get protocol directories
ProtocolInfo = bst_get('ProtocolInfo');
% Get subject/condition names
[oldSubj, oldCond] = bst_fileparts(oldPath, 1);
[newSubj, newCond] = bst_fileparts(newPath, 1);
% Check that old condition exists
if ~isdir(bst_fullfile(ProtocolInfo.STUDIES, oldPath))
    bst_error(['Folder "' oldPath '" does not exists.'], 'Rename', 0);
    return;
end
% Check that new condition does not exist
if isMove && isdir(bst_fullfile(ProtocolInfo.STUDIES, newPath))
    bst_error(['Folder "' newPath '" already exists.'], 'Rename', 0);
    return;
end

%% ===== RENAME FOLDER =====
% Move the folder
if isMove
    isOk = file_move(bst_fullfile(ProtocolInfo.STUDIES, oldPath), bst_fullfile(ProtocolInfo.STUDIES, newPath));
    if ~isOk
        bst_error(['Error: Could not rename "' oldPath '" to "' newPath '".'], 'Rename', 0);
        return;
    end
end

%% ===== UPDATE ALL THE FILES =====
% Get study structure
[sStudy, iStudy] = bst_get('StudyWithCondition', oldPath);
% Update study
if isUpdateStudyPath
    sStudy = replaceStruct(sStudy, 'FileName',        oldPath, newPath);
    sStudy = replaceStruct(sStudy, 'Name',            oldPath, newPath);
    sStudy = replaceStruct(sStudy, 'BrainStormSubject', oldSubj, newSubj);
    sStudy.Condition = {newCond};
end

% === CHANNEL ===
for i = 1:length(sStudy.Channel)
    sStudy.Channel(i) = replaceStruct(sStudy.Channel(i), 'FileName', oldPath, newPath);
end
% === DATA ===
for i = 1:length(sStudy.Data)
    sStudy.Data(i) = replaceStruct(sStudy.Data(i), 'FileName', oldPath, newPath);
end
% === HEADMODEL ===
for i = 1:length(sStudy.HeadModel)
    oldFile = sStudy.HeadModel(i).FileName;
    sStudy.HeadModel(i) = replaceStruct(sStudy.HeadModel(i), 'FileName', oldPath, newPath);
    % Update file
    if isMove
        fileFull = bst_fullfile(ProtocolInfo.STUDIES, sStudy.HeadModel(i).FileName);
    else
        fileFull = bst_fullfile(ProtocolInfo.STUDIES, oldFile);
    end
    if file_exist(fileFull)
        fileMat = load(fileFull);
        [fileMat, isModified] = replaceStruct(fileMat, 'SurfaceFile', oldSubj, newSubj);
        if isModified
            bst_save(fileFull, fileMat, 'v7');
        end
    end
end
% === RESULT ===
for i = 1:length(sStudy.Result)
    oldFile = sStudy.Result(i).FileName;
    sStudy.Result(i) = replaceStruct(sStudy.Result(i), 'FileName', oldPath, newPath);
    sStudy.Result(i) = replaceStruct(sStudy.Result(i), 'DataFile', oldPath, newPath);
    % Regular files
    if ~sStudy.Result(i).isLink
        % Update file
        if isMove
            fileFull = bst_fullfile(ProtocolInfo.STUDIES, sStudy.Result(i).FileName);
        else
            fileFull = bst_fullfile(ProtocolInfo.STUDIES, oldFile);
        end
        if file_exist(fileFull)
            fileMat = load(fileFull);
            [fileMat, isModified1] = replaceStruct(fileMat, 'DataFile',      oldPath, newPath);
            [fileMat, isModified2] = replaceStruct(fileMat, 'HeadModelFile', oldPath, newPath);
            [fileMat, isModified3] = replaceStruct(fileMat, 'SurfaceFile',   oldSubj, newSubj);
            if isModified1 || isModified2 || isModified3
                bst_save(fileFull, fileMat, 'v6');
            end
        end
    end
end
% === STAT ===
for i = 1:length(sStudy.Stat)
    oldFile = sStudy.Stat(i).FileName;
    sStudy.Stat(i) = replaceStruct(sStudy.Stat(i), 'FileName', oldPath, newPath);
    % Update file
    if isMove
        fileFull = bst_fullfile(ProtocolInfo.STUDIES, sStudy.Stat(i).FileName);
    else
        fileFull = bst_fullfile(ProtocolInfo.STUDIES, oldFile);
    end
    if file_exist(fileFull)
        fileMat = load(fileFull);
        [fileMat, isModified] = replaceStruct(fileMat, 'SurfaceFile', oldSubj, newSubj);
        if isModified
            bst_save(fileFull, fileMat, 'v6');
        end
    end
end
% === IMAGE ===
for i = 1:length(sStudy.Image)
    sStudy.Image(i) = replaceStruct(sStudy.Image(i), 'FileName', oldPath, newPath);
end
% === NOISECOV ===
for i = 1:length(sStudy.NoiseCov)
    sStudy.NoiseCov(i) = replaceStruct(sStudy.NoiseCov(i), 'FileName', oldPath, newPath);
end
% === DIPOLES ===
for i = 1:length(sStudy.Dipoles)
    oldFile = sStudy.Dipoles(i).FileName;
    sStudy.Dipoles(i) = replaceStruct(sStudy.Dipoles(i), 'FileName', oldPath, newPath);
    sStudy.Dipoles(i) = replaceStruct(sStudy.Dipoles(i), 'DataFile', oldPath, newPath);
    % Update file
    if isMove
        fileFull = bst_fullfile(ProtocolInfo.STUDIES, sStudy.Dipoles(i).FileName);
    else
        fileFull = bst_fullfile(ProtocolInfo.STUDIES, oldFile);
    end
    if file_exist(fileFull)
        fileMat = load(fileFull);
        [fileMat, isModified1] = replaceStruct(fileMat, 'DataFile', oldPath, newPath);
        [fileMat, isModified2] = replaceStruct(fileMat, 'SurfaceFile', oldPath, newPath);
        if isModified1 || isModified2
            bst_save(fileFull, fileMat, 'v7');
        end
    end
end
% === TIMEFREQ ===
for i = 1:length(sStudy.Timefreq)
    oldFile = sStudy.Timefreq(i).FileName;
    sStudy.Timefreq(i) = replaceStruct(sStudy.Timefreq(i), 'FileName', oldPath, newPath);
    sStudy.Timefreq(i) = replaceStruct(sStudy.Timefreq(i), 'DataFile', oldPath, newPath);
    % Update file
    if isMove
        fileFull = bst_fullfile(ProtocolInfo.STUDIES, sStudy.Timefreq(i).FileName);
    else
        fileFull = bst_fullfile(ProtocolInfo.STUDIES, oldFile);
    end
    if file_exist(fileFull)
        fileMat = load(fileFull);
        [fileMat, isModified1] = replaceStruct(fileMat, 'DataFile', oldPath, newPath);
        [fileMat, isModified2] = replaceStruct(fileMat, 'SurfaceFile', oldPath, newPath);
        [fileMat, isModified3] = replaceStruct(fileMat, 'HeadModelFile', oldPath, newPath);
        if isModified1 || isModified2 || isModified3
            bst_save(fileFull, fileMat, 'v6');
        end
    end
end
% === MATRIX ===
for i = 1:length(sStudy.Matrix)
    sStudy.Matrix(i) = replaceStruct(sStudy.Matrix(i), 'FileName', oldPath, newPath);
end

%% ===== UPDATE DATABASE =====
% Update condition in database
bst_set('Study', iStudy, sStudy);
% Close progress bar
bst_progress('stop');
% Not moving, not reloading
if isMove
    % Update tree display
    panel_protocols('UpdateTree');
end

end

%% ===== HELPER FUNCTIONS =====
function [s, isModified] = replaceStruct(s, field, oldPath, newPath)
    isModified = 0;
    if isempty(s) || ~isfield(s, field) || isempty(s.(field))
        return
    end
    % Get new and old subject name
    oldSubj = bst_fileparts(oldPath);
    newSubj = bst_fileparts(newPath);
    % Link: two filenames to modify
    if (length(s.(field)) > 4) && strcmpi(s.(field)(1:4), 'link')
        % Split the link
        splitLink = str_split(s.(field), '|');
        % Replace two strings
        for i = 2:3
            % Get file path
            [fPath,fBase,fExt] = bst_fileparts(splitLink{i});
            % Replace "oldPath" with "newPath"
            if ~isempty(strfind(file_win2unix(fPath), file_win2unix(oldPath)))
                fPath = strrep(file_win2unix(fPath), file_win2unix(oldPath), file_win2unix(newPath));
            % If the subject was renamed: change the subject directly (in the case of link files to @default_study)
            elseif ~strcmpi(oldSubj, newSubj)
                fPath = strrep(file_win2unix(fPath), [file_win2unix(oldSubj) '/'], [file_win2unix(newSubj) '/']);
            end
            % Update file path
            splitLink{i} = bst_fullfile(fPath, [fBase, fExt]);
        end
        % Rebuild full link
        s.(field) = ['link|' file_win2unix(splitLink{2}) '|' file_win2unix(splitLink{3})];
    % Regular filename
    else
        % Get file path
        [fPath,fBase,fExt] = bst_fileparts(s.(field));
        % Replace "oldPath" with "newPath"
        if ~isempty(strfind(file_win2unix(fPath), file_win2unix(oldPath)))
            fPath = strrep(file_win2unix(fPath), file_win2unix(oldPath), file_win2unix(newPath));
        % If the subject was renamed: change the subject directly (in the case of link files to @default_study)
        elseif ~strcmpi(oldSubj, newSubj)
            fPath = strrep([file_win2unix(fPath) '/'], [file_win2unix(oldSubj) '/'], [file_win2unix(newSubj) '/']);
            fPath = fPath(1:end-1);
        end
        % Update file path
        s.(field) = bst_fullfile(fPath, [fBase, fExt]);
        s.(field) = file_win2unix(s.(field));
    end
    isModified = 1;
end

author: brainstorm-tools


text:
function [ v, d ] = pts_pca( x, kernel )
% PTS_PCA  Linear and Kernel Principal Component Analysis (PCA and KPCA)
%
% [ V, D ] = PTS_PCA(X)
%
%   Center X and compute eigenvectors V and eigenvalues D using linear PCA.
%
%   X is a matrix where each column is a vector. Each row is a sample of a
%   variable.
%
%   The eigenproblem solved is
%
%     S * v = d * v
%
%   where v is each eigenvector, and d each corresponding eigenvalue. S is
%   the biased covariance matrix
%
%     S = 1/M * Xc * Xc'
%
%   where Xc is the matrix of centered vectors X.
%
%   NOTE: If there are less vectors than variables, or kernel PCA is
%   selected, MDS is used to speed up computation of eigenvectors.
%
%   NOTE: Tiny eigenvalues can give numerical errors, so eigenvalues <
%   1e-13 and the corresponding eigenvectors are removed from the output.
%   Negative eigenvalues are assumed to arise from numerical errors, and
%   removed too.
%
%   NOTE: Matlab's function cov() computes the _unbiased_ covariance.
%
% [ A, DA ] = PTS_PCA(X, KERNEL)
%
%   Generalization to any type of kernel. KERNEL is a struct with the
%   kernel description
%
%     Linear PCA, as above
%
%       KERNEL.type = 'linear'
%
%     Kernel PCA, Polynomial homogeneous order q
%
%       KERNEL.type = 'polyh'
%       KERNEL.params = q
%
%     Kernel PCA, Gaussian with standard deviation sigma
%
%       KERNEL.type = 'gauss'
%       KERNEL.params = sigma
%
%   Because in general we want to avoid any explicit computations in
%   feature space, we solve the eigenproblem [1]
%
%     Kc * a = da * a
%
%   where a is each "coefficient eigenvector", and da = m * d, where m is
%   the number of vectors X. Kc is the kernel matrix of centered feature
%   vectors phi(X).
%
%   If you want to obtain the eigenvectors and eigenvalues in feature space
%   (for those cases where feature space is finite dimensional), you need
%   to compute the feature space training vectors phi(X), center them to
%   phic, and do
%
%     v = phic * a * diag(1./sqrt(da))
%     d = da / m
%
%   where m is the number of training vectors. (Note that v computed as
%   above is equal to v computed from the feature vectors save possibly the
%   sign).
%
% [1] B. Schölkopf, A.J. Smola and K.-R. Müller. "Kernel Principal
% Component Analysis". In "Advances in Kernel Methods - SV Learning", B.
% Schölkopf, C.J.C. Burges and A.J. Smola, eds. pp. 327--352, MIT Press,
% 1999.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2009-2011 University of Oxford
% Version: 1.1.0
% $Rev: 434 $
% $Date: 2011-06-03 17:16:33 +0100 (Fri, 03 Jun 2011) $
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error( nargchk( 1, 2, nargin, 'struct' ) );
error( nargoutchk( 0, 2, nargout, 'struct' ) );

if ndims( x ) > 2
    error( 'X must be a matrix, not a volume' )
end

% get sizes
%   N: vector length
%   M: number of training vectors
[ N, M ] = size( x );

% defaults
if ( nargin < 2 || isempty( kernel ) )
    kernel.type = 'linear';
end

switch kernel.type

    case 'linear'

        % center training vectors
        xmean = mean( x, 2 );
        x = x - repmat( xmean, 1, M );

        % compute matrix for eigenproblem
        if ( M < N ) % "kernel matrix" trick to speed up computation
            x_cov = x' * x;
        else % normal case, covariance matrix
            x_cov = x * x' / M;
        end

        % compute eigenvalues and eigenvectors
        [ v, d ] = eig( x_cov );

        % save memory by reducing eigenvalues to a vector
        d = diag( d );

        % reoder in decreasing order of signed modulus value
        [~, idx] = sort(abs(d).*sign(d), 1, 'descend');
        d = d(idx);
        v = v(:, idx);

        % if "kernel matrix" trick was used, it is necessary to convert the
        % coefficient eigenvectors and eigenvalues to feature space
        % eigenvectors and eigenvalues
        if ( M < N )

            % scale eigenvectors
            v = x * ( v * diag( 1 ./ sqrt( d ) ) );

            % scale eigenvalues
            d = d / M;

        end

    otherwise

        % for kernel PCA, we use a generalization of the speed up trick
        % above
        k = pts_kmat( kernel, x );

        % center kernel matrix (this is the same as centering the feature
        % space training vectors and then computing the kernel matrix)
        onesm = ones( M ) / M;
        k = k - onesm * k - k * onesm + onesm * k * onesm;

        % k must be symmetric in order to avoid complex eigenvalues, and in
        % theory it should be, but there are small errors due to numeric
        % precision when the centered kernel matrix is computed
        k = ( k + k' ) / 2;

        % compute eigenvalues and eigenvectors
        [ v, d ] = eig( k );

        % save memory by reducing eigenvalues to a vector
        d = diag( d );

        % reoder in decreasing order of signed modulus value
        [~, idx] = sort(abs(d).*sign(d), 1, 'descend');
        d = d(idx);
        v = v(:, idx);

end

author: tomdoel

repo: pulmonarytoolkit


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Q = INVERSEKINEMATIC_IRB7600_400_255_m2000(robot, T)
%   Solves the inverse kinematic problem for the ABB IRB7600_400_255_m2000 robot
%   where:
%   robot stores the robot parameters.
%   T is an homogeneous transform that specifies the position/orientation
%   of the end effector.
%
%   A call to Q=INVERSEKINEMATIC_IRB7600_400_255_m2000 returns 8 possible solutions, thus,
%   Q is a 6x8 matrix where each column stores 6 feasible joint values.
%
%
%   Example code:
%
%   abb=load_robot('ABB', 'IRB7600_400_255_m2000');
%   q = [0 0 0 0 0 0];
%   T = directkinematic(abb, q);
%   %Call the inversekinematic for this robot
%   qinv = inversekinematic(abb, T);
%   check that all of them are feasible solutions!
%   and every Ti equals T
%   for i=1:8,
%        Ti = directkinematic(abb, qinv(:,i))
%   end
%	See also DIRECTKINEMATIC.
%
%   Author: Arturo Gil Aparicio
%           Universitas Miguel Hernandez, SPAIN.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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 q = inversekinematic_irb7600_400_255_m2000(robot, T)

%initialize q,
%eight possible solutions are generally feasible
q=zeros(6,8);

% %Evaluate the parameters
% theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
L6=d(6);

%T= [ nx ox ax Px;
%     ny oy ay Py;
%     nz oz az Pz];
Px=T(1,4);
Py=T(2,4);
Pz=T(3,4);

%Compute the position of the wrist, being W the Z component of the end effector's system
W = T(1:3,3);

% Pm: wrist position
Pm = [Px Py Pz]' - L6*W;

%first joint, two possible solutions admited:
% if q(1) is a solution, then q(1) + pi is also a solution
q1=atan2(Pm(2), Pm(1));

%solve for q2
q2_1=solve_for_theta2(robot, [q1 0 0 0 0 0 0], Pm);
%the other possible solution is q1 + pi
q2_2=solve_for_theta2(robot, [q1+pi 0 0 0 0 0 0], Pm);

%solve for q3
q3_1=solve_for_theta3(robot, [q1 0 0 0 0 0 0], Pm);
%solver for q3 for both cases
q3_2=solve_for_theta3(robot, [q1+pi 0 0 0 0 0 0], Pm);

%the next matrix doubles each column. For each two columns, two different
%configurations for theta4, theta5 and theta6 will be computed. These
%configurations are generally referred as wrist up and wrist down solution
q = [q1         q1         q1        q1       q1+pi   q1+pi   q1+pi   q1+pi;
     q2_1(1)    q2_1(1)    q2_1(2)   q2_1(2)  q2_2(1) q2_2(1) q2_2(2) q2_2(2);
     q3_1(1)    q3_1(1)    q3_1(2)   q3_1(2)  q3_2(1) q3_2(1) q3_2(2) q3_2(2);
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0;
     0          0          0         0         0      0       0       0];

 %leave only the real part of the solutions
q=real(q);

%Note that in this robot, the joint q3 has a non-simmetrical range. In this
%case, the joint ranges from 60 deg to -219 deg, thus, the typical normalizing
%step is avoided in this angle (the next line is commented). When solving
%for the orientation, the solutions are normalized to the [-pi, pi] range
%only for the theta4, theta5 and theta6 joints.

%normalize q to [-pi, pi]
q(1,:) = normalize(q(1,:));
q(2,:) = normalize(q(2,:));

% solve for the last three joints
% for any of the possible combinations (theta1, theta2, theta3)
for i=1:2:size(q,2),
    % use solve_spherical_wrist2 for the particular orientation
    % of the systems in this ABB robot
    % use either the geometric or algebraic method.
    % the function solve_spherical_wrist2 is used due to the relative
    % orientation of the last three DH reference systems.

    %use either one algebraic method or the geometric
    %qtemp = solve_spherical_wrist2(robot, q(:,i), T, 1, 'geometric'); %wrist up
    qtemp = solve_spherical_wrist2(robot, q(:,i), T, 1,'algebraic'); %wrist up
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i)=qtemp;

    %qtemp = solve_spherical_wrist2(robot, q(:,i), T, -1, 'geometric'); %wrist down
    qtemp = solve_spherical_wrist2(robot, q(:,i), T, -1, 'algebraic'); %wrist down
    qtemp(4:6)=normalize(qtemp(4:6));
    q(:,i+1)=qtemp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for second joint theta2, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q2 = solve_for_theta2(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=a(2);
L3=d(4);
A2=a(3);

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%The inverse kinematic problem can be solved as in the IRB 140 (for example)

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

beta = atan2(-p1(2), p1(1));
gamma = real(acos((L2^2+r^2-L4^2)/(2*r*L2)));

%return two possible solutions
%elbow up and elbow down
%the order here is important and is coordinated with the function
%solve_for_theta3
q2(1) = pi/2 - beta - gamma; %elbow up
q2(2) = pi/2 - beta + gamma; %elbow down

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve for third joint theta3, two different
% solutions are returned, corresponding
% to elbow up and down solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function q3 = solve_for_theta3(robot, q, Pm)

%Evaluate the parameters
theta = eval(robot.DH.theta);
d = eval(robot.DH.d);
a = eval(robot.DH.a);
alpha = eval(robot.DH.alpha);

%See geometry
L2=a(2);
L3=d(4);
A2=a(3);

%See geometry of the robot
%compute L4
L4 = sqrt(A2^2 + L3^2);

%the angle phi is fixed
phi=acos((A2^2+L4^2-L3^2)/(2*A2*L4));

%given q1 is known, compute first DH transformation
T01=dh(robot, q, 1);

%Express Pm in the reference system 1, for convenience
p1 = inv(T01)*[Pm; 1];

r = sqrt(p1(1)^2 + p1(2)^2);

eta = real(acos((L2^2 + L4^2 - r^2)/(2*L2*L4)));

%return two possible solutions
%elbow up and elbow down solutions
%the order here is important
q3(1) = pi - phi- eta;
q3(2) = pi - phi + eta;

repo: ARTE


text:
classdef UniversalElementSetVariable < AbstractOrbitModelVariable
    %UniversalElementSetVariable Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) UniversalElementSet = UniversalElementSet.getDefaultElements();

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

        varC3(1,1) logical = false;
        varRp(1,1) logical = false;
        varInc(1,1) logical = false;
        varRaan(1,1) logical = false;
        varArg(1,1) logical = false;
        varTau(1,1) logical = false;
    end

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

            obj.id = rand();
        end

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

            if(obj.varC3)
                x(end+1) = obj.varObj.c3;
            end

            if(obj.varRp)
                x(end+1) = obj.varObj.rP;
            end

            if(obj.varInc)
                x(end+1) = obj.varObj.inc;
            end

            if(obj.varRaan)
                x(end+1) = obj.varObj.raan;
            end

            if(obj.varArg)
                x(end+1) = obj.varObj.arg;
            end

            if(obj.varTau)
                x(end+1) = obj.varObj.tau;
            end
        end

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

            lb = obj.lb(useTf);
            ub = obj.ub(useTf);
        end

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

        function setBndsForVariable(obj, lb, ub)
            if(length(lb) == 6 && length(ub) == 6)
                obj.lb = lb;
                obj.ub = ub;
            else
                useTf = obj.getUseTfForVariable();

                obj.lb(useTf) = lb;
                obj.ub(useTf) = ub;
            end
        end

        function useTf = getUseTfForVariable(obj)
            useTf = [obj.varC3        obj.varRp        obj.varInc        obj.varRaan        obj.varArg        obj.varTau];
        end

        function setUseTfForVariable(obj, useTf)
            obj.varC3  = useTf(1);
            obj.varRp = useTf(2);
            obj.varInc  = useTf(3);
            obj.varRaan = useTf(4);
            obj.varArg  = useTf(5);
            obj.varTau = useTf(6);
        end

        function updateObjWithVarValue(obj, x)
            xInd = 1;

            if(obj.varC3)
                obj.varObj.c3 = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varRp)
                obj.varObj.rP = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varInc)
                obj.varObj.inc = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varRaan)
                obj.varObj.raan = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varArg)
                obj.varObj.arg = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varTau)
                obj.varObj.tau = x(xInd);
                xInd = xInd + 1; %#ok<NASGU>
            end
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
            if(evtNum > 0)
                subStr = sprintf('Event %i',evtNum);
            else
                subStr = varLocType;
            end

            nameStrs = {sprintf('%s C3 Energy', subStr), ...
                        sprintf('%s Radius of Periapsis', subStr), ...
                        sprintf('%s Inclination', subStr), ...
                        sprintf('%s RAAN', subStr), ...
                        sprintf('%s Argument of Periapsis', subStr), ...
                        sprintf('%s Time Past Periapsis', subStr)};

            nameStrs = nameStrs(obj.getUseTfForVariable());
        end

        function varsStoredInRad = getVarsStoredInRad(obj)
            varsStoredInRad = [false false true true true false];
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%This is a way to do a circular cross correlation of the columns of an
%array.
%
%In this file we take the FFT of the array, the FFT works on the columns
%only, then we make a second array, which is the conjugate of the first, we
%multiply these 2 arrays in an element by element manner, then do the
%inverse FFT.
%After this we circularly shift the columns of the conjugate array, go
%through the process again, then put the results nextto the previous ones.
%
%The first results will be autocorrelation, the subsequent ones will be
%circular cross correlation
%
%This goes pretty quickly if the output array is preallocated, so the method
%was developed to allow preallocation.
%
%The input variables are:
%NRows => The number of Rows of the input array
%NColumns=> guess
%Signal => the matrix with the time data in columns, for instance if you
%       are looking at the cross correlation of orthogonal codes, each code
%       stream would be a column, of course this is NRows by NColumns
%
%You can change this m-file to a function if you feel like it

CircXCorrData=zeros(NRows,NColumns^2);%Preallocate space
h = waitbar(0,'Please wait...');
SignalFFT=fft(Signal);
ConjSignalFFT=conj(SignalFFT);
for index= 1 : NColumns : ( (NColumns-1)*NColumns+1 )
    waitbar(index/( (NColumns-1)*NColumns+1 ))
    CircXCorrData( :,index : index+(NColumns-1) )=ifft(SignalFFT.*ConjSignalFFT);
    ConjSignalFFT=circshift(ConjSignalFFT,[0,1]);
end
close(h)

%Comment the plot out if you want
plot(max(abs(CircXCorrData)))
grid on

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function Y = as_l(tr,X)
%
%  Solves linear systems with the real, symmetric, negative definite matrix A,
%  i.e., Y = inv(A)*X.
%
%  The Cholesky factor of -A is provided as global data. This data must
%  be generated by calling 'as_l_i' before calling this routine!
%
%  Calling sequence:
%
%    Y = as_l(tr,X)
%
%  Input:
%
%    tr        is not referenced;
%    X         matrix of proper size.
%
%  Output:
%
%    Y         the solution matrix.
%
%
%   LYAPACK 1.0 (Thilo Penzl, May 1999)

if nargin~=2
  error('Wrong number of input arguments.');
end

global LP_U

if ~length(LP_U)
  error('This routine needs global data which must be generated by calling ''as_l_i'' first.');
end

Y = -LP_U\(LP_U'\X);      % Note the minus!

author: Sable

repo: mcbench-benchmarks


text:
% LegendrePoly.m by David Terr, Raytheon, 5-10-04

% Given nonnegative integer n, compute the
% Legendre polynomial P_n. Return the result as a vector whose mth
% element is the coefficient of x^(n+1-m).
% polyval(LegendrePoly(n),x) evaluates P_n(x).

function pk = LegendrePoly(n)

if n==0
    pk = 1;
elseif n==1
    pk = [1 0]';
else

    pkm2 = zeros(n+1,1);
    pkm2(n+1) = 1;
    pkm1 = zeros(n+1,1);
    pkm1(n) = 1;

    for k=2:n

        pk = zeros(n+1,1);

        for e=n-k+1:2:n
            pk(e) = (2*k-1)*pkm1(e+1) + (1-k)*pkm2(e);
        end

        pk(n+1) = pk(n+1) + (1-k)*pkm2(n+1);
        pk = pk/k;

        if k<n
            pkm2 = pkm1;
            pkm1 = pk;
        end

    end

end


text:
function cs = load(fname,varargin)
% import crystal symmetry from phl or cif files
%
% if cif file not found and input name is a valid COD entry, this function
% tries to download the file from
% <http://www.crystallography.net/cif/ http://www.crystallography.net/cif/>
%
% Syntax
%   cs = crystalSymmetry.load('5000035.cif') % load from cif file
%   cs = crystalSymmetry.load('copper.phl') % load from phl file
%   cs = crystalSymmetry.load(5000035) % lookup online
%
% See also
% crystalSymmetry/crystalSymmetry
%

[~, ~, ext] = fileparts(char(fname));

if strcmpi(ext,'.phl')
  cs = loadPHL(fname,varargin{:});
else
  cs = loadCIF(fname,varargin{:});
end

author: mtex-toolbox

repo: mtex

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


text:
function n = normal(c, unit)
%NORMAL   Normal to a complex-valued CHEBFUN.
%   N = NORMAL(C) returns the normal vector to the curve C as a CHEBFUN with two
%   columns. The vector has the same magntiude as the curve's tangent vector.
%
%   N = NORMAL(C, 'unit') returns the unit normal vector to the curve C.

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

% [TODO]:  Reconsider how this function should behave in the presence of cusps
% once we have singfun in place.

n = -1i*diff(c);

if ( nargin > 1 )
    if ( strcmpi(unit, 'unit') )
        nrmn = norm(n);
        if ( nrmn == 0 )
            error('CHEBFUN:CHEBFUN:normal:zero', 'Normal vector is zero.');
        else
            n = n./nrmn;
        end
    else
        error('CHEBFUN:CHEBFUN:normal:args', ...
            'Second argument is not recognised.');
    end
end

% Return an array-valued CHEBFUN:
n = [real(n), imag(n)];

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% RSPFUNC - sub-function used by RSGET
%
% Usage: res = rspfunc(pval, l, rval)
%
% Input:
%   pval - p-value to optimize
%   l    - [l1 l2 l3 l4] l-values for Ramberg-Schmeiser distribution
%   rval - expected r-value
%
% Output:
%   res  - residual
%
% Author: Arnaud Delorme, SCCN, 2003
%
% See also: RSGET
%
% Reference: Ramberg, J.S., Tadikamalla, P.R., Dudewicz E.J., Mykkytka, E.F.
%            A probability distribution and its uses in fitting data.
%            Technimetrics, 1979, 21: 201-214.

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

function rp = rspfunc( pval, l, rval);

    % for fiting rp with fminsearch
    % -----------------------------
    rp   = l(1) + (pval.^l(3) - (1-pval).^l(4))/l(2);
    rp   = abs(rval-rp);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function cuboidhyp = cuboid_pts2hyp(pts, vp)

junc3dummy = struct('type',[],'pt',[],'regionid',[]);
cuboidhyp.junc3(8) = junc3dummy; % initialize

for i = 1:8
    cuboidhyp.junc3(i).type = i;
    cuboidhyp.junc3(i).pt = pts(i,:);
    cuboidhyp.junc3(i).regionid = getregionid(pts(i,1), pts(i,2), vp);
end

author: zouchuhang

repo: LayoutNet


text:
function ori = circle3dOrigin(varargin)
%CIRCLE3DORIGIN Return the first point of a 3D circle
%
%   P = circle3dOrigin([XC YC ZC R THETA PHI])
%   P = circle3dOrigin([XC YC ZC R THETA PHI PSI])
%   Returns the origin point of the circle, i.e. the first point used for
%   drawing circle.
%
%   See also:
%   circles3d, points3d, circle3dPosition
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 21/02/2005
%

%   HISTORY

% get center and radius
circle = varargin{1};
xc = circle(:,1);
yc = circle(:,2);
zc = circle(:,3);
r  = circle(:,4);

% get angle of normal
theta   = circle(:,5);
phi     = circle(:,6);

% get roll
if size(circle, 2)==7
    psi = circle(:,7);
else
    psi = zeros(size(circle, 1), 1);
end

% create origin point
pt0 = [r 0 0];

% compute transformation from local basis to world basis
trans   = localToGlobal3d(xc, yc, zc, theta, phi, psi);

% transform the point
ori = transformPoint3d(pt0, trans);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function order = level_to_order_open ( dim_num, level )

%*****************************************************************************80
%
%% LEVEL_TO_ORDER converts a level to an order for open rules.
%
%  Discussion:
%
%    Sparse grids can naturally be nested.  A natural scheme is to use
%    a series of one-dimensional rules arranged in a series of "levels"
%    whose order roughly doubles with each step.
%
%    The arrangement described here works naturally for the Fejer Type 1,
%    Fejer Type 2, Newton Cotes Open, Newton Cotes Half Open,
%    and Gauss-Patterson rules.  It also can be used, partially, to describe
%    the growth of Gauss-Legendre rules.
%
%    The idea is that we start with LEVEL = 0, ORDER = 1 indicating the single
%    point at the center, and for all values afterwards, we use the relationship
%
%      ORDER = 2**(LEVEL+1) - 1.
%
%    The following table shows how the growth will occur:
%
%    Level    Order
%
%    0          1
%    1          3 =  4 - 1
%    2          7 =  8 - 1
%    3         15 = 16 - 1
%    4         31 = 32 - 1
%    5         63 = 64 - 1
%
%    For the Fejer Type 1, Fejer Type 2, Newton Cotes Open,
%    Newton Cotes Open Half, and Gauss-Patterson rules, the point growth is
%    nested.  If we have ORDER points on a particular LEVEL, the next level
%    includes all these old points, plus ORDER+1 new points, formed in the
%    gaps between successive pairs of old points plus an extra point at each
%    end.
%
%    Level    Order = New + Old
%
%    0          1   =  1  +  0
%    1          3   =  2  +  1
%    2          7   =  4  +  3
%    3         15   =  8  +  7
%    4         31   = 16  + 15
%    5         63   = 32  + 31
%
%    If we use a series of Gauss-Legendre rules, then there is almost no
%    nesting, except that the central point is shared.  If we insist on
%    producing a comparable series of such points, then the "nesting" behavior
%    is as follows:
%
%    Level    Order = New + Old
%
%    0          1   =  1  +  0
%    1          3   =  2  +  1
%    2          7   =  6  +  1
%    3         15   = 14  +  1
%    4         31   = 30  +  1
%    5         63   = 62  +  1
%
%    Moreover, if we consider ALL the points used in such a set of "nested"
%    Gauss-Legendre rules, then we must sum the "NEW" column, and we see that
%    we get roughly twice as many points as for the truly nested rules.
%
%    In this routine, we assume that a vector of levels is given,
%    and the corresponding orders are desired.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL(DIM_NUM), the nesting level.
%
%    Output, integer ORDER(DIM_NUM), the order (number of points) of the rule.
%
  for dim = 1 : dim_num

    if ( level(dim) < 0 )
      order(dim) = -1;
    elseif ( level(dim) == 0 )
      order(dim) = 1;
    else
      order(dim) = 2^( level(dim) + 1 ) - 1;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% runstats.m
% This script collects statistics about positions of each component,
% creates matrix S 4xM [xmin_0, ..., xmin_{M-1}; ymin_0, ..., ymin_{M-1};
% xmin_0, ..., xmax_{M-1}; ymax_0, ..., ymax_{M-1}] each column defines
% bounding box for one component.
%
% Parameters in structure options are now set to cover approx. 97.5% of all
% tested images (about 12200 images, original database was pruned - file
% ./MAT/eye_on_lfw_pruned.mat - so now it does not contain images with
% multiple faces annotated or detected. Also images with incomplete
% annotations were discarded.)
%
% 15-07-10 Michal Uricar
% 26-02-11 Michal Uricar

% clc;
clearvars; close all;

%% Timestamp

fprintf(1,'Started on %s\n\n', datestr(now));

%% Load data (images with detected frames)

addpath('./Functions');
addpath('../../matlab_toolbox/mex/');

load('./MAT/paths.mat');
load('./MAT/options.mat');

% load('./MAT/eyefd_on_lfw_pruned.mat');
annotation_struct = read_lfw_annotation_file(['../' DB_ANN]);

%% Print parameters

fprintf('Parameters: \nbw = [%d; %d]\n', options.bw);
fprintf('bw_margin = [%d; %d]\n', options.bw_margin);
fprintf('components= \n'); disp(options.components);
fprintf('M = %d\n', options.M);
fprintf('image_path = ''%s''\n', options.image_path);
fprintf('....................................\n');
fprintf('\nStarting test...\n');

%% run test

tic
[cnt, S, bad_idx] = paramStats(annotation_struct, options);
% [cnt, S, bad_idx] = paramStats(annotation_struct, options, 1);
toc

percent = cnt/annotation_struct.N * 100;
fprintf('Passed %d images. That is %f%%.\n', cnt, percent);

%% create new database - only images that passed through test

idx = 1:annotation_struct.N;
idx(bad_idx) = [];
% image = image(idx);
annotation_struct.bbox = annotation_struct.bbox(idx, :);
annotation_struct.names = annotation_struct.names(idx);
annotation_struct.eye_r = annotation_struct.eye_r(idx, :);
annotation_struct.eye_l = annotation_struct.eye_l(idx, :);
annotation_struct.canthus_rr = annotation_struct.canthus_rr(idx, :);
annotation_struct.canthus_rl = annotation_struct.canthus_rl(idx, :);
annotation_struct.canthus_lr = annotation_struct.canthus_lr(idx, :);
annotation_struct.canthus_ll = annotation_struct.canthus_ll(idx, :);
annotation_struct.mouth = annotation_struct.mouth(idx, :);
annotation_struct.mouth_corner_r = annotation_struct.mouth_corner_r(idx, :);
annotation_struct.mouth_corner_l = annotation_struct.mouth_corner_l(idx, :);
annotation_struct.nose = annotation_struct.nose(idx, :);
annotation_struct.N = numel(annotation_struct.names);

fprintf('Creating new database (db_good.mat)...\n');
save('./MAT/db_good.mat', 'annotation_struct');

%% Save data

fprintf('Saving statistics...\n');
save('./results/statistics_orig.mat', 'cnt', 'S', 'percent', 'options');

 %% shrink S (S boundary must represent positions of components center)

load('./results/statistics_orig.mat', 'S');

for i = 1 : options.M
    S(1, i) = S(1, i) + options.components(1, i)/2;
    S(2, i) = S(2, i) + options.components(2, i)/2;
    S(3, i) = S(3, i) - options.components(1, i)/2;
    S(4, i) = S(4, i) - options.components(2, i)/2;
end;

save('./results/statistics_orig.mat', 'cnt', 'S', 'percent', 'options');

%% Show computed S

[Iframe, Annotation] = getImageFrame(options, 1, annotation_struct);
Points = Annotation.P;
Points = prepareS0gt(Annotation.P, options);    % transform nose to the center of face
Points = Points(:, options.comselect);          % extract relevant points only (name list in options.compnames)
Points(:, options.M) = Annotation.P(:, 10);             % copy back original nose position

colors = colormap(hsv(options.M)); close gcf;
% names = fieldnames(Annotation.face);
names = options.compnames;

scrsz = get(0,'ScreenSize');

figure;
imshow(Iframe, [], 'Border', 'tight'); hold on;
set(gcf, 'OuterPosition', [scrsz(1) scrsz(2) scrsz(3) scrsz(4)]);
% set(gcf, 'Position', [100, 100, 800, 600]);
for k = 1 : options.M
%     plot(Points(1, k), Points(2, k), 'r.', 'LineWidth', 2, 'MarkerSize', 10);
    plot(Points(1, k), Points(2, k), '.', 'color', colors(k, :), 'LineWidth', 2, 'MarkerSize', 10);
    bb = makeAABB(S(:, k));
    line([bb(1,:) bb(1,1)], [bb(2,:) bb(2,1)], 'color', colors(k, :));
    bb = makeAABB([Points(1, k) - options.components(1, k)/2 Points(2, k) - options.components(2, k)/2 ...
                   Points(1, k) + options.components(1, k)/2 Points(2, k) + options.components(2, k)/2 ]);
%     text(Points(1, k)-1, Points(2, k)+1, names(k+3), 'color', colors(k, :));
    text(Points(1, k)-1, Points(2, k)+1, names(k), 'color', colors(k, :));
    line([bb(1,:) bb(1,1)], [bb(2,:) bb(2,1)], 'color', 'y');
end;
hold off;

if (~exist('./img/', 'dir'))
    mkdir('./img');
end;

saveas(gcf, './img/S.png');

%% Show S

S = ceil(S);

colors = colormap(hsv(options.M)); close gcf;
% names = fieldnames(Annotation.face);
names = options.compnames;
i = 1;

figure;
subplot(1, 2, 1);
imshow(Iframe, []); hold on;
for j = 1 : options.M;
    bb = makeAABB(S(:, j));
    line([bb(1,:) bb(1,1)], [bb(2,:) bb(2,1)], 'color', colors(j, :));
    text(Points(1, j), Points(2, j), names(j), 'color', colors(j, :));
end;
line([20 20], [0 41], 'color', 'm', 'LineWidth', 2);
hold off;
saveas(gcf, './img/S.png');

%% make each frame of component in S symetric (this works only on lfw database with 4 components!)

S_new = S;

% canthus - inner
S_new(2, 2) = min(S(2, 2), S(2, 3)); S_new(2, 3) = min(S(2, 2), S(2, 3));   % ymin
S_new(4, 2) = max(S(4, 2), S(4, 3)); S_new(4, 3) = max(S(4, 2), S(4, 3));   % ymax
S_new(3, 2) = options.bw(1)/2 + max(S(3, 2) - options.bw(1)/2, options.bw(1)/2 - S(1, 3));     % x_max
S_new(1, 3) = options.bw(1)/2 - max(S(3, 2) - options.bw(1)/2, options.bw(1)/2 - S(1, 3));
S_new(1, 2) = min(S(1, 2), options.bw(1) - S(3, 3));
S_new(3, 3) = options.bw(1) - min(S(1, 2), options.bw(1) - S(3, 3));
% canthus - outer
S_new(2, 6) = min(S(2, 6), S(2, 7)); S_new(2, 7) = min(S(2, 6), S(2, 7));   % ymin
S_new(4, 6) = max(S(4, 6), S(4, 7)); S_new(4, 7) = max(S(4, 6), S(4, 7));   % ymax
S_new(3, 6) = options.bw(1)/2 + max(S(3, 6) - options.bw(1)/2, options.bw(1)/2 - S(1, 7));     % x_max
S_new(1, 7) = options.bw(1)/2 - max(S(3, 6) - options.bw(1)/2, options.bw(1)/2 - S(1, 7));
S_new(1, 6) = min(S(1, 6), options.bw(1) - S(3, 7));
S_new(3, 7) = options.bw(1) - min(S(1, 6), options.bw(1) - S(3, 7));
% mouth corners
S_new(2, 4) = min(S(2, 4), S(2, 5)); S_new(2, 5) = min(S(2, 4), S(2, 5));   % ymin
S_new(4, 4) = max(S(4, 4), S(4, 5)); S_new(4, 5) = max(S(4, 4), S(4, 5));   % ymax
S_new(3, 4) = options.bw(1)/2 + max(S(3, 4) - options.bw(1)/2, options.bw(1)/2 - S(1, 5));     % x_max
S_new(1, 5) = options.bw(1)/2 - max(S(3, 4) - options.bw(1)/2, options.bw(1)/2 - S(1, 5));
S_new(1, 4) = min(S(1, 4), options.bw(1) - S(3, 5));
S_new(3, 5) = options.bw(1) - min(S(1, 4), options.bw(1) - S(3, 5));
% center of face
S_new(3, 1) = options.bw(1)/2 + max(S(3, 1) - options.bw(1)/2, options.bw(1)/2 - S(1, 1));     % x_max
S_new(1, 1) = options.bw(1)/2 - max(S(3, 1) - options.bw(1)/2, options.bw(1)/2 - S(1, 1));
% nose
S_new(3, 8) = options.bw(1)/2 + max(S(3, 8) - options.bw(1)/2, options.bw(1)/2 - S(1, 8));     % x_max
S_new(1, 8) = options.bw(1)/2 - max(S(3, 8) - options.bw(1)/2, options.bw(1)/2 - S(1, 8));

subplot(1, 2, 2);
imshow(Iframe); hold on;
for j = 1 : options.M;
    bb = makeAABB(S_new(:, j));
    line([bb(1,:) bb(1,1)], [bb(2,:) bb(2,1)], 'color', colors(j, :));
    text(Points(1, j), Points(2, j), names(j), 'color', colors(j, :));
end;
line([20 20], [0 41], 'color', 'm', 'LineWidth', 2);
hold off;
saveas(gcf, './img/S_symetric.png');

%% Save symetric S for 4 components
S2 = S_new;
% save('./results/symetric_S_orig.mat', 'S');

% set new options.S
options.S = S2;

fname = './results/options.mat';
fprintf('Saving new options struct to file %s...', fname);
save(fname, 'options');
fprintf(' done.\n');

%% Check for TRN, VAL and TST annotation files

% if (exist('./MAT/LFW_TRN_annotation.txt', 'file') && exist('./MAT/LFW_VAL_annotation.txt', 'file') && exist('./MAT/LFW_TST_annotation.txt', 'file'))
%     annotation_struct = read_lfw_annotation_file('./MAT/LFW_TRN_annotation.txt');
%     save('./MAT/TRN.mat', 'annotation_struct');
%     annotation_struct = read_lfw_annotation_file('./MAT/LFW_VAL_annotation.txt');
%     save('./MAT/VAL.mat', 'annotation_struct');
%     annotation_struct = read_lfw_annotation_file('./MAT/LFW_TST_annotation.txt');
%     save('./MAT/TST.mat', 'annotation_struct');
% end;
if (~isempty(TRN_EXAMPLES)  && ~isempty(VAL_EXAMPLES) && ~isempty(TST_EXAMPLES))
    if (exist(['../' TRN_EXAMPLES], 'file') && exist(['../' VAL_EXAMPLES], 'file') && exist(['../' TST_EXAMPLES], 'file'))
        annotation_struct = read_lfw_annotation_file(['../' TRN_EXAMPLES]);
        save('./MAT/TRN.mat', 'annotation_struct');
        annotation_struct = read_lfw_annotation_file(['../' VAL_EXAMPLES]);
        save('./MAT/VAL.mat', 'annotation_struct');
        annotation_struct = read_lfw_annotation_file(['../' TST_EXAMPLES]);
        save('./MAT/TST.mat', 'annotation_struct');
    end;
end;

%% Timestamp

fprintf(1,'\nFinished on %s\n\n', datestr(now));

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function u = SB_ATV(g,mu)
% Split Bregman Anisotropic Total Variation Denoising
%
%   u = arg min_u 1/2||u-g||_2^2 + mu*ATV(u)
%
%   g : noisy image
%   mu: regularisation parameter
%   u : denoised image
%
% Refs:
%  *Goldstein and Osher, The split Bregman method for L1 regularized problems
%   SIAM Journal on Imaging Sciences 2(2) 2009
%  *Micchelli et al, Proximity algorithms for image models: denoising
%   Inverse Problems 27(4) 2011
%
% Benjamin Trémoulhéac
% University College London
% b.tremoulheac@cs.ucl.ac.uk
% April 2012
g = g(:);
n = length(g);
[B Bt BtB] = DiffOper(sqrt(n));
b = zeros(2*n,1);
d = b;
u = g;
err = 1;k = 1;
tol = 1e-3;
lambda = 1;
while err > tol
    fprintf('it. %g ',k);
    up = u;
    [u,~] = cgs(speye(n)+BtB, g-lambda*Bt*(b-d),1e-5,100);
    Bub = B*u+b;
    d = max(abs(Bub)-mu/lambda,0).*sign(Bub);
    b = Bub-d;
    err = norm(up-u)/norm(u);
    fprintf('err=%g \n',err);
    k = k+1;
end
fprintf('Stopped because norm(up-u)/norm(u) <= tol=%.1e\n',tol);
end

function [B Bt BtB] = DiffOper(N)
D = spdiags([-ones(N,1) ones(N,1)], [0 1], N,N+1);
D(:,1) = [];
D(1,1) = 0;
B = [ kron(speye(N),D) ; kron(D,speye(N)) ];
Bt = B';
BtB = Bt*B;
end


text:
function xmat = ptrend(p,nobs)
% PURPOSE: produce an explanatory variables matrix
%          containing a polynomial time-trend
% ----------------------------------------------------
% USAGE: xmat = ptrend(p,nobs);
% where: p = order of the time-trend polynomial
%            p < 0, xmat = iota (nobs x 1) vector of ones
%            p = 1, xmat = time trend
%            p > 1, xmat = higher order polynomial in time
%     nobs = size of the matrix
% ----------------------------------------------------
% RETURNS: xmat = matrix containing polynomial trend model data set
% ----------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jlesage@spatia-econometrics.com

     u = ones(nobs,1) ;
     if p > 0
      timep = zeros(nobs,p) ;
      t = 1:nobs;
      t = (t')/nobs;
       m        = 1 ;
       while (m <= p)
          timep(:,m) = t.^m ;
          m = m + 1 ;
       end;
       xmat = [u timep];
     else,
       xmat = u ;
     end;

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function [favar]=favar_olsfevd(IRFperiods,gamma,favar,n,IRFt,strctident)

% function [fevd_estimates]=bear.olsfevd(irf_estimates,IRFperiods,gamma,n,endo,datapath)
% computes and displays fevd values for the OLS VAR model
% inputs:  - cell 'irf_estimates': lower bound, point estimates, and upper bound for the IRFs
%          - integer 'IRFperiods': number of periods for IRFs
%          - matrix 'gamma': structural disturbance variance-covariance matrix (defined p 48 of technical guide)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - cell 'endo': list of endogenous variables of the model
%          - string 'datapath': user-supplied path to excel data spreadsheet
% outputs: - cell 'fevd_estimates': lower bound, point estimates, and upper bound for the FEVD

% preliminary tasks
npltX=favar.npltX;

if IRFt==1||IRFt==2||IRFt==3
    identified=n; % fully identified
elseif IRFt==4 || IRFt==6 %if the model is identified by sign restrictions or sign restrictions+IV
    identified=size(strctident.signreslabels_shocks,1); % count the labels provided in the sign res sheet (+IV)
elseif IRFt==5
    identified=1; % one IV shock
end

% load IRFs
favar_irf_estimates=favar.IRF.favar_irf_estimates;

% relevant loadings of restricted information variables
L=favar.L(favar.plotX_index,:);

% scale gamma, irf estimates are already scaled
        for ii=1:npltX
            for ll=1:identified
                favar_gamma{ii}(:,ll)=L(ii,ll)*gamma(:,ll);
            end
        end

% create the first cell
temp=cell(npltX,identified+1);

% start by filling the first column of every Tij matrix in the cell
% loop over rows of temp
for jj=1:npltX
   % loop over columns of temp
   for ii=1:identified
   % square each element
   temp{jj,ii}(:,1)=favar_irf_estimates{jj,ii}(:,1).^2;
   end
end
% fill all the other entries of the Tij matrices
% loop over rows of temp
for jj=1:npltX
   % loop over columns of temp
   for ii=1:identified
      % loop over remaining columns
      for kk=2:IRFperiods
      % define the column as the square of the corresponding column in orthogonalised_irf_record
      % additioned to the value of the preceeding columns, which creates the cumulation
      temp{jj,ii}(:,kk)=favar_irf_estimates{jj,ii}(:,kk).^2+temp{jj,ii}(:,kk-1);
      end
   end
end
% multiply each matrix in the cell by the variance of the structural shocks
% loop over rows of temp
for ii=1:npltX
% loop over columns of temp
   for jj=1:identified
   % multiply column jj of the matrix by the variance of the structural shock
   temp{ii,jj}(1,:)=temp{ii,jj}(1,:)*favar_gamma{ii}(jj,jj);
   end
end

% obtain now the values for Ti, the (n+1)th matrix of each row
% loop over rows of temp
for ii=1:npltX
% start the summation over Tij matrices
temp{ii,identified+1}=temp{ii,1};
   % sum over remaining columns
   for jj=2:identified
   temp{ii,identified+1}=temp{ii,identified+1}+temp{ii,jj};
   end
end

% create the output cell fevd_record, scale the shocks with R2 in spirit of BBE (2005)
favar_fevd_estimates=cell(npltX,identified);
R2=bear.favar_R2(favar.X(:,favar.plotX_index),favar.FY,favar.L,favar.plotX_index);
% fill the cell
% loop over rows of fevd_estimates
for ii=1:npltX
    % load the R2 to determine the "true" share of variance explained
    scale=R2(ii);
    shocks=[];
   % loop over columns of fevd_estimates
   for jj=1:identified
   % define the matrix Vfij as the division (pairwise entry) of Tfij by Tfj
   shock=(temp{ii,jj}./temp{ii,identified+1})*scale;
   favar_fevd_estimates{ii,jj}=shock;
   % save shocks to compute residual
   shocks(:,:,jj)=shock;
   end
   % finally add the idiosyncratic component (residual)
   favar_fevd_estimates{ii,jj+1}=1-sum(shocks,3);
end

%save output
favar.FEVD.favar_fevd_estimates=favar_fevd_estimates;

author: european-central-bank

repo: BEAR-toolbox


text:
function [inpoly onboundary] = insidepoly(varargin)
% [inpoly onboundary] = insidepoly(X, Y, PX, PY)
%
% Check if (X,Y) are inside the interior of a 2D polygon delimited by the
% polygon vertices (PX,PY).
%
% INPUTS:
%   - X, Y: arrays of same size, coordinates of N data points
%   - PX, PY: arrays of same size, coordinates of M vertices
%   - Provide optionally EDGE: insidepoly(x, y, Px, Py, edges)
%   EDGES are (m x 2) linear indexes of (Px,Py) that will
%   be considered as vertices of the polygon. This allows user to provide
%   polygons with multiple connexed boundaries (e.g., having holes)
%
% Alternate calling:
%   >> inpoly = insidepoly(XY, P) % or
%   >> inpoly = insidepoly(XY, P, EDGES)
%   - XY: (n x 2) array of n points, ranged by row
%   - P: (m x 2) array of m vertices, ranged by row
%   - EDGES (if provided in third argument) correponds to the
%   first-dimension indexed of P: P(EDGES,:) is the polygone boundary
%
% Advanced control the algorithm:
%   Optionally, user can provide parameters to control the algorithm
%   >> insidepoly(..., 'property1', value1, ...)
% Valid properties (case sensitive) are
%   - 'tol': ['auto'] or positive value: eulidian distance tolerance to
%      detect boundary points. When 'tol' is is set to 'auto', the
%      tolerance value of 1e-9*max(dPx,dPy) is used, where dPx, dPy are
%      respectively the horizontal and vertical size of the polygon.
%   - 'presortflag':  ['auto'], 0, or 1:
%     For large number of vertices and points, INSIDEPOLY uses a sorting
%     strategy to reduce the data to be scanned by each polygonal edge.
%     This reduction comes at the cost of additional processing. The
%     complexity of the sorting step is O((N + M)*log(N+M)). If the numbers
%     of vertices and/or data are small, the sorting is relatively expensive
%     with respect to the complexity of O(M*N) needed for the calculation
%     for the polygonal region determitaion. It might be preferable to
%     disable this step. By default, the presorting is enabled when:
%                      M>32 and M*N>25000.
%     Set the presortflag to 0 or 1 to enable/disable the sorting step.
%
% OUTPUTS:
%   - INPOLY: logical array same size as (X,Y), TRUE if the point is
%   inside or on the boundary of the polygon, FALSE otherwise
%   - ONBOUNDARY: logical array same size as (X,Y) to determine points
%   on the boundary (within numercial tolerance)
%
% This functions has C-MEX engine to speed up the calculation. These Mex
% files must be compiled by lauching insidepoly_install.m
%
% See also: inpolygon, inpoly, insidepoly_dblengine, insidepoly_sglengine
%
% Acknowlegment: The idea of sorting coordinates is first implemented by
% Darren Engwirda, http://www.mathworks.com/matlabcentral/fileexchange/10391
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% History:
%     Original: 06-Jun-2010

% Processing the argument list

% Default options
options = struct('tol', 'auto', ...
                 'presortflag', 'auto');

optionloc = cellfun('isclass',varargin,'char');
if any(optionloc)
    optionloc = find(optionloc,1,'first');
    mainargin = varargin(1:optionloc-1);
    % retreive property/value pairs
    for k=optionloc:2:nargin-1
        options.(varargin{k})= varargin{k+1};
    end
else
    mainargin = varargin;
end

edges = NaN;
if length(mainargin)<=3
    [xy P] = deal(varargin{1:2});
    x = xy(:,1);
    y = xy(:,2);
    if length(mainargin)==3 % edge provided
        edges = varargin{3};
    end
    Px = P(:,1);
    Py = P(:,2);
else
    [x y Px Py] = deal(varargin{1:4});
     if length(mainargin)==5 % edge provided
        edges = varargin{5};
     end
end

edgeprovided = ~isnan(edges);

% Orginal size of the data
sz = size(x);

% Reshape in columns
Px = Px(:);
Py = Py(:);
x = x(:);
y = y(:);

% Cast to double of one of them is (required by Mex)
isxdbl = isa(x,'double');
isydbl = isa(y,'double');
isPxdbl = isa(Px,'double');
isPydbl = isa(Py,'double');
doubleengine = isxdbl || isydbl || isPxdbl|| isPydbl;
if doubleengine
    if ~isxdbl, x = double(x); end
    if ~isydbl, y = double(y); end
    if ~isPxdbl, Px = double(Px); end
    if ~isPydbl, Py = double(Py); end
end

Pxmin = min(Px);
Pxmax = max(Px);
Pymin = min(Py);
Pymax = max(Py);

% Select the tolerance for determining on-boundary points
if isequal(options.tol,'auto')
    ontol = 1e-9*max(Pxmax-Pxmin,Pymax-Pymin);
else
    ontol = options.tol;
end

if ~edgeprovided
    % We don't want duplicate end vertices
    if Px(1)==Px(end) && Py(1)==Py(end)
        Px(end) = [];
        Py(end) = [];
    end
    % Linear wrap around the points
    Px1 = Px;
    Py1 = Py;
    Px2 = Px([2:end 1]);
    Py2 = Py([2:end 1]);
else
    Px1 = Px(edges(:,1),:);
    Py1 = Py(edges(:,1),:);
    Px2 = Px(edges(:,2),:);
    Py2 = Py(edges(:,2),:);
end

% Filter data outside the rectangular box
inpoly = x>=Pxmin-ontol & x<=Pxmax+ontol & ...
         y>=Pymin-ontol & y<=Pymax+ontol;
x = x(inpoly);
y = y(inpoly);

if isequal(options.presortflag,'auto')
    n = size(x,1);
    m = size(Px1,1);
    % We don't do presorting for small size polygon or small size data
    % Empirical law from experimental tests (Bruno)
    presortflag = (m>32) && (n*m > 25000);
else
    presortflag = options.presortflag;
end

if presortflag
    % Sort the array in x, and find the brackets. This is used to find
    % easily which data points have abscissa fall into the abcissa bracket
    % of each edge of the polygon.
    [x ix first last] = presort(Px1, Px2, x);
    y = y(ix); % arrange y in the same order

    % Call mex engine
    if doubleengine
        [in on] = insidepoly_dblengine(x, y, Px1, Py1, Px2, Py2, ontol, ...
                                       first, last);
    else % single arrays
        [in on] = insidepoly_sglengine(x, y, Px1, Py1, Px2, Py2, ontol, ...
                                       first, last);
    end

    % Restore the original order
    in(ix) = in;
    on(ix) = on;
else % No presorting
    % Call mex engine, without presorting
    if doubleengine
        [in on] = insidepoly_dblengine(x, y, Px1, Py1, Px2, Py2, ontol);
    else % single arrays
        [in on] = insidepoly_sglengine(x, y, Px1, Py1, Px2, Py2, ontol);
    end
end

in = in | on;

if nargout>=2
    onboundary = inpoly;
    onboundary(inpoly) = on;
    % Reshape to original size
    onboundary = reshape(onboundary, sz);
end

inpoly(inpoly) = in;
% Reshape to original size
inpoly = reshape(inpoly, sz);

end % insidepoly

%%
function [xsorted ix first last] = presort(Px1, Px2, x)
% (Px1, Px2) abscissa of vertices, x abscissa of data, they are supposed
% to be ranged in column.
% Return:
%   xsorted as sort(x) = x(ix)
%   "first" & "last" indexes such that, for each vertice point
%       min(Px1,Px2) <= xsorted(first:last) <= max(Px1,Px2)

% left and right brackets of the segment
Pmin = min(Px1,Px2);
Pmax = max(Px1,Px2);

nvertices = size(Px1,1);

% We seek to see how x interveaves with Pmin by sorting the ensemble
[trash is] = sort([Pmin; x],1); %#ok
isdata = is>nvertices; % tail index, i.e., belong to data abscissa x
anchor = find(~isdata);
% Get the sorted data alone
ix = is(isdata)-nvertices;
xsorted = x(ix); % sorted x
% Index of the first element in xsorted such that
%   xsorted(first)>=sort(Pmin)
first = anchor-(0:nvertices-1).';
% Rearrange first corresponds to the original order
ip = is(anchor);
first(ip) = first;

% determine how Pmax interleaves with xsorted, i.e.,
% index of the last element in xsorted such that xsorted(first)<=Pmax
% Note: in case of draw in binning edges, HISTC must return the last edge
[trash last] = histc(Pmax, [xsorted; inf]); %#ok

end % presort

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function calpak_test39 ( )

%*****************************************************************************80
%
%% CALPAK_TEST39 tests MONTH_TO_MONTH_NAME_COMMON.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST39\n' );
  fprintf ( 1, '  For the Common calendar,\n' );
  fprintf ( 1, '  MONTH_TO_MONTH_NAME_COMMON names the months:\n' );
  fprintf ( 1, '\n' );

  y = 1;
  months = year_length_months_common ( y );

  for m = 1 : months
    month_name = month_to_month_name_common ( m );
    fprintf ( 1, '  %2d  %s\n', m, month_name );
  end

  return
end


text:
function [p_, KLdiv, exitflag] = optimizeEntropy(p, A, b, Aeq, beq, options)

% number of inequality constraints
K_ = size(A, 1);
% number of equality constraints
K  = size(Aeq, 1);
A_ = A';
%b_ = b';
Aeq_ = Aeq';
beq_ = beq';
x0   = zeros(K_ + K, 1);
InqMat = -eye(K_ + K);
InqMat(K_ + 1:end, :) = [];
InqVec = zeros(K_, 1);

lnp = log(p);

% uses analytic gradient and Hessian during optimization
if nargin < 6
    options = optimset(...
        'Display', 'iter', ...
        'GradObj', 'on', ...
        'Hessian', 'on', ...
        'TolFun', 1e-10, ...
        'TolX', 1e-10, ...
        'TolCon', 1e-10, ...
        'MaxFunEvals', 5000, ...
        'MaxIter', 5000);
end

% minimize negative Lagrange function (i.e. maximize Lagrange function)
if ~K_
    % equality constraints only
    [v, dummy, exitflag] = fminunc(@nestedfunU, x0, options);
    if exitflag > 0
        lnp_ = lnp - 1 - Aeq_ * v;
    else
        lnp_ = lnp;
    end
else
    % inequality (and equality) constraints
    [lv, dummy, exitflag] = fmincon(@nestedfunC, x0, InqMat, InqVec, [], [], [], [], [], options);
    if exitflag > 0
        % inequality Lagrange multipliers
        l = lv(1:K_);
        % equality Lagrange multipliers
        v = lv(K_ + 1:end);
        lnp_ = lnp - 1 - A_ * l - Aeq_ * v;
    else
        lnp_ = lnp;
    end
end

p_ = exp(lnp_);

if nargout > 1
    if exitflag > 0
        KLdiv = p_' * (lnp_ - lnp);
    else
        KLdiv = NaN;
    end
end

    % sub-function -- equality case
    function [mL, g, H] = nestedfunU(v)
        lnx = lnp - 1 - Aeq_ * v;
        % robustificaton
        lnx = max(lnx, -150.0);
        x = exp(lnx);
        % Lagrange dual function
        L = x' * (lnx - lnp + Aeq_ * v) - beq_ * v;
        % take neg values since we want to maximize
        mL = -L;

        % gradient and Hessian
        g = beq - Aeq * x;
        H = Aeq * ((x * ones(1, K)) .* Aeq_);
    end

    % sub-function -- inequality case
    function [mL, g, H] = nestedfunC(lv)
        % inequality Lagrange multiplier
        l = lv(1:K_);
        % equality Lagrange multiplier
        v = lv(K_ + 1:end);
        lnx = lnp - 1 - A_ * l - Aeq_ * v;
        % robustification
        lnx = max(lnx, -150.0);
        x = exp(lnx);
        % Lagrange dual function
        L = x' * (lnx - lnp) + l' * (A * x - b) + v' * (Aeq * x - beq);
        % take neg values since we want to maximize
        mL = -L;

        % gradient and Hessian
        g = [b - A * x; beq - Aeq * x];
        H = [A * ((x * ones(1, K_)) .* A_), A * ((x * ones(1, K)) .* Aeq_);...
            Aeq * ((x * ones(1, K_)) .* A_), Aeq * ((x * ones(1, K)) .* Aeq_)];
    end
end


text:
function varargout = infoto(varargin)
% INFOTO M-file for infoto.fig
%      INFOTO, by itself, creates a new INFOTO or raises the existing
%      singleton*.
%
%      H = INFOTO returns the handle to a new INFOTO or the handle to
%      the existing singleton*.
%
%      INFOTO('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in INFOTO.M with the given input arguments.
%
%      INFOTO('Property','Value',...) creates a new INFOTO or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before infoto_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to infoto_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help infoto

% Last Modified by GUIDE v2.5 09-Jul-2012 15:08:20

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @infoto_OpeningFcn, ...
                   'gui_OutputFcn',  @infoto_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before infoto is made visible.
function infoto_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to infoto (see VARARGIN)

% Choose default command line output for infoto
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes infoto wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = infoto_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes during object creation, after setting all properties.
function axes1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to axes1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes1
axes (hObject)
imshow('infoto.jpg')


text:
%CHEBFUN3V   Class constructor for CHEBFUN3V objects.
%
%   CHEBFUN3V(F, G) constructs a CHEBFUN3V with two components from the
%   function handles F and G. F and G can also be CHEBFUN3 objects or any
%   other object that the CHEBFUN3 constructor accepts. Each component is
%   represented as a CHEBFUN3.
%
%   CHEBFUN3V(F, G, H) constructs a CHEBFUN3V with three components from
%   the function handles F, G, and H.  F, G, and H can also be CHEBFUN3
%   objects or any other objects that the CHEBFUN3 constructor accepts.
%
%   CHEBFUN3V(F, G, [A B C D E K]) constructs a CHEBFUN3V object from F and
%   G on the domain [A B] x [C D] x [E K].
%
%   CHEBFUN3V(F, G, H, [A B C D E K]) constructs a CHEBFUN3V object from F,
%   G, and H on the domain [A B] x [C D] x [E K].
%
% See also CHEBFUN3.

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

classdef chebfun3v

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS PROPERTIES:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    properties ( Access = public )
        components   % Array of CHEBFUN3 objects.
        nComponents  % Number of components
        isTransposed % transposed?
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function F = chebfun3v(varargin)
            % The main CHEBFUN3V constructor!

            % Return an empty CHEBFUN3V:
            if ( (nargin == 0) || isempty(varargin) )
                return
            end

            % This function calls the CHEBFUN3 constructor once for each
            % non-zero component because a CHEBFUN3V is just vector of
            % CHEBFUN3 objects.

            % If argument is a CHEBFUN3V, nothing to do:
            if ( isa(varargin{1}, 'chebfun3v') )
                F = varargin{1};
                return
            end

            % Find the domain:
            domain = [-1 1 -1 1 -1 1];
            for jj = 1:numel(varargin)
               if ( isa(varargin{jj}, 'double') && numel(varargin{jj}) == 6 )
                   domain = varargin{jj};
                   varargin(jj) = [];
               elseif ( isa( varargin{jj}, 'chebfun3') )
                   domain = varargin{jj}.domain;
               end
            end

            % Pick up vectorize flag:
            vectorize = 0;
            for jj = 1:numel(varargin)
                if ( strcmpi( varargin{jj}, 'vectorize' ) )
                    vectorize = 1;
                    varargin(jj) = [];
                end
            end

            % Unwrap input arguments;
            component = 1;
            for jj = 1:numel(varargin)
                if ( iscell(varargin{jj}) )
                    for kk = 1:numel(varargin{jj})
                        fh{component} = varargin{jj}{kk};
                        component = component + 1;
                    end
                else
                    fh{component} = varargin{jj};
                    component = component + 1;
                end
            end
            varargin = fh;

            % Convert all function handles to chebfun3 objects:
            for jj = 1:numel(varargin)
                if ( isa( varargin{jj}, 'function_handle') )
                    if ( ~vectorize )
                        newcheb = chebfun3(varargin{jj}, domain);
                    else
                        newcheb = chebfun3(varargin{jj}, domain, 'vectorize');
                    end
                    fh{jj} = newcheb;
                elseif ( isa(varargin{jj}, 'chebfun3') )
                    fh{jj} = varargin{jj};
                elseif ( isa(varargin{jj}, 'double') )
                    fh{jj} = chebfun3(varargin{jj}, domain);
                end
            end

            % Stop if there are too many components
            if ( numel( fh ) > 3 )
                error('CHEBFUN:CHEBFUN3V:arrayValued', ...
                          'More than three components is not supported.')
            end

            % Stop if there are no components:
            if ( numel(fh) == 0 )
                error('CHEBFUN:CHEBFUN3V:empty', ...
                ['The Chebfun3 constructor needs to be given function ' ...
                     'handles or chebfun3 objects.'])
            end

            % Check the domains of all the chebfun3 objects are the same:
            pass = zeros(numel(fh)-1, 1);
            for jj = 2:numel(fh)
               pass(jj-1) = domainCheck(fh{1}, fh{jj});
            end

            if ( ~all(pass) )
                error('CHEBFUN:CHEBFUN3V:domainCheck', ...
                    'All chebfun3 objects need to have the same domain.');
            end

            % Assign to the Chebfun3v object:
            F.components = fh;
            F.nComponents = numel(fh);
            F.isTransposed = 0;

        end
    end

end


text:
function p = bspeval(d,c,k,u)
%
% Function Name:
%
%   bspeval - Evaluate a univariate B-Spline.
%
% Calling Sequence:
%
%   p = bspeval(d,c,k,u)
%
% Parameters:
%
%   d	: Degree of the B-Spline.
%
%   c	: Control Points, matrix of size (dim,nc).
%
%   k	: Knot sequence, row vector of size nk.
%
%   u	: Parametric evaluation points, row vector of size nu.
%
%   p	: Evaluated points, matrix of size (dim,nu)
%
% Description:
%
%   Evaluate a univariate B-Spline. This function provides an interface to
%   a toolbox 'C' routine.
nu = numel(u);
[mc,nc] = size(c);
                                                %   int bspeval(int d, double *c, int mc, int nc, double *k, int nk, double *u,int nu, double *p){
                                                %   int ierr = 0;
                                                %   int i, s, tmp1, row, col;
                                                %   double tmp2;
                                                %
                                                %   // Construct the control points
                                                %   double **ctrl = vec2mat(c,mc,nc);
                                                %
                                                %   // Contruct the evaluated points
p = zeros(mc,nu);                               %   double **pnt = vec2mat(p,mc,nu);
                                                %
                                                %   // space for the basis functions
N = zeros(d+1,1);                               %   double *N = (double*) mxMalloc((d+1)*sizeof(double));
                                                %
                                                %   // for each parametric point i
for col=1:nu                                    %   for (col = 0; col < nu; col++) {
                                                %     // find the span of u[col]
    s = findspan(nc-1, d, u(col), k);           %     s = findspan(nc-1, d, u[col], k);
    N = basisfun(s,u(col),d,k);                 %     basisfun(s, u[col], d, k, N);
                                                %
    tmp1 = s - d + 1;                           %     tmp1 = s - d;
    for row=1:mc                                %     for (row = 0; row < mc; row++)  {
        tmp2 = 0;                               %       tmp2 = 0.0;
        for i=0:d                               %       for (i = 0; i <= d; i++)
           tmp2 = tmp2 + N(i+1)*c(row,tmp1+i);  % 	tmp2 += N[i] * ctrl[tmp1+i][row];
        end                                     %
        p(row,col) = tmp2;                      %       pnt[col][row] = tmp2;
    end                                         %     }
end                                             %   }
                                                %
                                                %   mxFree(N);
                                                %   freevec2mat(pnt);
                                                %   freevec2mat(ctrl);
                                                %
                                                %   return ierr;
                                                %   }


text:
function [ftH2O] = MPa2ftH2O(MPa)
% Convert pressure from megapascals to feet of water column.
% Chad Greene 2012
ftH2O = MPa*334.553;


text:
function cvx_pop( p, pid )

global cvx___
try
    np = length( cvx___.problems );
catch
    np = 0;
end
if nargin == 0,
    p = np;
    if p == 0,
        if np == 0, return; end
    elseif p > np,
        return
    else
        pstr = cvx___.problems( p );
        pid = pstr.id;
    end
elseif p > np
    pid = 0;
elseif p == 0
    pid = -1;
else
    pstr = cvx___.problems( p );
    if pstr.id ~= pid, pid = 0; end
end
s1 = evalin( 'caller', 'who' );
s2 = evalin( 'caller', 'cellfun(@eval,who,''UniformOutput'',false)' );
s2 = cellfun( @cvx_id, s2 );
do_timer = false;
if pid <= 0,
    do_erase = ~isempty( cvx___.problems );
    p = +do_erase;
    tt = s2 >= 0;
elseif ~pstr.finished
    tt = s2 >= pid;
    do_erase = true;
elseif isnan( pstr.direction )
    tt = s2 == pid;
    do_erase = false;
elseif pstr.checkpoint(4) >= pstr.checkpoint(1),
    tt = s2 > pid;
    temp = sprintf( '%s,', s1{tt} );
    if ~isempty( temp ),
        temp = temp(1:end-1);
        evalin( 'caller', sprintf( '[%s]=cvx_values(%s);', temp, temp ) );
    end
    do_timer = true;
    tt = s2 == pid;
    do_erase = true;
else
    tt = s2 >= pid;
    do_erase = false;
    if ~isempty( pstr.objective ),
        tt = tt & ( s2 ~= cvx_id( pstr.objective ) );
    end
end
s1 = s1(tt);
if ~isempty(s1)
    temp = sprintf( ' %s', s1{:} );
    evalin( 'caller', [ 'clear', temp ] );
end
if do_erase
    cvx_erase( p );
end
if p > 1
    cvx___.problems( p : end ) = [];
else
    cvx___.problems = [];
    cvx_clearpath( 1 );
end
if pid == 0
    cvx_throw( 'Internal CVX data corruption. Please rebuild your model.' );
end
if do_timer
    tfin = tic;
    try
        tval = tfin - pstr.tictime(1);
        cvx___.timers(2) = cvx___.timers(2) + tval;
        cvx___.timers(3) = cvx___.timers(3) + ( tfin - pstr.tictime(2) );
        for q = 1 : p - 1,
            cvx___.problems(q).tictime = cvx___.problems(q).tictime + tval;
        end
    catch
        tval = double(tfin) - double(pstr.tictime(1));
        cvx___.timers(2) = double(cvx___.timers(2)) + double(tval);
        cvx___.timers(3) = double(cvx___.timers(3)) + ( double(tfin) - double(pstr.tictime(2)) );
        for q = 1 : p - 1,
            cvx___.problems(q).tictime = double(cvx___.problems(q).tictime) + tval;
        end
    end
end

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

author: yu-jiang

repo: radpbook

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