text:
function Create_data_train()

load '../models/pdm/pdm_68_aligned_menpo';
load '../models/tri_68.mat';

% This script uses the same format used for patch expert training, and
% expects the data to be there (this can be found in
% https://github.com/TadasBaltrusaitis/CCNF)

% Replace with your location of training data
dataset_loc = 'C:\Users\tbaltrus\Documents\CCNF\patch_experts\data_preparation/prepared_data/';
addpath('../PDM_helpers/');
addpath('./paw_helpers/');

% Collect Menpo, Multi-PIE and 300W data for training the validator
scale = '0.5';
prefix_menpo= 'menpo_train_';
prefix_mpie_300W = 'combined_';

% Find the available positive training data
data_files = dir(sprintf('%s/%s%s*.mat', dataset_loc, prefix_menpo, scale));
data_files_c = dir(sprintf('%s/%s%s*.mat', dataset_loc, prefix_mpie_300W, scale));

centres_all = [];
for i=1:numel(data_files)

    % Load the orientation of the training data
    load([dataset_loc, '/', data_files(i).name], 'centres');
    centres_all = cat(1, centres_all, centres);

end

% Construct mirror indices (which views need to be flipped to create other
% profile training data)
mirror_inds = zeros(size(centres_all,1), 1);
for i=1:numel(data_files)

    % mirrored image has inverse yaw
    mirrored_centre = centres_all(i,:);
    mirrored_centre(2) = -mirrored_centre(2);

    % if mirrored version has same orientation, do not need mirroring
    if(~isequal(mirrored_centre, centres_all(i,:)))

        centres_all = cat(1, centres_all, mirrored_centre);
        mirror_inds = cat(1, mirror_inds, i);

    end

end

% Replace with your location of training data
outputLocation = 'D:\Datasets/detection_validation/prep_data/';

num_more_neg = 10;

% Make sure same data generated all the time
rng(0);

% Negative samples from teh INRIAPerson dataset
neg_image_loc = 'D:\Datasets\INRIAPerson\INRIAPerson\Train\neg/';

neg_images = cat(1,dir([neg_image_loc, '/*.jpg']),dir([neg_image_loc, '/*.png']));

max_img_used = 8000;

% do it separately for centers due to memory limitations
for r=1:size(centres_all,1)

    a_mod = 0.4;

    mirror = false;

    if(mirror_inds(r) ~= 0 )
        mirror = true;
        label_mirror_inds = [1,17;2,16;3,15;4,14;5,13;6,12;7,11;8,10;18,27;19,26;20,25;21,24;22,23;...
                  32,36;33,35;37,46;38,45;39,44;40,43;41,48;42,47;49,55;50,54;51,53;60,56;59,57;...
                  61,65;62,64;68,66];
        load([dataset_loc, '/', data_files_c(mirror_inds(r)).name]);
        all_images_t = all_images;
        landmark_locations_t = landmark_locations;
        visiIndex_t = visiIndex;

        load([dataset_loc, '/', data_files(mirror_inds(r)).name]);

        % Combining Menpo + MPIE + 300W
        all_images = cat(1, all_images, all_images_t);
        landmark_locations = cat(1, landmark_locations, landmark_locations_t);
        % Taking a subset of visibilities from all the datasets
        visiIndex = visiIndex_t & visiIndex;
    else
        load([dataset_loc, '/', data_files_c(r).name]);
        all_images_t = all_images;
        landmark_locations_t = landmark_locations;
        visiIndex_t = visiIndex;

        load([dataset_loc, '/', data_files(r).name]);
        all_images = cat(1, all_images, all_images_t);
        landmark_locations = cat(1, landmark_locations, landmark_locations_t);
        visiIndex = visiIndex_t & visiIndex;
    end

    visiCurrent = logical(visiIndex);

    if(mirror)
        centres = [centres(1), -centres(2), -centres(3)];
        tmp1 = visiCurrent(label_mirror_inds(:,1));
        tmp2 = visiCurrent(label_mirror_inds(:,2));
        visiCurrent(label_mirror_inds(:,2)) = tmp1;
        visiCurrent(label_mirror_inds(:,1)) = tmp2;
    end

    visibleVerts = 1:numel(visiCurrent);
    visibleVerts = visibleVerts(visiCurrent)-1;

    % Correct the triangulation to take into account the vertex
    % visibilities
    triangulation = [];

    shape = a_mod * Euler2Rot(centres * pi/180) * reshape(M, numel(M)/3, 3)';
    shape = shape';

    for i=1:size(T,1)
        visib = 0;
        for j=1:numel(visibleVerts)
            if(T(i,1)==visibleVerts(j))
                visib = visib+1;
            end
            if(T(i,2)==visibleVerts(j))
                visib = visib+1;
            end
            if(T(i,3)==visibleVerts(j))
                visib = visib+1;
            end
        end

        % Only if all three of the vertices are visible
        if(visib == 3)

            % Also want to remove triangles facing the wrong way (self occluded)
            v1 = [shape(T(i,1)+1,1), shape(T(i,1)+1,2), shape(T(i,1)+1,3)];
            v2 = [shape(T(i,2)+1,1), shape(T(i,2)+1,2), shape(T(i,2)+1,3)];
            v3 = [shape(T(i,3)+1,1), shape(T(i,3)+1,2), shape(T(i,3)+1,3)];
            normal = cross((v2-v1), v3 - v2);
            normal = normal / norm(normal);
            direction = normal * [0,0,1]';

            % And only if the triangle is facing the camera
            if(direction > 0)
                triangulation = cat(1, triangulation, T(i,:));
            end
        end
    end

    % Initialise the warp
    [ alphas, betas, triX, mask, minX, minY, nPix ] = InitialisePieceWiseAffine(triangulation, shape);

    mask = logical(mask);

    imgs_to_use = randperm(size(landmark_locations, 1));

    if(size(landmark_locations, 1) > max_img_used)
        imgs_to_use = imgs_to_use(1:max_img_used);
    end

    % Extracting relevant filenames
    examples = zeros(numel(imgs_to_use) * (num_more_neg+1), nPix);
    errors = zeros(numel(imgs_to_use) * (num_more_neg+1), 1);

    unused_pos = 0;

    curr_filled = 0;

    for j=imgs_to_use

        labels = squeeze(landmark_locations(j,:,:));

        img = squeeze(all_images(j,:,:));

        if(mirror)
            img = fliplr(img);
            imgSize = size(img);
            flippedLbls = labels;
            flippedLbls(:,1) = imgSize(1) - flippedLbls(:,1) + 1;
            tmp1 = flippedLbls(label_mirror_inds(:,1),:);
            tmp2 = flippedLbls(label_mirror_inds(:,2),:);
            flippedLbls(label_mirror_inds(:,2),:) = tmp1;
            flippedLbls(label_mirror_inds(:,1),:) = tmp2;
            labels = flippedLbls;
        end

        % If for some reason some of the labels are not visible in the
        % current sample skip this label
        non_existent_labels = labels(:,1)==0 | labels(:,2)==0;
        non_existent_inds = find(non_existent_labels)-1;
        if(numel(intersect(triangulation(:), non_existent_inds)) > 0)
           unused_pos = unused_pos + 1;
           continue;
        end

        % Centering the pixel so that 0,0 is center of the top left pixel
        labels = labels - 1;

        curr_filled = curr_filled + 1;
        [features] = ExtractFaceFeatures(img, labels, triangulation, triX, mask, alphas, betas, nPix, minX, minY);
%         sample_img = zeros(size(mask));sample_img(mask) = features;imagesc(sample_img)
        examples(curr_filled,:) = features;
        errors(curr_filled,:) = 0;

        % Extract the correct PDM parameters for the model (we will perturb
        % them for some negative examples)
        [ a_orig, R_orig, trans_orig, ~, params_orig] = fit_PDM_ortho_proj_to_2D(M, E, V, labels);
        eul_orig = Rot2Euler(R_orig);

        % a slightly perturbed example, too tight
        % from 0.3 to 0.9
        a_mod = a_orig * (0.6 + (randi(7) - 4)*0.1);
        p_global = [a_mod; eul_orig'; trans_orig];

        labels_mod = GetShapeOrtho(M, V, params_orig, p_global);
        labels_mod = labels_mod(:,1:2);

        [features] = ExtractFaceFeatures(img, labels_mod, triangulation, triX, mask, alphas, betas, nPix, minX, minY);
%         sample_img = zeros(size(mask));sample_img(mask) = features;imagesc(sample_img)

        curr_filled = curr_filled + 1;
        examples(curr_filled,:) = features;

        % Compute the badness of fit
        error = norm(labels_mod(:) - labels(:)) / (max(labels(:,2))-min(labels(:,2)));
        errors(curr_filled,:) = error;

        % a slightly perturbed example, too broad
        % from 1.2 to 0.6
        a_mod = a_orig * (1.4 + (randi(5) - 3)*0.1);
        p_global = [a_mod; eul_orig'; trans_orig];

        labels_mod = GetShapeOrtho(M, V, params_orig, p_global);
        labels_mod = labels_mod(:,1:2);

        [features] = ExtractFaceFeatures(img, labels_mod, triangulation, triX, mask, alphas, betas, nPix, minX, minY);
%         sample_img = zeros(size(mask));sample_img(mask) = features;imagesc(sample_img)

        curr_filled = curr_filled + 1;
        examples(curr_filled,:) = features;

        error = norm(labels_mod(:) - labels(:)) / (max(labels(:,2))-min(labels(:,2)));
        errors(curr_filled,:) = error;

        % A somewhat offset example

        trans_mod = trans_orig + randn(2,1) * 20;
        p_global = [a_orig; eul_orig'; trans_mod];

        labels_mod = GetShapeOrtho(M, V, params_orig, p_global);
        labels_mod = labels_mod(:,1:2);

        [features] = ExtractFaceFeatures(img, labels_mod, triangulation, triX, mask, alphas, betas, nPix, minX, minY);

        curr_filled = curr_filled + 1;
        examples(curr_filled,:) = features;

        error = norm(labels_mod(:) - labels(:)) / (max(labels(:,2))-min(labels(:,2)));
        errors(curr_filled,:) = error;

        % A rotated sample
        eul_mod = eul_orig + randn(1,3)*0.3;
        p_global = [a_orig; eul_mod'; trans_orig];

        labels_mod = GetShapeOrtho(M, V, params_orig, p_global);
        labels_mod = labels_mod(:,1:2);

        [features] = ExtractFaceFeatures(img, labels_mod, triangulation, triX, mask, alphas, betas, nPix, minX, minY);

        curr_filled = curr_filled + 1;
        examples(curr_filled,:) = features;

        error = norm(labels_mod(:) - labels(:)) / (max(labels(:,2))-min(labels(:,2)));
        errors(curr_filled,:) = error;

        % A sample with modified shape parameters
        p_global = [a_orig; eul_orig'; trans_orig];
        params_mod = params_orig + randn(size(params_orig)).*sqrt(E);
        labels_mod = GetShapeOrtho(M, V, params_mod, p_global);
        labels_mod = labels_mod(:,1:2);

        [features] = ExtractFaceFeatures(img, labels_mod, triangulation, triX, mask, alphas, betas, nPix, minX, minY);

        curr_filled = curr_filled + 1;
        examples(curr_filled,:) = features;

        error = norm(labels_mod(:) - labels(:)) / (max(labels(:,2))-min(labels(:,2)));
        errors(curr_filled,:) = error;

        % pick a random image from negative inriaperson dataset, use original location if
        % first, otherwhise resize it to fit
        for n=6:num_more_neg
            n_img = randi(numel(neg_images));

            neg_image = imread([neg_image_loc, neg_images(n_img).name]);

            if(size(neg_image,3) == 3)
                neg_image = rgb2gray(neg_image);
            end

            [h_neg, w_neg] = size(neg_image);

            % if the current labels fit just use them, if not, then resize
            % to fit
            max_x = max(labels(:,1));
            max_y = max(labels(:,2));

            if(max_x > w_neg || max_y > h_neg)
               neg_image = imresize(neg_image, [max_y, max_x]);
            end

            [features] = ExtractFaceFeatures(neg_image, labels, triangulation, triX, mask, alphas, betas, nPix, minX, minY);

            curr_filled = curr_filled + 1;
            examples(curr_filled,:) = features;

            % Set high error to 3
            errors(curr_filled,:) = 3;
        end

        if(mod(curr_filled, 10) == 0)
            fprintf('%d/%d done\n', curr_filled/(num_more_neg+1), numel(imgs_to_use));
        end
        % add the pos example to the background

    end

    examples = examples(1:curr_filled,:);
    errors = errors(1:curr_filled);

    filename = sprintf('%s/face_validator_train_%d.mat', outputLocation, r);
    save(filename, 'examples', 'errors', 'alphas', 'betas', 'triangulation', 'minX', 'minY', 'nPix', 'shape', 'triX', 'mask', 'centres');

end

end

function [features] = ExtractFaceFeatures(img, labels, triangulation, triX, mask, alphas, betas, nPix, minX, minY)

    % Make sure labels are within range
    [hRes, wRes] = size(img);
    labels(labels(:,1) < 0,1) = 0;
    labels(labels(:,2) < 0,2) = 0;

    labels(labels(:,1) > wRes-1,1) = wRes-1;
    labels(labels(:,2) > hRes-1,2) = hRes-1;

    crop_img = Crop(img, labels, triangulation, triX, mask, alphas, betas, nPix, minX, minY);
    crop_img(isnan(crop_img)) = 0;

    % vectorised version
    features = reshape(crop_img(logical(mask)), 1, nPix);

    % normalisations
    features = (features - mean(features));
    norms = std(features);
    if(norms==0)
        norms = 1;
    end
    features = features / norms;

end


text:
function fem_basis_t6_display ( prefix )

%*****************************************************************************80
%
%% FEM_BASIS_T6_DISPLAY displays a finite element T6 basis.
%
%  Discussion:
%
%    This program reads a data file defining a set of nodes, and a
%    data file defining the triangulation of those nodes using 6 node triangles.
%
%    It is assumed that the 6 node triangles have straight sides, that is,
%    that each midside node lies on the line segment connecting the
%    corresponding vertices.
%
%    The program then asks the user interactively to select one of the
%    nodes.  It computes the basis function associated with that node
%    and displays it over the entire mesh.  Of course, the basis function
%    will only be nonzero over a small number of the elements, but it
%    is instructive to see the entire mesh.
%
%    The display is initially "flat", but by using the manipulator
%    on the graphics menu, the user can easily get some dramatic images
%    of the basis function.
%
%  Usage:
%
%    fem_basis_t6_display ( 'prefix' )
%
%    where 'prefix' is the common prefix for the FEM files:
%
%    * 'prefix'_nodes.txt,    the node coordinates.
%    * 'prefix'_elements.txt, the nodes that make up each element;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 March 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM_BASIS_T6_DISPLAY:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Display basis functions associated with \n' );
  fprintf ( 1, '  a finite element grid of quadratic triangles ("T6").\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This program reads two files:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * node_file,     the node file,\n' );
  fprintf ( 1, '  * element_file,  the element file,\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The user specifies a basis function by node index,\n' );
  fprintf ( 1, '  and the program displays a surface plot of that basis function.\n' );
  fprintf ( 1, '  (Use the 3D ROTATE option to see the full picture!\n' );
%
%  Get the prefix if missing.
%
  if ( nargin < 1 )
    prefix = input ( '  Enter the common file prefix:  ' );
  end
%
%  Construct the file names.
%
  node_filename = strcat ( prefix, '_nodes.txt' );
  element_filename = strcat ( prefix, '_elements.txt' );
%
%  Read the nodes.
%
  [ dim_num, node_num ] = r8mat_header_read ( node_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension = %d\n', dim_num );
  fprintf ( 1, '  Number of nodes = %d\n', node_num );

  if ( dim_num ~= 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FEM_BASIS_T6_DISPLAY - Fatal error!\n' );
    fprintf ( 1, '  The spatial dimension is not 2.\n' );
    error ( 'FEM_BASIS_T6_DISPLAY - Fatal error!' );
  end

  node_xy = r8mat_data_read ( node_filename, dim_num, node_num );
%
%  Read the elements.
%
  [ element_order, element_num ] = i4mat_header_read ( element_filename );

  if ( element_order ~= 6 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FEM_BASIS_T6_DISPLAY\n' );
    fprintf ( 1, '  This program requires that the elements have order 6.\n' );
    fprintf ( 1, '  Your elements have order %d\n', element_order );
    error ( 'FEM_BASIS_T6_DISPLAY - Fatal error!' );
  end

  element_node = i4mat_data_read ( element_filename, element_order, ...
    element_num );
%
%  Set up the graph.
%
  x_min = min ( node_xy(1,:) );
  x_max = max ( node_xy(1,:) );
  y_min = min ( node_xy(2,:) );
  y_max = max ( node_xy(2,:) );

  node_min = min ( min ( element_node ) );
  node_max = max ( max ( element_node ) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Every basis function is associated with a node.\n' );
  fprintf ( 1, '  To chooose a basis function, you specify a node.\n' );
  fprintf ( 1, '  Nodes range in value from %d to %d.\n', node_min, node_max );

  while ( 1 )

    fprintf ( 1, '\n' );
    prompt = sprintf ( 'Enter a node between %d and %d: ', node_min, node_max );
    node_index = input ( prompt );

    if ( node_index < node_min | node_max < node_index )
      break
    end
%
%  Clear the graphics page.
%
    clf

    fprintf ( 1, '\n' );

    z_min = 0.0;
    z_max = 0.0;

    for element = 1 : element_num

      pt(1,1:6) = node_xy(1,element_node(1:6,element));
      pt(2,1:6) = node_xy(2,element_node(1:6,element));

      local = 0;
      for j = 1 : element_order
        if ( element_node(j,element) == node_index )
          local = j;
          fprintf ( 1, '  Node %d occurs as local node %d in element %d.\n', ...
            node_index, j, element );
        end
      end
%
%  If the node is not part of this element, we can draw the zero patch easily.
%
      if ( local == 0 )

        x = zeros(2,2);
        y = zeros(2,2);
        z = zeros(2,2);
        x(1,1) = pt(1,1);
        y(1,1) = pt(2,1);
        x(2,1) = pt(1,1);
        y(2,1) = pt(2,1);
        x(1,2) = pt(1,2);
        y(1,2) = pt(2,2);
        x(2,2) = pt(1,3);
        y(2,2) = pt(2,3);
%
%  If the node is part of this element, then we want to use an 11x11 grid of points
%  to sample the (reference) element and generate data for a smooth surface.
%
      else

        x = zeros(11,11);
        y = zeros(11,11);
        z = zeros(11,11);

        for i = 0 : 10
%
%  To use SURFACE, we need to define a logical rectangular grid of points.
%  But our region is a triangle, so the easiest thing to do is to identify
%  two corners of the grid.
%  We try to do this in such a way that the same thing happens in each
%  element, and the resulting grid looks less haphazard.
%
          if ( local == 1 | local == 5 )
            xlo = ( ( 10 - i ) * pt(1,1) + i * pt(1,2) ) / 10;
            ylo = ( ( 10 - i ) * pt(2,1) + i * pt(2,2) ) / 10;
            xhi = ( ( 10 - i ) * pt(1,1) + i * pt(1,3) ) / 10;
            yhi = ( ( 10 - i ) * pt(2,1) + i * pt(2,3) ) / 10;
          elseif ( local == 2 | local == 6 )
            xlo = ( ( 10 - i ) * pt(1,2) + i * pt(1,1) ) / 10;
            ylo = ( ( 10 - i ) * pt(2,2) + i * pt(2,1) ) / 10;
            xhi = ( ( 10 - i ) * pt(1,2) + i * pt(1,3) ) / 10;
            yhi = ( ( 10 - i ) * pt(2,2) + i * pt(2,3) ) / 10;
          elseif ( local == 3 | local == 4 )
            xlo = ( ( 10 - i ) * pt(1,3) + i * pt(1,2) ) / 10;
            ylo = ( ( 10 - i ) * pt(2,3) + i * pt(2,2) ) / 10;
            xhi = ( ( 10 - i ) * pt(1,3) + i * pt(1,1) ) / 10;
            yhi = ( ( 10 - i ) * pt(2,3) + i * pt(2,1) ) / 10;
          end

          for j = 0 : 10

            pp(1) = ( ( 10 - j ) * xlo + j * xhi ) / 10;
            pp(2) = ( ( 10 - j ) * ylo + j * yhi ) / 10;

            zz = basis_11_t6 ( pt, local, pp );

            z_min = min ( z_min, zz );
            z_max = max ( z_max, zz );

            x(i+1,j+1) = pp(1);
            y(i+1,j+1) = pp(2);
            z(i+1,j+1) = zz;

          end
        end

      end

      caxis ( 'auto' )

      surface ( x, y, z, 'FaceColor', 'interp' )

    end
%
%  Using SHADING FACETED draws nice little boundary lines between each patch.
%  However, the way we draw the basis function means that that surface is simply
%  covered in black lines, losing all the color information.
%
%  SHADING INTERP loses the little boundary lines, but allows the color information
%  to be seen.
%
%  You have your clumsy choice of ONE or THE OTHER.
%
    shading faceted;
%   shading interp;

    axis ( [ x_min, x_max, y_min, y_max, z_min, z_max ] );
    axis equal

    xlabel ( '--X axis--' )
    ylabel ( '--Y axis--' )
    zlabel ( '--Z axis--' )

    title_string = sprintf ( 'T6 basis function for node %d', node_index );
    title ( title_string );

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM_BASIS_T6_DISPLAY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ bi, dbidx, dbidy ] = basis_11_t6 ( t, i, p )

%*****************************************************************************80
%
%% BASIS_11_T6: one basis at one point for the T6 element.
%
%  Discussion:
%
%    The routine is given the coordinates of the nodes of a triangle.
%
%           3
%          / \
%         6   5
%        /     \
%       1---4---2
%
%    It evaluates the quadratic basis function B(I)(X,Y) associated with
%    node I, which has the property that it is a quadratic function
%    which is 1 at node I and zero at the other five nodes.
%
%    This routine assumes that the sides of the triangle are straight,
%    so that the midside nodes fall on the line between two vertices.
%
%    This routine relies on the fact that each basis function can be
%    written as the product of two linear factors, which are easily
%    computed and normalized.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,6), the coordinates of the nodes.
%
%    Input, integer I, the index of the desired basis function.
%    I should be between 1 and 6.
%
%    Input, real P(2), the coordinates of a point at which the basis
%    function is to be evaluated.
%
%    Output, real BI, DBIDX, DBIDY, the values of the basis function
%    and its X and Y derivatives.
%
  if ( i < 1 | 6 < i )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BASIS_11_T6 - Fatal error!\n' );
    fprintf ( 1, '  Basis index I is not between 1 and 6.\n' );
    fprintf ( 1, '  I = ', i );
    error ( 'BASIS_11_T6 - Fatal error!' );
  end
%
%  Determine the pairs of nodes.
%
  if ( i <= 3 )
    j1 = i4_wrap ( i + 1, 1, 3 );
    j2 = i4_wrap ( i + 2, 1, 3 );
    k1 = i + 3;
    k2 = i4_wrap ( i + 5, 4, 6 );
  else
    j1 = i - 3;
    j2 = i4_wrap ( i - 3 + 2, 1, 3 );
    k1 = i4_wrap ( i - 3 + 1, 1, 3 );
    k2 = i4_wrap ( i - 3 + 2, 1, 3 );
  end
%
%  Evaluate the two linear factors GF and HF,
%  and their normalizers GN and HN.
%
  gf = ( p(1)    - t(1,j1) ) * ( t(2,j2) - t(2,j1) ) ...
     - ( t(1,j2) - t(1,j1) ) * ( p(2)    - t(2,j1) );

  gn = ( t(1,i)  - t(1,j1) ) * ( t(2,j2) - t(2,j1) ) ...
     - ( t(1,j2) - t(1,j1) ) * ( t(2,i)  - t(2,j1) );

  hf = ( p(1)    - t(1,k1) ) * ( t(2,k2) - t(2,k1) ) ...
     - ( t(1,k2) - t(1,k1) ) * ( p(2)    - t(2,k1) );

  hn = ( t(1,i)  - t(1,k1) ) * ( t(2,k2) - t(2,k1) ) ...
     - ( t(1,k2) - t(1,k1) ) * ( t(2,i)  - t(2,k1) );
%
%  Construct the basis function and its derivatives.
%
  bi =        ( gf                  / gn ) * (   hf                  / hn );

  dbidx =   ( ( t(2,j2) - t(2,j1) ) / gn ) * (   hf                  / hn ) ...
          + (   gf                  / gn ) * ( ( t(2,k2) - t(2,k1) ) / hn );

  dbidy = - ( ( t(1,j2) - t(1,j1) ) / gn ) * (   hf                  / hn ) ...
          - (   gf                  / gn ) * ( ( t(1,k2) - t(1,k1) ) / hn );

  return
end
function column_num = file_column_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_COLUMN_COUNT counts the columns in the first line of a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%    Most lines of the file are presumed to consist of COLUMN_NUM words,
%    separated by spaces.  There may also be some blank lines, and some
%    comment lines, which have a "#" in column 1.
%
%    The routine tries to find the first non-comment non-blank line and
%    counts the number of words in that line.
%
%    If all lines are blanks or comments, it goes back and tries to analyze
%    a comment line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the file.
%
%    Output, integer COLUMN_NUM, the number of columns in the file.
%
  FALSE = 0;
  TRUE = 1;
%
%  Open the file.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_COLUMN_COUNT - Error!' );
    column_num = -1;
    return;
  end
%
%  Read one line, but skip blank lines and comment lines.
%  Use FGETL so we drop the newline character!
%
  got_one = FALSE;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( s_len_trim ( line ) == 0 )

    elseif ( line(1) == '#' )

    else
      got_one = TRUE;
      break;
    end

  end

  fclose ( input_unit );

  if ( got_one == FALSE )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_COLUMN_COUNT - Warning!\n' );
    fprintf ( 1, '  The file does not seem to contain any data.\n' );
    column_num = -1;
    return;
  end

  column_num = s_word_count ( line );

  return
end
function row_num = file_row_count ( input_file_name )

%*****************************************************************************80
%
%% FILE_ROW_COUNT counts the number of row records in a file.
%
%  Discussion:
%
%    Each input line is a "RECORD".
%
%    The records are divided into three groups:
%
%    * BLANK LINES (nothing but blanks)
%    * COMMENT LINES (begin with a '#')
%    * DATA RECORDS (anything else)
%
%    The value returned by the function is the number of data records.
%
%    By the way, if the MATLAB routine FGETS is used, instead of
%    FGETL, then the variable LINE will include line termination
%    characters, which means that a blank line would not actually
%    have zero characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Output, integer ROW_NUM, the number of rows found.
%
  input_unit = fopen ( input_file_name );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_ROW_COUNT - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'FILE_ROW_COUNT - Error!' );
    row_num = -1;
    return;
  end

  blank_num = 0;
  comment_num = 0;
  row_num = 0;

  record_num = 0;

  while ( 1 )

    line = fgetl ( input_unit );

    if ( line == -1 )
      break;
    end

    record_num = record_num + 1;
    record_length = s_len_trim ( line );

    if ( record_length <= 0 )
      blank_num = blank_num + 1;
    elseif ( line(1) == '#' )
      comment_num = comment_num + 1;
    else
      row_num = row_num + 1;
    end

  end

  fclose ( input_unit );

  return
end
function value = i4_modp ( i, j )

%*****************************************************************************80
%
%% I4_MODP returns the nonnegative remainder of I4 division.
%
%  Discussion:
%
%    If
%      NREM = I4_MODP ( I, J )
%      NMULT = ( I - NREM ) / J
%    then
%      I = J * NMULT + NREM
%    where NREM is always nonnegative.
%
%    The MOD function computes a result with the same sign as the
%    quantity being divided.  Thus, suppose you had an angle A,
%    and you wanted to ensure that it was between 0 and 360.
%    Then mod(A,360) would do, if A was positive, but if A
%    was negative, your result would be between -360 and 0.
%
%    On the other hand, I4_MODP(A,360) is between 0 and 360, always.
%
%  Example:
%
%        I     J     MOD  I4_MODP    Factorization
%
%      107    50       7       7    107 =  2 *  50 + 7
%      107   -50       7       7    107 = -2 * -50 + 7
%     -107    50      -7      43   -107 = -3 *  50 + 43
%     -107   -50      -7      43   -107 =  3 * -50 + 43
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 1999
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number to be divided.
%
%    Input, integer J, the number that divides I.
%
%    Output, integer VALUE, the nonnegative remainder when I is
%    divided by J.
%
  if ( j == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_MODP - Fatal error!\n' );
    fprintf ( 1, '  Illegal divisor J = %d\n', j );
    error ( 'I4_MODP - Fatal error!' );
  end

  value = mod ( i, j );

  if ( value < 0 )
    value = value + abs ( j );
  end

  return
end
function value = i4_wrap ( ival, ilo, ihi )

%*****************************************************************************80
%
%% I4_WRAP forces an integer to lie between given limits by wrapping.
%
%  Example:
%
%    ILO = 4, IHI = 8
%
%    I   Value
%
%    -2     8
%    -1     4
%     0     5
%     1     6
%     2     7
%     3     8
%     4     4
%     5     5
%     6     6
%     7     7
%     8     8
%     9     4
%    10     5
%    11     6
%    12     7
%    13     8
%    14     4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer IVAL, an integer value.
%
%    Input, integer ILO, IHI, the desired bounds for the integer value.
%
%    Output, integer I4_WRAP, a "wrapped" version of IVAL.
%
  jlo = min ( ilo, ihi );
  jhi = max ( ilo, ihi );

  wide = jhi - jlo + 1;

  if ( wide == 1 )
    value = jlo;
  else
    value = jlo + i4_modp ( ival - jlo, wide );
  end

  return
end
function table = i4mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% I4MAT_DATA_READ reads data from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns in the data.
%
%    Output, integer TABLE(M,N), the point coordinates.
%
  table = [];
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %d' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the input file.\n' );
    error ( 'I4MAT_DATA_READ - Error!' );
    return;
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'I4MAT_DATA_READ - Error!\n' );
      fprintf ( 1, '  End of input while reading data.\n' );
      error ( 'I4MAT_DATA_READ - Error!' );
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = i4mat_header_read ( input_filename )

%*****************************************************************************80
%
%% I4MAT_HEADER_READ reads the header from an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function table = r8mat_data_read ( input_filename, m, n )

%*****************************************************************************80
%
%% R8MAT_DATA_READ reads data from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Input, integer M, N, the number of rows and columns of data.
%
%    Output, real TABLE(M,N), the point coordinates.
%
  table = [];
%
%  Build up the format string for reading M real numbers.
%
  string = ' ';

  for i = 0 : m
    string = strcat ( string, ' %f' );
  end

  input_unit = fopen ( input_filename );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_DATA_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file.\n' );
    error ( 'R8MAT_DATA_READ - Error!' );
    return;
  end

  i = 0;

  while ( i < n )

    line = fgets ( input_unit );

    if ( line == -1 )
      break;
    end

    if ( line(1) == '#' )

    elseif ( s_len_trim ( line ) == 0 )

    else

      [ x, count ] = sscanf ( line, string );

      if ( count == m )
        i = i + 1;
        table(1:m,i) = x(1:m);
      end

    end

  end

  fclose ( input_unit );

  return
end
function [ m, n ] = r8mat_header_read ( input_filename )

%*****************************************************************************80
%
%% R8MAT_HEADER_READ reads the header from an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILENAME, the name of the input file.
%
%    Output, integer M, the spatial dimension.
%
%    Output, integer N, the number of points.
%
  m = file_column_count ( input_filename );

  if ( m <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data columns in\n' );
    fprintf ( 1, '  the file %s.\n', input_filename );
  end

  n = file_row_count ( input_filename );

  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  There was some kind of I/O problem while trying\n' );
    fprintf ( 1, '  to count the number of data rows in\n' );
    fprintf ( 1, '  the file %s\n', input_filename );
  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LENGTH, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  FALSE = 0;
  TRUE = 1;

  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = TRUE;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = TRUE;
    elseif ( blank == TRUE )
      word_num = word_num + 1;
      blank = FALSE;
    end

  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [d,dt]=formatts(d)
%
% Usage: [d,dt]=formatts(d)
%
% Helper function for CWT,XWT,WTC
%
% Brings a timeseries into a shape so that it has two columns: [time, value].
%
%
% (C) Aslak Grinsted 2002-2014
% http://www.glaciology.net/wavelet-coherence

% -------------------------------------------------------------------------
%The MIT License (MIT)
%
%Copyright (c) 2014 Aslak Grinsted
%
%Permission is hereby granted, free of charge, to any person obtaining a copy
%of this software and associated documentation files (the "Software"), to deal
%in the Software without restriction, including without limitation the rights
%to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
%copies of the Software, and to permit persons to whom the Software is
%furnished to do so, subject to the following conditions:
%
%The above copyright notice and this permission notice shall be included in
%all copies or substantial portions of the Software.
%
%THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
%IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
%FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
%AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
%LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
%OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
%THE SOFTWARE.
%---------------------------------------------------------------------------

if (ndims(d)>2)
    error('Input time series should be 2 dimensional.');
end
if (numel(d)==length(d))
    d=[(1:length(d))' d(:)];
end
if size(d,1)<size(d,2)
    d=d';
end
if (size(d,2)~=2)
    error('Time series must have 2 columns.')
end

if (d(2,1)-d(1,1)<0)
    d=flipud(d);
end

dt=diff(d(:,1));
if any(abs(dt-dt(1))>1e-1*dt(1))
    error('Time step must be constant.');
end
if (dt==0)
    error('Time step must be greater than zero.')
end

dt=dt(1);

author: grinsted

repo: wavelet-coherence


text:
function [solution2] = computeMin2Norm_HH(model)
% This function computes the min2Norm for the whole-body metabolic model.
% As the underlying quadratic programming problem may not be proven to
% yield an optiomal solution due to numercial difficulties, we compute the
% min2Norm twice.
% In the first QP, we take the input model as is, and compute the feasible
% but not necessarily optimal flux distribution. In a second step, we
% rescale the bounds on the model's reaction to ensure numerical stability.
% The computed solution from the second QP is rescaled accordingly
% afterwards.
% Note that this function assumes that the infinity constraints are represented by -1,000,000 or 1,000,000.
%
% INPUT
% model      model structure (whole-body metabolic model)
%
% OUTPUT
% solution2  optimal min2Norm solution
%
% Ines Thiele 01/2018

model.c = zeros(length(model.c),1);

% First QP on vanilla model
[solution,LPProblem]=solveCobraLPCPLEX(model,0,0,0,[],1e-6,'tomlab_cplex');

% Prepare 2nd QP
model2=model;
% all high flux values in the solution vectore of the first QP retain a high bound
model2.lb(find(solution.full<-1e4))=-500000;
model2.ub(find(solution.full>1e4))=500000;

% reduced the "infinity" bounds on all other reactions.
% note this step does not affect any non-infinity bounds set on the
% whole-body metabolic model
model2.lb(find(model2.lb==-1000000))=-10000; % reduce the effective unbound constraints to lower number, representing inf
model2.ub(find(model2.ub==1000000))=10000;% reduce the effective unbound constraints to lower number, representing inf

% we then rescale all bounds on the model reactions by a factor of 1/1000,
% which proven to result in an optimal QP solution
model2.lb=model2.lb/1000;
model2.ub=model2.ub/1000;

if solution.origStat == 6
    % we solve the 2nd QP
    [solution2,LPProblem]=solveCobraLPCPLEX(model2,1,0,0,[],1e-6,'tomlab_cplex');

    % we resale the computed solution by the factor of 1000
    solution2.full = solution2.full*1000; % rescale solution by the factor of 1000
    solution2.obj = solution2.obj*1000; % rescale solution by the factor of 1000
else
    solution2 = solution;
end

author: opencobra

repo: cobratoolbox


text:
% Create a canlab_dataset object from a Matlab Table object
%
% The canlab_dataset object provides a standard data format, and allows the use of more methods,
% such as writing the data object to text for backup/export to other
% software, plotting, analysis tools, and more.
%
% Input ::
%
%   StudyTable          - n x d table object of text and numeric data to
%                           store in dataset object.
%
%   'labels'            - (optional arg) followed by a n x 1 vector of
%                           subject labels. If provided DAT.Event_Level
%                           will be populated with data from StudyTable as
%                           well, with one cell per subject.
%
% Output ::
%
%   DAT                 - canlab_dataset object
%
function DAT = table2canlab_dataset(StudyTable, varargin)

labels = [];

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case 'labels'
                if length(varargin{i+1}) ~= size(StudyTable,1);
                    error('''labels'' length does not match StudyTable rows.');
                else
                    labels = varargin{i+1};
                end
            otherwise
                warning('Did not understand argument %s. Ignoring.',varargin{i});
        end
    end
end

% Create canlab_dataset object
% -------------------------------------------------------------------------

DAT = canlab_dataset;

% Fill in Experiment-level information
% ------------------------------------------------------------------------

DAT.Description.Experiment_Name = [];
DAT.Description.Missing_Values = NaN;

DAT.Description.Subj_Level = {};

DAT.Description.Subj_Level{1} = 'Imported from structure' ;

% Names of subject-level data we want to store
% -------------------------------------------------------------------------

DAT.Subj_Level.names = StudyTable.Properties.VariableNames;

% Data from table
% -------------------------------------------------------------------------

studyCell = table2cell(StudyTable);
studyRow = studyCell(1,:);
numeric = ones(size(studyRow));
character = ones(size(studyRow));
for i = 1:length(studyRow)
    if ~isnumeric(studyRow{i});
        numeric(i) = 0;
    end
    if ~ischar(studyRow{i})
        character(i) = 0;
    end
end

if any(numeric==0) % upload numeric data to rows where available, nans elsewhere, upload text to textdata.
    warning(sprintf('Table is a cell array. Normally entire table gets assigned to Subj_Level.data. Will additionally assign text data to Subj_Level.textdata. Unexpected things may happen downstream.'));
    if sum(numeric) + sum(character) ~= size(studyRow,2)
        error('Non numeric and non character data is no supported by table2canlab_dataset. Please remove nontext and nonnumeric classes from table and try again.');
    end
    DAT.Subj_Level.data = zeros(size(studyCell));
    DAT.Subj_Level.data(:) = nan;
    DAT.Subj_Level.data(:,logical(numeric)) =  cell2mat(studyCell(:,logical(numeric)));

    DAT.Subj_Level.textdata = cell(size(studyCell));
    DAT.Subj_Level.textdata(:,logical(character)) = studyCell(:,logical(character));
else
    DAT.Subj_Level.data = table2array(StudyTable);
end
DAT.Subj_Level.descrip{1} = 'Imported from Table object. Add descriptions.';
DAT.Subj_Level.type = repmat({'numeric'}, 1, length(DAT.Subj_Level.names));

% Subject ID - put in cell array of strings

DAT.Subj_Level.id = {};
nobs = size(DAT.Subj_Level.data, 1);

for i = 1:nobs
    DAT.Subj_Level.id{i} = sprintf('Obs%3.0f', i);
end

%% populate Event_Level property
if ~isempty(labels)
    DAT.Event_Level.names = DAT.Subj_Level.names;
    DAT.Event_Level.descrip = 'Imported from table.';
    DAT.Event_Level.units = repmat({'Imported from table.'},1,length(DAT.Event_Level.names));
    DAT.Event_Level.type = cell(1,length(numeric));
    DAT.Event_Level.type(logical(numeric)) = repmat({class(DAT.Subj_Level.data)},1,sum(numeric));
    DAT.Event_Level.type(logical(character)) = repmat({'char'},1,sum(character));
    uniq_lbls = unique(labels);
    n = length(uniq_lbls);
    for i = 1:n
        this_lbl = uniq_lbls(i);
        lbl_idx = find(labels == this_lbl);
        DAT.Event_Level.data{end+1} = DAT.Subj_Level.data(lbl_idx,:);
        if ~isempty(DAT.Subj_Level.textdata)
            DAT.Event_Level.textdata{end+1} = DAT.Subj_Level.textdata(lbl_idx,:);
        end
    end
end

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function g = negNoiseGradientParam(params, model, prior)

% NEGNOISEGRADIENTPARAM Wrapper function for calling noise gradients.

% NOISE

model.noise = noiseExpandParam(model.noise, params);
fhandle = str2func([model.noise.type 'GradientParam']);
g = - fhandle(model);

if prior
  g = g + params;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function f=vimpelplot(amplitude)

% Uppbyggnad och utritning av vimpel.
global vimpel

i=[-2 -1 1 2];
sinstartx=3.545; sinstopx=6;
Xv=[]; Yv=[]; Zv=[];
for x=1:0.05:2
  scaledx=(sinstopx-sinstartx)*(x-1)+sinstartx;
  Xv=[Xv;
      x*ones(size(i))];
  Yv=[Yv;
      amplitude*sin(0.5*scaledx^2)/(scaledx)*ones(size(i))];
  Zv=[2.6-i*(2-x)/20
      Zv];
end
Cv=zeros(size(Zv))+1;
Cv(:,2)=Cv(:,2)*0;
vimpel=surf(Xv,Yv,Zv,Cv);

author: Sable

repo: mcbench-benchmarks


text:
function element_size = p00_element_size ( test )

%*****************************************************************************80
%
%% P00_ELEMENT_SIZE returns a typical element size for a problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TEST, the index of the test problem
%
%    Input, ELEMENT_SIZE, a typical element size.
%
  if ( test == 1 )
    element_size = p01_element_size ( );
  elseif ( test == 2 )
    element_size = p02_element_size ( );
  elseif ( test == 3 )
    element_size = p03_element_size ( );
  elseif ( test == 4 )
    element_size = p04_element_size ( );
  elseif ( test == 5 )
    element_size = p05_element_size ( );
  elseif ( test == 6 )
    element_size = p06_element_size ( );
  elseif ( test == 7 )
    element_size = p07_element_size ( );
  elseif ( test == 8 )
    element_size = p08_element_size ( );
  elseif ( test == 9 )
    element_size = p09_element_size ( );
  elseif ( test == 10 )
    element_size = p10_element_size ( );
  elseif ( test == 11 )
    element_size = p11_element_size ( );
  elseif ( test == 12 )
    element_size = p12_element_size ( );
  elseif ( test == 13 )
    element_size = p13_element_size ( );
  elseif ( test == 14 )
    element_size = p14_element_size ( );
  elseif ( test == 15 )
    element_size = p15_element_size ( );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_ELEMENT_SIZE - Fatal error!\n' );
    fprintf ( 1, '  Input value of TEST is out of range.\n' );
    error ( 'P00_ELEMENT_SIZE - Fatal error!' );
  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Setting up advection-diffusion solver
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C = please_Update_Adv_Diff_Concentration_Split(C,dt,dx,dy,uX,uY,k)

% C:     concentration
% dt:    time-step
% dx,dy: spatial steps in x and y, respectively
% uX:    x-Component of Velocity
% uY:    y-Component of Velocity
% k:     diffusion coefficient

% Compute Necessary Derivatives for x-Advection
Cx = give_Necessary_Derivative(C,dx,uX,'x');
Cxx = DD(C,dx,'x');

% Time-step #1 (give auxillary)
C = C + dt * ( k*(Cxx) - uX.*Cx );

% Compute Necessary Derivatives for y-Advection
Cy = give_Necessary_Derivative(C,dy,uY,'y');
Cyy = DD(C,dy,'y');

% Time-step #2 (give next iteration)
C = C + dt * ( k*(Cyy) - uY.*Cy );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Computes derivative based on sign of Velocity, u, using UPWIND
% approach
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C_z = give_Necessary_Derivative(C,dz,uZ,string)

C_z = zeros(size(C));
len = length(uZ(:,1));
signs = sign(uZ);

if strcmp(string,'x')

    %For periodicity on ends w/ UPWIND
    for i=1:len

        %left side of grid
        if signs(i,1) <= 0
            C_z(i,1) =  ( C(i,2) - C(i,1) ) / (dz);
        else
            C_z(i,1) =  ( C(i,1) - C(i,len) ) / (dz);
        end

        %right side of grid
        if signs(i,len) <= 0
            C_z(i,len) =  ( C(i,1) - C(i,len) ) / (dz);
        else
            C_z(i,len) =  ( C(i,len) - C(i,len-1) ) / (dz);
        end

    end

    %Standard Upwind
    for i=1:len
        for j=2:len-1
            if signs(i,j) <= 0
                C_z(i,j) = ( C(i,j+1) - C(i,j) ) / (dz);
            else
                C_z(i,j) = ( C(i,j) - C(i,j-1) ) / (dz);
            end
        end
    end

    % Ends x-Direction calculation %

elseif strcmp(string,'y')

    %For periodicity on ends w/ UPWIND
    for i=1:len

        %bottom of grid
        if signs(1,i) <= 0
            C_z(1,i) =  ( C(2,i) - C(1,i) ) / (dz);
        else
            C_z(1,i) =  ( C(1,i) - C(len,i) ) / (dz);
        end

        %top of grid
        if signs(len,i) <= 0
            C_z(len,i) =  ( C(1,i) - C(len,i) ) / (dz);
        else
            C_z(len,i) =  ( C(len,i) - C(len-1,i) ) / (dz);
        end

    end

    %Standard Upwind
    for i=2:len-1
        for j=1:len
            if signs(i,j) <= 0
                C_z(i,j) = ( C(i+1,j) - C(i,j) ) / (dz);
            else
                C_z(i,j) = ( C(i,j) - C(i-1,j) ) / (dz);
            end
        end
    end

    % Ends y-Direction calculation %

else

    fprintf('\n\n\n ERROR IN FUNCTION D FOR COMPUTING 1ST DERIVATIVE\n');
    fprintf('Need to specify which desired derivative, x or y.\n\n\n');

end

clear signs; clear len;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Finds CENTERED finite difference approximation to 2ND
% DERIVATIVE in z direction, specified by input and 'string'
% Note: It automatically accounts for periodicity of the domain.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u_zz = DD(u,dz,string)

% u:      velocity
% dz:     spatial step in "z"-direction
% string: specifies which 2ND derivative to take (to enforce periodicity)

len = length(u(:,1));

if strcmp(string,'x')

    %For periodicity on ends
    u_zz(:,1) =  ( u(:,2) - 2*u(:,1)   + u(:,len) )   / (dz^2);
    u_zz(:,len)= ( u(:,1) - 2*u(:,len) + u(:,len-1) ) / (dz^2);

    %Standard Upwind Scheme (Centered Difference)
    for j=2:len-1
        u_zz(:,j) = ( u(:,j+1) - 2*u(:,j) + u(:,j-1) ) / (dz^2);
    end

elseif strcmp(string,'y')

    %For periodicity on ends
    u_zz(1,:) =  ( u(2,:) - 2*u(1,:)   + u(len,:) )   / (dz^2);
    u_zz(len,:)= ( u(1,:) - 2*u(len,:) + u(len-1,:) ) / (dz^2);

    %Standard Upwind Scheme (Centered Difference)
    for j=2:len-1
        u_zz(j,:) = ( u(j+1,:) - 2*u(j,:) + u(j-1,:) ) / (dz^2);
    end

else

    fprintf('\n\n\n ERROR IN FUNCTION DD FOR COMPUTING 2ND DERIVATIVE\n');
    fprintf('Need to specify which desired derivative, x or y.\n\n\n');

end

author: nickabattista


text:
%% Example of using KST class for interfacing with KUKA iiwa robots

% Controlling KUKA iiwa robot using 3D mouse.

% Important note: when operating the robot with 3D mouse
% do not drive the robot near joints limits
% or near work space limits, the author is not liable for any sort of
% injury or loss of property due to the use of the provided programs.

% This example works with Sunrise application version KST_1.7  and higher.
% Copyright Mohammad SAFEEA, 26th-June-2018

close all;clear;clc;
instrreset;
warning('off')
%% Create the robot object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_elektrisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
return;
end
pause(1);
disp('Control iiwa using 3D mouse')
disp('Do not move the robot near joint limits')
disp('Do not move the robot near work space limits')

% Initial configuration
jPos={0.0, pi / 180 * 30, 0, -pi / 180 * 60, 0,...
                    pi / 180 * 90, 0};
relVel=0.15;
iiwa.movePTPJointSpace(jPos, relVel); % move to initial configuration

% Start direct servo in joint space
iiwa.realTime_startDirectServoJoints();

% IK solver parameters
numberOfIterations=10;
lambda=0.1;
%% Tool transform matrix
TefTool=eye(4);
TefTool(1,4)=0.0;
TefTool(2,4)=0.0;
TefTool(3,4)=0.04;

%% initial configuration vector
qin=zeros(7,1);
daINDEX=0;
for i=1:7
qin(i)=jPos{i};
end

%% setup the mouse
% setup the mouse, normally it should be in USB1
ID = 'USB1';
MOUSE = [];
try
% try to create the space mouse object MOUSE
MOUSE = vrspacemouse(ID);
catch ME
fprintf('Unable to initialize the Space Mouse on port %s.\n', ID);
end

MOUSE.PositionSensitivity = 1e-2*10;
MOUSE.RotationSensitivity = 1e-5*10000;
vector = [1;2;3;4;5;6];

%% Variables Initiation
[Tt,j]=iiwa.directKinematics(qin); % EEF frame transformation amtrix
filterVar=0.2;
v_control=zeros(3,1);
w_control=zeros(3,1);
firstExecution=0;
linearVelocityScale=4500;
angularVelocityScale=4000;
% Transform the measurments, change it accouding to your own 3D mouse
% configuration
darV=[1 0 0;
0 -1 0;
0 0 -1];

darW=[0 -1 0;
1 0 0;
0 0 0];

pause(0.1);

disp('Control is activeted');
disp('Use space mouse to control the robot');
%% Control loop
while true

% Calculate the elapsed time
if(firstExecution==0)
    firstExecution=1;
    a=datevec(now);
    timeNow=a(6)+a(5)*60+a(4)*60*60; % calculate time at this instant
    dt=0; % elapsed is zero at first excution
    time0=timeNow;
else
    a=datevec(now);
    timeNow=a(6)+a(5)*60+a(4)*60*60; % calculate time at this instant
    dt=timeNow-time0; % calculate elapsed time interval
    time0=timeNow;
end

% Read input data from mouse
 mouseSpeed = speed(MOUSE, vector);

%% Calculate the cartesian motion command
%------------------------------------------------------------------------
normalizedAngularMotionInput=[mouseSpeed(4);
     mouseSpeed(5);
     mouseSpeed(6)]/angularVelocityScale; % normalized angular motion command,
% calculated from mouse measurments
normalizedLinearMotionInput=[mouseSpeed(1);
    mouseSpeed(2);
    mouseSpeed(3)]/linearVelocityScale;% normalized linear motion command,
% calculated from mouse measurments

w_control=filterVar*w_control+(1-filterVar)*darV*normalizedAngularMotionInput;
v_control=filterVar*v_control+(1-filterVar)*darV*normalizedLinearMotionInput;

% Calculate target transform using command, and current transform
[Tt,notImportant]=iiwa.directKinematics(qin);
Tt=updateTransform2(Tt,w_control,v_control,dt);

qt= iiwa.gen_InverseKinematics(qin, Tt,numberOfIterations,lambda );
% A condition to break the loop, when first space mouse button is
% pressed
if button(MOUSE, 1) == 1
    break;
end

for i=1:7
    jPos{i}=qt(i);
end
%% Send reference joints to robot
iiwa.sendJointsPositions(jPos);
%% update initial positions, target transform
qin=qt;
end

%% Turn off the direct servo
iiwa.realTime_stopDirectServoJoints();

%% turn Off blue light
iiwa.net_turnOffServer();

repo: KST-Kuka-Sunrise-Toolbox