text:
function output = ...
  EstimateInitialF0AndDeviation(probability_structure, trajectory, option)
% estimate F0 initial value from instantaneous frequency map, residual map
% and given trajectory
%
% output = ...
%   EstimateInitialF0AndDeviation(probability_structure, trajectory)
%
% output = ...
%   EstimateInitialF0AndDeviation(probability_structure, trajectory, option)
%
% Arguments
%
%   probability_structure: structure variable, following fields are used
%     inst_freq_map: matrix, time-frequency map of instantaneous frequency
%     residual_map: matrix, time-frequency map of residuals
%     center_frequency_list: vector, list of filter center frequency
%     channels_in_octave: scalar, number of channels in one octave
%  trajectory: vector: sequence of F0 trajectory in channel ID number
%  option: structure variable with the following fields
%    residual_margin: scalar, safeguard for residual selection
%    frequency_margin: scalar, safeguard for limiting instantaneous
%      frequency
%
% Return value
%
%  output: structure variable with the following fields
%    f0_initial: vector, initial estimate of F0 sequence (Hz)
%    sd_initial: vector, initial estimate of relative deviation
%    temporal_positions: vector, frame center location (s)
%    residual_margin: scalar, actual value used
%    frequency_margin: scalar, actual value used
%    sd_on_trajectory: vector, residual value on give trajectory
%    elapsed_time: scalar, time needed to process this function

% Copyright 2016 Google Inc. All Rights Reserved
% Author: hidekik@google.com (Hideki Kawahara)

start_tic = tic;
residual_margin = 1.5; % default value
frequency_margin = 1.15; % default value
narginchk(2, 3);
if nargin == 3
  if isfield(option, 'residual_margin')
    residual_margin = option.residual_margin;
  end;
  if isfield(option, 'frequency_margin')
    frequency_margin = option.frequency_margin;
  end;
end;
inst_freq_map = probability_structure.inst_freq_map;
residual_map = probability_structure.residual_map;
fc_list = probability_structure.center_frequency_list;
n_channels = size(inst_freq_map,1);
n_frames = size(inst_freq_map, 2);
f0_initial = zeros(n_frames, 1);
sd_initial = zeros(n_frames, 1);
sample_f0_floor = ...
  probability_structure.center_frequency_list(min(trajectory));
sample_f0_ceiling = ...
  probability_structure.center_frequency_list(max(trajectory));
sd_on_trajectory = zeros(n_frames, 1); % for debug
n_in_octave = probability_structure.channels_in_octave;
scan_base = -ceil(n_in_octave / 2 + 2):ceil(n_in_octave / 4);
for ii = 1:n_frames
  base_index = trajectory(ii) + scan_base;
  base_index = base_index(base_index > 0 & base_index <= n_channels);
  reference_residual = residual_map(trajectory(ii), ii);
  sd_on_trajectory(ii) = sqrt(reference_residual);
  base_index = base_index(residual_map(base_index, ii) < ...
    reference_residual * residual_margin);
  reference_if = fc_list(trajectory(ii));
  tmp_residual = residual_map(base_index, ii);
  tmp_if = inst_freq_map(base_index, ii);
  tmp_if = max(reference_if / frequency_margin, ...
    min(reference_if * frequency_margin, tmp_if));
  if length(base_index) > 1
    best_weight = GetBestMixingWeight(sqrt(tmp_residual));
  else
    best_weight = 1;
  end;
  f0_initial(ii) = sum(best_weight .* tmp_if);
  sd_initial(ii) = sqrt(sum(best_weight.^2 .* tmp_residual));
end;
output = ...
  struct('f0_initial', max(sample_f0_floor, ...
                           min(sample_f0_ceiling, f0_initial)), ...
         'sd_initial', sd_initial, ...
         'temporal_positions', probability_structure.temporal_positions, ...
         'sample_f0_floor', sample_f0_floor, ...
         'sample_f0_ceiling', sample_f0_ceiling, ...
         'residual_margin', residual_margin, ...
         'frequency_margin', frequency_margin, ...
         'sd_on_trajectory', sd_on_trajectory, ...
         'elapsed_time', toc(start_tic));
end

author: google


text:
function [Y, lbls, Ytest, lblstest] = lvmLoadData(dataset, seedVal)

% LVMLOADDATA Load a latent variable model dataset.
% FORMAT
% DESC loads a data set for a latent variable modelling problem.
% ARG dataset : the name of the data set to be loaded. Currently
% the possible names are 'robotWireless',
% 'robotTwoLoops', 'robotTraces', 'robotTracesTest', 'cmu35gplvm',
% 'cmu35Taylor', 'cmu35walkJog', 'vowels', 'stick', 'brendan',
% 'digits', 'twos', 'oil', 'oilTest', 'oilValid', 'oil100',
% 'swissRoll', 'missa'.
% RETURN Y : the training data loaded in.
% RETURN lbls : a set of labels for the data (if there are no
% labels it is empty).
% RETURN Ytest : the test data loaded in. If no test set is
% available it is empty.
% RETURN lblstest : a set of labels for the test data (if there are
% no labels it is empty).
%
% SEEALSO : mapLoadData, datasetsDirectory
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2008, 2009

% DATASETS

  if nargin > 1
    randn('seed', seedVal)
    rand('seed', seedVal)
  end

  % get directory

  baseDir = datasetsDirectory;
  dirSep = filesep;
  lbls = [];
  lblstest = [];
  switch dataset
   case 'YaleSubset6_1'
          load([baseDir 'YaleSubset6.mat']);
          lbls = [height;width];
   case 'missa'
       load([baseDir 'miss-americaHD.mat']);
       lbls = [288 ;360];
   case 'movielens'
    try
      load([baseDir 'movielens.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');
        fileName = [baseDir dirSep 'movielens' dirSep 'large' dirSep 'ratings.dat'];
        [users, films, ratings, timeStamp] = textread(fileName, '%n::%n::%n::%n');
        ind = randperm(size(users, 1));
        users = users(ind, :);
        films = films(ind, :);
        ratings = ratings(ind, :);
        numUsers = max(users);
        numFilms = max(films);

        numRatings = size(users, 1);
        numTrainRatings = ceil(0.8*numRatings);
        Y = spalloc(numFilms, numUsers, numTrainRatings);
        Ytest = spalloc(numFilms, numUsers, numRatings-numTrainRatings);
        indTrain = sub2ind(size(Y), films(1:numTrainRatings), users(1:numTrainRatings));
        indTest = sub2ind(size(Ytest), films(numTrainRatings+1:numRatings), users(numTrainRatings+1:numRatings));
        Y(indTrain) = ratings(1:numTrainRatings);
        Ytest(indTest) = ratings(numTrainRatings+1:numRatings);
        save([baseDir 'movielens.mat'], 'Y', 'lbls', 'Ytest', 'lblstest');
      else
        error(lasterr);
      end
    end

   case {'movielensSmall1', 'movielensSmall2', 'movielensSmall3', 'movielensSmall4', 'movielensSmall5'}

    partNo = str2num(dataset(end));
    uTrain = load([baseDir dirSep 'movielens' dirSep 'small' dirSep 'u' num2str(partNo) '.base']);
    numUsers = max(uTrain(:, 1));
    numFilms = max(uTrain(:, 2));
    numRatings = size(uTrain, 1);
    Y = spalloc(numFilms, numUsers, numRatings);

    for i = 1:size(uTrain, 1);
      Y(uTrain(i, 2), uTrain(i, 1)) = uTrain(i, 3);
    end
    meanY = mean(Y(find(Y)));
    Y(find(Y)) = (Y(find(Y))-meanY);
    uTest = load([baseDir dirSep 'movielens' dirSep 'small' dirSep 'u' num2str(partNo) '.test']);
    numTestRatings = size(uTest, 1);
    Ytest = spalloc(numFilms, numUsers, numTestRatings);
    for i = 1:size(uTest, 1);
      Ytest(uTest(i, 2), uTest(i, 1)) = uTest(i, 3);
    end
    Ytest(find(Ytest)) = (Ytest(find(Ytest))-meanY);
    stdY = sqrt(var(Y(find(Y))));
    Y = Y/stdY;
    Ytest = Ytest/stdY;

   case 'crabs'
    fid = fopen([baseDir 'crabs.dat']);
    if fid == -1
      error(['No such file name ' fileName])
    end
    readLine = fgets(fid);
    counter = 0;
    data = [];
    while readLine ~= -1
      readLine = fgets(fid);
      counter = counter + 1;
      parts = stringSplit(readLine, ' ');
      if parts{1} == 'B' && parts{2} == 'F'
        lbls(counter, :) = [0 0 0 1];
      elseif parts{1} == 'B' && parts{2} == 'M'
        lbls(counter, :) = [0 0 1 0];
      elseif parts{1} == 'O' && parts{2} == 'F'
        lbls(counter, :) = [0 1 0 0];
      else
        lbls(counter, :) = [1 0 0 0];
      end
      counter2 = 0;
      for i = 3:size(parts, 2)
        if ~isempty(str2num(parts{i}))
          counter2 = counter2 + 1;
          data(counter, counter2) = str2num(parts{i});
        end
      end
    end
    fclose(fid);
    data = data(:, 2:end);
    Y = data./repmat(sum(data, 2), 1, size(data, 2));

   case 'spellman'
    try
      load([baseDir 'spellman.mat']);

    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');
        fileName = ([baseDir 'combined.csv']);
        [geneName, ...
         y1, y2, y3, y4, y5, y6, y7, y8, y9, ...
         y10, y11, y12, y13, y14, y15, y16, y17, y18, ...
         y19, y20, y21, y22, y23, y24, y25, y26, y27, ...
         y28, y29, y30, y31, y32, y33, y34, y35, y36, ...
         y37, y38, y39, y40, y41, y42, y43, y44, y45, ...
         y46, y47, y48, y49, y50, y51, y52, y53, y54, ...
         y55, y56, y57, y58, y59, y60, y61, y62, y63, ...
         y64, y65, y66, y67, y68, y69, y70, y71, y72, ...
         y73, y74, y75, y76, y77] = ...
            textread(fileName, ['%q ' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f %f %f %f %f' ...
                            '%f %f %f %f %f'], ...
                     'headerlines', 1, 'whitespace', '', 'delimiter', ',');

        Y = [y1 y2 y3 y4 y5 y6 y7 y8 y9 ...
             y10, y11, y12, y13, y14, y15, y16, y17, y18, ...
             y19, y20, y21, y22, y23, y24, y25, y26, y27, ...
             y28, y29, y30, y31, y32, y33, y34, y35, y36, ...
             y37, y38, y39, y40, y41, y42, y43, y44, y45, ...
             y46, y47, y48, y49, y50, y51, y52, y53, y54, ...
             y55, y56, y57, y58, y59, y60, y61, y62, y63, ...
             y64, y65, y66, y67, y68, y69, y70, y71, y72, ...
             y73, y74, y75, y76, y77]';
        lbls = [repmat([1 0 0 0 0], 4, 1);
                repmat([0 1 0 0 0], 18,1);
                repmat([0 0 1 0 0], 24,1);
                repmat([0 0 0 1 0], 17,1);
                repmat([0 0 0 0 1], 14,1)];
        save([baseDir 'spellman.mat'], 'Y', 'lbls');
      else
        error(lasterr);
      end
    end
   case 'robotWireless'
    Ydat = parseWirelessData([baseDir 'uw-floor.txt']);
    Ydat = (Ydat + 85)/15;
    lbls = 'connect';
    Y = Ydat(1:215, :);
    Ytest = Ydat(216:end, :);

   case 'robotTwoLoops'
    Y = csvread([baseDir 'TwoLoops.slam'], 1, 0);
    Y = Y(1:floor(end/2), 4:end);
    Y(find(Y==-100))=-NaN;
    Y = (Y + 85)/15;
    lbls = 'connect';

   case 'robotTraces'
    Y = csvread([baseDir 'Trace-3rdFloor-01.uwar.slam'], 1, 0);
    Y = Y(1:floor(end/2), 4:end);
    Y(:, [3 4 38]) = []; % Remove columns of missing data.
    Y(find(Y==-100))=NaN;
    Y = (Y + 85)/15;
    lbls = 'connect';

   case 'robotTracesTest'
    Y = csvread([baseDir 'Trace-3rdFloor-01.uwar.slam'], 1, 0);
    Y = Y(ceil(end/2):end, 4:end);
    Y(:, [3 4 38]) = []; % Remove columns of missing data.
    Y(find(Y==-100))=NaN;
    Y = (Y + 85)/15;
    lbls = 'connect';

   case 'osr'
    t = linspace(-pi, pi, 256);
    window = (cos(t) + 1)/2;
    [x, fs] = wavread([baseDir 'osr' dirSep 'OSR_uk_000_0020_8k.wav']);
    [r, theta] = stft(x, fs, 256, window);
    Y = [r(:, 1:32) theta(:, 1:32)];

   case 'cmu49BalanceArm'
    [Y, lbls, Ytest, lblstest] = lvmLoadData('cmu49Balance');
    % Extract left arm only
    ind = [41:47 49:50];
    Y = Y(:, ind);
    Ytest = Ytest(:, ind);

   case 'cmu49Balance'

    try
      load([baseDir 'cmu49Balance.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');
        sampleEvery = 32;
        skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '49' dirSep '49.asf']);
        examples = {'18', '19'};
        testExamples = {'20'};
        % Label differently for each sequence
        exlbls = eye(length(examples));
        testexlbls = eye(length(testExamples));
        totLength = 0;
        totTestLength = 0;
        for i = 1:length(examples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '49' dirSep '49_' ...
                              examples{i} '.amc'], skel);
          tY{i} = tmpchan(1:sampleEvery:end, :);
          tlbls{i} = repmat(exlbls(i, :), size(tY{i}, 1), 1);
          totLength = totLength + size(tY{i}, 1);
        end
        Y = zeros(totLength, size(tY{1}, 2));
        lbls = zeros(totLength, size(tlbls{1}, 2));
        endInd = 0;
        for i = 1:length(tY)
          startInd = endInd + 1;
          endInd = endInd + size(tY{i}, 1);
          Y(startInd:endInd, :) = tY{i};
          lbls(startInd:endInd, :) = tlbls{i};
        end
        for i = 1:length(testExamples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '49' dirSep '49_' ...
                              testExamples{i} '.amc'], skel);
          tYtest{i} = tmpchan(1:sampleEvery:end, :);
          tlblstest{i} = repmat(testexlbls(i, :), size(tYtest{i}, 1), 1);
          totTestLength = totTestLength + size(tYtest{i}, 1);
        end
        Ytest = zeros(totTestLength, size(tYtest{1}, 2));
        lblstest = zeros(totTestLength, size(tlblstest{1}, 2));
        endInd = 0;
        for i = 1:length(tYtest)
          startInd = endInd + 1;
          endInd = endInd + size(tYtest{i}, 1);
          Ytest(startInd:endInd, :) = tYtest{i};
          lblstest(startInd:endInd, :) = tlblstest{i};
        end
        save([baseDir 'cmu49Balance.mat'], 'Y', 'lbls', 'Ytest', 'lblstest');
      else
        error(lasterr);
      end
    end

   case 'cmu07WalkFeet'
    [Y, lbls, Ytest, lblstest] = lvmLoadData('cmu07Walk');
    % Extract left arm only
    ind = [7:12 14:19];
    Y = Y(:, ind);
    Ytest = Ytest(:, ind);

   case 'cmu07Walk'

    try
      load([baseDir 'cmu07Walk.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');
        sampleEvery = 4;
        skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '07' dirSep '07.asf']);
        examples = {'01', '02'};
        testExamples = {'03'};
        % Label differently for each sequence
        exlbls = eye(length(examples));
        testexlbls = eye(length(testExamples));
        totLength = 0;
        totTestLength = 0;
        for i = 1:length(examples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '07' dirSep '07_' ...
                              examples{i} '.amc'], skel);
          tY{i} = tmpchan(1:sampleEvery:end, :);
          tlbls{i} = repmat(exlbls(i, :), size(tY{i}, 1), 1);
          totLength = totLength + size(tY{i}, 1);
        end
        Y = zeros(totLength, size(tY{1}, 2));
        lbls = zeros(totLength, size(tlbls{1}, 2));
        endInd = 0;
        for i = 1:length(tY)
          startInd = endInd + 1;
          endInd = endInd + size(tY{i}, 1);
          Y(startInd:endInd, :) = tY{i};
          lbls(startInd:endInd, :) = tlbls{i};
        end
        for i = 1:length(testExamples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '07' dirSep '07_' ...
                              testExamples{i} '.amc'], skel);
          tYtest{i} = tmpchan(1:sampleEvery:end, :);
          tlblstest{i} = repmat(testexlbls(i, :), size(tYtest{i}, 1), 1);
          totTestLength = totTestLength + size(tYtest{i}, 1);
        end
        Ytest = zeros(totTestLength, size(tYtest{1}, 2));
        lblstest = zeros(totTestLength, size(tlblstest{1}, 2));
        endInd = 0;
        for i = 1:length(tYtest)
          startInd = endInd + 1;
          endInd = endInd + size(tYtest{i}, 1);
          Ytest(startInd:endInd, :) = tYtest{i};
          lblstest(startInd:endInd, :) = tlblstest{i};
        end
        save([baseDir 'cmu07Walk.mat'], 'Y', 'lbls', 'Ytest', 'lblstest');
      else
        error(lasterr);
      end
    end

   case 'cmu35gplvm'
    [Y, lbls, Ytest, lblstest] = lvmLoadData('cmu35WalkJog');
    skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35.asf']);
    [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35_01.amc'], skel);

    Ytest = Ytest(find(lblstest(:, 2)), :);
    lblstest = lblstest(find(lblstest(:, 2)), 2);

    %left indices
    xyzInd = [2];
    xyzDiffInd = [1 3];
    rotInd = [4 6];
    rotDiffInd = [5];
    generalInd = [7:38 41:47 49:50 53:59 61:62];

    jointAngles  = asin(sin(pi*Y(:, generalInd)/180));
    jointAnglesTest  = asin(sin(pi*Ytest(:, generalInd)/180));

    endInd = [];
    for i = 1:size(lbls, 2)
      endInd = [endInd max(find(lbls(:, i)))];
    end
    catJointAngles = [];
    xyzDiff = [];
    catSinCos = [];
    startInd = 1;
    for i = 1:length(endInd)
      ind1 = startInd:endInd(i)-1;
      ind2 = startInd+1:endInd(i);
      catJointAngles = [catJointAngles; ...
                        jointAngles(ind2, :)];
      xyzDiff = [xyzDiff;
                 Y(ind1, xyzDiffInd) - Y(ind2, xyzDiffInd) ...
                 Y(ind2, xyzInd)];
      catSinCos = [catSinCos; ...
                   sin(pi*Y(ind2, rotInd)/180) ...
                   sin(pi*Y(ind1, rotDiffInd)/180)-sin(pi*Y(ind2, rotDiffInd)/180) ...
                   cos(pi*Y(ind2, rotInd)/180) ...
                   cos(pi*Y(ind1, rotDiffInd)/180)-cos(pi*Y(ind2, rotDiffInd)/180)];
      startInd = endInd(i)+1;
    end
    Y = [catJointAngles xyzDiff catSinCos];
    lbls = [];

    endInd = [];
    for i = 1:size(lblstest, 2)
      endInd = [endInd max(find(lblstest(:, i)))];
    end
    catJointAnglesTest = [];
    xyzDiffTest = [];
    catSinCosTest = [];
    startInd = 1;
    for i = 1:length(endInd)
      ind1 = startInd:endInd(i)-1;
      ind2 = startInd+1:endInd(i);
      catJointAnglesTest = [catJointAnglesTest; ...
                          jointAnglesTest(ind2, :)];
      xyzDiffTest = [xyzDiffTest;
                     Ytest(ind1, xyzDiffInd) - Ytest(ind2, xyzDiffInd) ...
                     Ytest(ind2, xyzInd)];
      catSinCosTest = [catSinCosTest; ...
                       sin(pi*Ytest(ind2, rotInd)/180) ...
                       sin(pi*Ytest(ind1, rotDiffInd)/180)-sin(pi*Ytest(ind2, rotDiffInd)/180) ...
                       cos(pi*Ytest(ind2, rotInd)/180) ...
                       cos(pi*Ytest(ind1, rotDiffInd)/180)-cos(pi*Ytest(ind2, rotDiffInd)/180)];
      startInd = endInd(i)+1;
    end
    Ytest = [catJointAnglesTest xyzDiffTest catSinCosTest];
    lblstest = [];

   case 'cmu35Taylor'
    % The CMU 35 data set as Graham Taylor et al. have it in their
    % NIPS 2006 paper.
    [rawY, lbls, rawYtest, lblstest] = lvmLoadData('cmu35WalkJog');
    rawYtest = rawYtest(45:end, :);
    lblstest = lblstest(45:end, 2);
    skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35.asf']);
    [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35_01.amc'], skel);
    posInd = [1 2 3];
    rotInd = [4 5 6];
    generalInd = [7:38 41:47 49:50 53:59 61:62];

    rawY = [rawY(:, posInd) sind(rawY(:, rotInd)) cosd(rawY(:, rotInd)) asind(sind(rawY(:, generalInd)))];
    rawYtest = [rawYtest(:, posInd) sind(rawYtest(:, rotInd)) cosd(rawYtest(:, rotInd)) asind(sind(rawYtest(:, generalInd)))];

    endInd = [];
    for i = 1:size(lbls, 2)
      endInd = [endInd max(find(lbls(:, i)))];
    end
    Y = [];
    startInd = 1;
    for i = 1:length(endInd)
      ind1 = startInd:endInd(i)-1;
      ind2 = startInd+1:endInd(i);
      Y = [Y; rawY(ind1, :) rawY(ind2, :)];
      startInd = endInd(i)+1;
    end
    lbls = [];

    endInd = [];
    for i = 1:size(lblstest, 2)
      endInd = [endInd max(find(lblstest(:, i)))];
    end
    Ytest = [];
    startInd = 1;
    for i = 1:length(endInd)
      ind1 = startInd:endInd(i)-1;
      ind2 = startInd+1:endInd(i);
      Ytest = [Ytest; rawYtest(ind1, :) rawYtest(ind2, :)];
      startInd = endInd(i)+1;
    end
    lblstest = [];

   case {'dur', 'cmp'}
    [Y, void, lbls] = synthLoadData(dataset);

   case {'cmpdur'}
    [Y1, void, lbls] = synthLoadData('cmp');
    [Y2, void, lbls] = synthLoadData('dur');
    Y = [Y1, Y2];

   case {'grid_vowels', ...
         'grid_consonants', ...
         'grid_approximants', ...
         'grid_uplosives', ...
         'grid_vplosives', ...
         'grid_plosives', ...
         'grid_nasals', ...
         'grid_ufricatives', ...
         'grid_vfricatives', ...
         'grid_fricatives', ...
         'grid_affricatives', ...
         'grid_silence', ...
         'grid_other', ...
         'grid_all'}
    try
      load([baseDir dataset '.mat']);
    catch

      fullList = {'iy2', 'iy4', 'iy5', 'iy6', ...     % Vowels (iy, ih, eh, ey, ae, aa, aw, ah, ao, ow, uw, ax, ia)
                'ih1', 'ih3', 'ih5', ...
                'eh1', 'eh2', 'eh4', 'eh5', 'eh6', ...
                'ey1', 'ey4', 'ey5', ...
                'ae3', ...
                'aa4', ...
                'aw6', ...
                'ay2', 'ay3', 'ay4', 'ay5', ...
                'ah5', ...
                'ao5', ...
                'ow4', 'ow5', ...
                'uw2', 'uw4', 'uw5', 'uw6', ...
                'ax4', 'ax6', ...
                'ia5', ...
                'w2', 'w3', 'w4', 'w5', ... % approximants -- kind of between a vowel and consonant (w, y, r, l)
                'y4', ...
                'r2', 'r5', ...
                'l1', 'l2', 'l4', 'l6', ...
                'b1', 'b2', 'b3', 'b4', ... % Consonants  - plosives, voiced (b, d, g)
                'd2', 'd4', ...
                'g2', 'g6', ...
                'p1', 'p4', 'p6', ... % Consonants - plosives, unvoiced (p, t, k)
                't1', 't2', 't3', 't4', 't5',...
                'k4', 'k5', ...
                'm4', ... % Consonants - nasals (m, n)
                'n1', 'n2', 'n3', 'n4', 'n5', 'n6', ...
                'f4', 'f5', ... % Consonants - fricatives, unvoiced (f, s, th, sh)
                's1', 's4', 's5', 's6', ...
                'th5', ...
                'v4', 'v5', ... % Consonant - fricatives, voiced (v, z, dh, zh)
                'z4', 'z5', 'z6', ...
                'dh3', ...
                ...%(but consonsant can also be grouped strong/weak:   strong = s, sh, z, zh;  weak=f,v,th,dh)
                'ch4', ...  % Consonant - affricative  (basically a blend between a plosive and a fricative) (ch  =  t + sh, jh = d + zh)
                'jh4', ...
                'sil', ... %silence
                'sp'}; % ask Jon

      vowelInd = 1:33;
      consonantInd = 34:84;
      approximantInd = 34:44;
      vplosivesInd = 45:52;
      uplosivesInd = 53:62;
      nasalsInd = 63:69;
      ufricativesInd = 70:76;
      vfricativesInd = 77:82;
      affricativesInd = 83:84;
      silenceInd = 85;
      otherInd = 86;
      constanants = {};
      all = {};
      switch dataset(6:end)
       case 'vowels'
        phoneList = fullList(vowelInd);
       case 'consonants'
        phoneList = fullList(consonantInd);
       case 'approximants'
        phoneList = fullList(approximantInd);
       case 'uplosives'
        phoneList = fullList(uplosivesInd);
       case 'vplosives'
        phoneList = fullList(vplosivesInd);
       case 'plosives'
        phoneList = fullList([uplosivesInd vplosivesInd]);
       case 'nasals'
        phoneList = fullList(nasalsInd);
       case 'ufricatives'
        phoneList = fullList(ufricativesInd);
       case 'vfricatives'
        phoneList = fullList(vfricativesInd);
       case 'fricatives'
        phoneList = fullList([ufricativesInd vfricativesInd]);
       case 'affricatives'
        phoneList = fullList(affricativesInd);
       case 'silence'
        phoneList = fullList(silenceInd);
       case 'other'
        phoneList = fullList(otherInd);
       case 'all'
        phoneList = fullList;
      end

      motherDir = '/local/data/datasets/synthesis/models/';
      fileLoc = '/ver3/cmp/monophoneC.mmf';
      numId = 34;
      i = 1;
      idStr = ['id' num2str(i)];
      fileName = [motherDir idStr fileLoc];
      [meanValCell, varValCell] = htkLoadMmf(phoneList, fileName);
      meanVals = [meanValCell{:, :}];
      varVals = [varValCell{:, :}];
      meanVals = zeros(numId, size(meanVals, 2));
      varVals = zeros(numId, size(varVals, 2));

      for i = 1:numId
        idStr = ['id' num2str(i)];
        fileName = [motherDir idStr fileLoc];
        [meanValCell, varValCell] = htkLoadMmf(phoneList, fileName);
        meanVals(i, :) = [meanValCell{:, :}];
        varVals(i, :) = [varValCell{:, :}];
      end
      meanVals = meanVals./sqrt(varVals);
      varVals = 0.5*log(varVals);
      meanMean = mean(meanVals);
      meanLogStd = mean(varVals);
      for i = 1:numId
        meanVals(i, :) = meanVals(i, :) - meanMean;
        varVals(i, :) = varVals(i, :) - meanLogStd;
      end
      varMeanVals = var(meanVals(:));
      varVarVals = var(varVals(:));
      Y = [meanVals/sqrt(varMeanVals) varVals/sqrt(varVarVals)];
      lbls = [1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, ...
              0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0]';
      lbls = [lbls ~lbls];
      save([baseDir dataset '.mat'], 'Y', 'lbls', 'meanMean', 'meanLogStd', 'varMeanVals', 'varVarVals');
    end
    if size(lbls, 2) == 1
      lbls = [lbls ~lbls];
      save([baseDir dataset '.mat'], 'Y', 'lbls', 'meanMean', 'meanLogStd', 'varMeanVals', 'varVarVals');
    end
    retLbls = cell(2,1);
    retLbls{1} = lbls;
    %/~retLbls{2} = {'JB', 'Martin Cooke', 'Stuart Cunningham', 'Sue Harding', ...
%                  'Jonny Laidler', 'Stuart Wriggley', 'Gillian', 'Mike Stannett', ...
%                  'John Karn', 'UG', 'Anna', 'UG', ...
%                  'UG', 'UG', 'Female', 'UG', ...
%                  'UG', 'UG',  'UG', 'UG', ...
%                  'UG', 'UG', 'UG', 'UG', ...
%                  'UG', 'James Carmichael', 'Rob Mill', 'Matt Gibson', ...
%                  'Sarah Simpson', 'Vinny Wan', 'Lucy', 'NDL', ...
%                  'UG', 'UG'};
                  %~/
                  retLbls{2} = {'JB', 'MPC', 'Stuart Cunningham', 'Sue Harding', ...
                                'Jonny Laidler', 'Stuart Wriggley', 'South Yorkshire G', 'MS', ...
                                'JK', 'UG 1', 'South Yorkshire A', 'UG 2', ...
                                'UG 3', 'UG 4', 'UG 5', 'UG 6', ...
                                'UG 7', 'UG 8',  'UG 9', 'UG 10', ...
                                'UG 11', 'UG 12', 'UG 13', 'UG14', ...
                                'UG', 'West Indies', 'Rob Mill', 'Scottish', ...
                                'Sarah Simpson', 'Vinny Wan', 'South Yorkshire L', 'NDL', ...
                                'UG 15', 'UG 16'};
    lbls = retLbls;

   case 'cmu35WalkJog'
    try
      load([baseDir 'cmu35WalkJog.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');
        skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35.asf']);
        examples = ...
            {'01', '02', '03', '04', '05', '06', '07', '08', '09', '10', ...
             '11', '12', '13', '14', '15', '16', '17', '19', '20', ...
             '21', '22', '23', '24', '25', '26', '28', '30', '31', '32', '33', '34'};
        testExamples = {'18', '29'};
        % Label differently for each sequence
        exlbls = eye(31);
        testexlbls = eye(2);
        totLength = 0;
        totTestLength = 0;
        for i = 1:length(examples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35_' ...
                              examples{i} '.amc'], skel);
          tY{i} = tmpchan(1:4:end, :);
          tlbls{i} = repmat(exlbls(i, :), size(tY{i}, 1), 1);
          totLength = totLength + size(tY{i}, 1);
        end
        Y = zeros(totLength, size(tY{1}, 2));
        lbls = zeros(totLength, size(tlbls{1}, 2));
        endInd = 0;
        for i = 1:length(tY)
          startInd = endInd + 1;
          endInd = endInd + size(tY{i}, 1);
          Y(startInd:endInd, :) = tY{i};
          lbls(startInd:endInd, :) = tlbls{i};
        end
        for i = 1:length(testExamples)
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '35' dirSep '35_' ...
                              testExamples{i} '.amc'], skel);
          tYtest{i} = tmpchan(1:4:end, :);
          tlblstest{i} = repmat(testexlbls(i, :), size(tYtest{i}, 1), 1);
          totTestLength = totTestLength + size(tYtest{i}, 1);
        end
        Ytest = zeros(totTestLength, size(tYtest{1}, 2));
        lblstest = zeros(totTestLength, size(tlblstest{1}, 2));
        endInd = 0;
        for i = 1:length(tYtest)
          startInd = endInd + 1;
          endInd = endInd + size(tYtest{i}, 1);
          Ytest(startInd:endInd, :) = tYtest{i};
          lblstest(startInd:endInd, :) = tlblstest{i};
        end
        save([baseDir 'cmu35WalkJog.mat'], 'Y', 'lbls', 'Ytest', 'lblstest');
      else
        error(lasterr);
      end
    end

   case 'vowels'
    load([baseDir 'jon_vowel_data']);
    Y = [a_raw; ae_raw; ao_raw; ...
         e_raw; i_raw; ibar_raw; ...
         o_raw; schwa_raw; u_raw];
    Y(:, [13 26]) = [];
    lbls = [];
    for i = 1:9
      lbl = zeros(1, 9);
      lbl(i) = 1;
      lbls = [lbls; repmat(lbl, size(a_raw, 1), 1)];
    end

   case 'stick'
    Y = mocapLoadTextData([baseDir 'run1']);
    Y = Y(1:4:end, :);
    lbls = 'connect';

   case 'brendan'
    load([baseDir 'frey_rawface.mat']);
    Y = double(ff)';

   case 'usps'

    % Fix seeds
    randn('seed', 1e5);
    rand('seed', 1e5);

    load([baseDir 'usps_train.mat']);
    [ALL_T, sortIndices] = sort(ALL_T);
    ALL_DATA = ALL_DATA(sortIndices(:), :);
    Y = ALL_DATA;
    lbls = zeros(size(ALL_T,1),10);
    for digit = 0:9;
      ind1 = min(find(ALL_T==digit));
      ind2 = max(find(ALL_T==digit));
      lbls(ind1:ind2,digit+1) = ones(size(ind1:ind2,2),1);
    end

    % Fix seeds
    randn('seed', 1e5);
    rand('seed', 1e5);

    load([baseDir 'usps_test.mat']);
    [ALL_T, sortIndices] = sort(ALL_T);
    ALL_DATA = ALL_DATA(sortIndices(:), :);
    Ytest = ALL_DATA;
    lblstest = zeros(size(ALL_T,1),10);
    for digit = 0:9;
      ind1 = min(find(ALL_T==digit));
      ind2 = max(find(ALL_T==digit));
      lblstest(ind1:ind2,digit+1) = ones(size(ind1:ind2,2),1);
    end

   case 'digits'

    % Fix seeds
    randn('seed', 1e5);
    rand('seed', 1e5);

    load([baseDir 'usps_train.mat']);
    % Extract 600 of digits 0 to 4
    [ALL_T, sortIndices] = sort(ALL_T);
    ALL_DATA = ALL_DATA(sortIndices(:), :);
    Y = [];
    lbls = [];
    numEachDigit = 600;
    for digit = 0:4;
      firstDigit = min(find(ALL_T==digit));
      Y = [Y; ALL_DATA(firstDigit:firstDigit+numEachDigit-1, :)];
      lbl = zeros(1, 5);
      lbl(digit+1) = 1;
      lbls = [lbls; repmat(lbl, numEachDigit, 1)];
    end

   case 'twos'
    % load data
    load([baseDir 'twos']);
    Y = 2*a-1;

   case 'oil'
    load([baseDir '3Class.mat']);
    Y = DataTrn;
    lbls = DataTrnLbls;

   case 'oilTest'
    load([baseDir '3Class.mat']);
    Y = DataTst;
    lbls = DataTstLbls;

   case 'oilValid'
    load([baseDir '3Class.mat']);
    Y = DataVdn;
    lbls = DataVdnLbls;

   case 'oil100'
    randn('seed', 1e5);
    rand('seed', 1e5);
    load([baseDir '3Class.mat']);
    Y = DataTrn;
    lbls = DataTrnLbls;
    indices = randperm(size(Y, 1));
    indices = indices(1:100);
    Y = Y(indices, :);
    lbls = lbls(indices, :);

   case 'swissRoll'
    load([baseDir 'swiss_roll_data']);
    Y = X_data(:, 1:1000)';

   case 'swissRollFull'
    load([baseDir 'swiss_roll_data']);
    Y = X_data(:, 1:3000)';

    %/~
   case 'fourWalks'
    try
      load([baseDir 'fourWalks.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile');

        subj = {'35', '10', '12', '16'};
        examples = {'02', '04', '01', '15'};
        exlbls = eye(4);
        startVals = [55, 222, 22, 62];
        endVals = [338, 499, 328, 342];
        diffInd = [1 3 5];
        totLength = 0;
        totTestLength = 0;
        for i = 1:length(examples)
          skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep subj{i} dirSep subj{i} ...
                              '.asf']);
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep subj{i} dirSep subj{i} '_' ...
                              examples{i} '.amc'], skel);
          tY{i} = tmpchan(startVals(i):4:endVals(i), :);
          % Make given indices differences.
          tY{i}(:, diffInd) = tmpchan(startVals(i)+1:4:endVals(i)+1, diffInd)...
              - tY{i}(:, diffInd);
          tlbls{i} = repmat(exlbls(i, :), size(tY{i}, 1), 1);
          totLength = totLength + size(tY{i}, 1);
        end
        Y = zeros(totLength, size(tY{1}, 2));
        lbls = zeros(totLength, size(tlbls{1}, 2));
        endInd = 0;
        for i = 1:length(tY)
          startInd = endInd + 1;
          endInd = endInd + size(tY{i}, 1);
          Y(startInd:endInd, :) = tY{i};
          lbls(startInd:endInd, :) = tlbls{i};
        end
        Ytest = [];
        lblstest = [];
        subj = {'35', '12', '16', '12', '07', '07', '08', '08'};
        examples = {'03', '02', '21', '03', '01', '02', '01', '02'};
        testexlbls = eye(8);
        startVals = ones(1, 8);
        endVals = 200*ones(1, 8);
        totLength = 0;
        totTestLength = 0;
        for i = 1:length(examples)
          skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep subj{i} dirSep subj{i} ...
                              '.asf']);
          [tmpchan, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep subj{i} dirSep subj{i} '_' ...
                              examples{i} '.amc'], skel);
          tYtest{i} = tmpchan(startVals(i):4:endVals(i), :);
          tYtest{i}(:, diffInd) = tmpchan(startVals(i)+1:4:endVals(i)+1, diffInd)...
              - tYtest{i}(:, diffInd);
          tlblstest{i} = repmat(testexlbls(i, :), size(tYtest{i}, 1), 1);
          totLength = totLength + size(tYtest{i}, 1);
        end
        Ytest = zeros(totTestLength, size(tYtest{1}, 2));
        lblstest = zeros(totTestLength, size(tlblstest{1}, 2));
        endInd = 0;
        for i = 1:length(tYtest)
          startInd = endInd + 1;
          endInd = endInd + size(tYtest{i}, 1);
          Ytest(startInd:endInd, :) = tYtest{i};
          lblstest(startInd:endInd, :) = tlblstest{i};
        end
        save([baseDir 'fourWalks.mat'], 'Y', 'lbls', 'Ytest', 'lblstest');
      else
        error(lasterr);
      end
    end

   case 'walkSitJog'
    Y = mocapLoadTextData([baseDir 'walkSitJog']);
    Y = Y(1:4:end, :)*200;

   case 'robotAtrium'
    Y = parseWirelessData([baseDir 'uw-atrium.txt']);
    Y = (Y + 85)/15;
    Y = Y(1:281, :);
    lbls = 'connect'

    Ytest = parseWirelessData([baseDir 'uw-atrium.txt']);
    Ytest = (Ytest + 85)/15;
    Ytest = Ytest(282:end, :);

   case 'accWalkSitJog'
    try
      load([baseDir 'walkSitJog.mat']);
    catch
      [void, errid] = lasterr;
      if strcmp(errid, 'MATLAB:load:couldNotReadFile')
        skel = acclaimReadSkel([baseDir 'mocap' dirSep 'cmu' dirSep '86' dirSep '86.asf']);
        [channels, skel] = acclaimLoadChannels([baseDir 'mocap' dirSep 'cmu' dirSep '86' dirSep '86_10.amc'], skel);
        % Remove root z & x pos channels.
        channels(:, [1 3]) = zeros(size(channels, 1), 2);
        Y = channels(1:4:end, :);
        save([baseDir 'walkSitJog.mat'], 'Y')
      else
        error(lasterr);
      end
    end
   case 'horse'
    load([baseDir 'horse.dat']);
    horse([133, 309], :) = [];
    ordinalAtr = [1 2 7 8 9 10 11 12 13 14 15 17 18 21];
    horse(:, ordinalAtr) = horse(:, ordinalAtr) - 1;

    % Set labels as outcome --- 1 lived, 2 died, 3 was euthanized
    lbls = zeros(size(horse, 1), 3);
    for i = 1:size(horse, 1)
      lbls(i, horse(i, 23)) = 1;
    end

    horse(:, 23:28) = [];
    horse(:, 3) = [];
    Y = horse;
    % Normalise Gaussian variables.
    for i = [3 4 5 15 18 19 21];
      va = var(Y(find(~isnan(Y(:, i))), i));
      Y(:, i) = Y(:, i)/va;
    end

   case 'pitch'
    load('-ascii', [baseDir 'bballpitch.mat']);
    Y = bballpitch(:,8:end);

   case 'ratemaps'
    % Fix seeds
    Y = loadRateMap([baseDir 'ratemaps'])';
    Y = Y(1:10000, :);

   case 'cepstral'
    Y = load([baseDir 'cepvecs']);
    Y = Y(1:5000, :);

   case 'vowels90'
    load([baseDir 'jon_vowel_data']);
    Y = [a_raw(1:30:end, :); ae_raw(1:30:end, :); ao_raw(1:30:end, :); ...
         e_raw(1:30:end, :); i_raw(1:30:end, :); ibar_raw(1:30:end, :); ...
         o_raw(1:30:end, :); schwa_raw(1:30:end, :); u_raw(1:30:end, :)];
    lbls = [];
    for i = 1:9
      lbl = zeros(1, 9);
      lbl(i) = 1;
      lbls = [lbls; repmat(lbl, 10, 1)];
    end
   case 'vowels3'
    load([baseDir 'jon_vowel_data']);
    Y = [a_raw(1:3:end, :); ae_raw(1:3:end, :); ao_raw(1:3:end, :); ...
         e_raw(1:3:end, :); i_raw(1:3:end, :); ibar_raw(1:3:end, :); ...
         o_raw(1:3:end, :); schwa_raw(1:3:end, :); u_raw(1:3:end, : ...
                                                      )];
    lbls = [];
    for i = 1:9
      lbl = zeros(1, 9);
      lbl(i) = 1;
      lbls = [lbls; repmat(lbl, 100, 1)];
    end
   case 'breakdance'
    Y = mocapLoadTextData([baseDir 'EricCamper04'], 0);
    Y = Y(1:4:end, :);

   case 'dance'
    Y = mocapLoadTextData([baseDir 'dancemotion'], 0);
    Y = Y(1:4:end, :);

   case 'danceReduced'
    Y = mocapLoadTextData([baseDir 'dancemotion'], 0);
    Y = Y(1:16:end, :);
    %~/

   otherwise
    error('Unknown data set requested.')

  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% OP_CURLU_CURLV_3D: assemble the matrix A = [a(i,j)], a(i,j) = (coeff curl u_j, curl v_i), with vector-valued curl.
%
%   mat = op_curlu_curlv_3d (spu, spv, msh, epsilon);
%   [rows, cols, values] = op_curlu_curlv_3d (spu, spv, msh, epsilon);
%
% INPUT:
%
%   spu:   structure representing the space of trial functions (see sp_vector/sp_evaluate_col)
%   spv:   structure representing the space of test functions  (see sp_vector/sp_evaluate_col)
%   msh:   structure containing the domain partition and the quadrature rule (see msh_cartesian/msh_evaluate_col)
%   coeff: physical parameter
%
% OUTPUT:
%
%   mat:    assembled matrix
%   rows:   row indices of the nonzero entries
%   cols:   column indices of the nonzero entries
%   values: values of the nonzero entries
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011, 2017 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function varargout = op_curlu_curlv_3d (spu, spv, msh, coeff)

  rows = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  cols = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);
  values = zeros (msh.nel * spu.nsh_max * spv.nsh_max, 1);

  jacdet_weights = msh.jacdet .* msh.quad_weights .* coeff;

  ncounter = 0;
  for iel = 1:msh.nel
    if (all (msh.jacdet(:,iel)))
      curlu_iel = reshape (spu.shape_function_curls(:, :, :, iel), 3, msh.nqn, 1, spu.nsh_max);
      curlv_iel = reshape (spv.shape_function_curls(:, :, :, iel), 3, msh.nqn, spv.nsh_max, 1);

      jacdet_iel = reshape (jacdet_weights(:,iel), [1,msh.nqn,1,1]);

      jacdet_curlu = bsxfun (@times, jacdet_iel, curlu_iel);
      tmp1 = sum (bsxfun (@times, jacdet_curlu, curlv_iel), 1);
      elementary_values = reshape (sum (tmp1, 2), spv.nsh_max, spu.nsh_max);

      [rows_loc, cols_loc] = ndgrid (spv.connectivity(:,iel), spu.connectivity(:,iel));
      indices = rows_loc & cols_loc;
      rows(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = rows_loc(indices);
      cols(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = cols_loc(indices);
      values(ncounter+(1:spu.nsh(iel)*spv.nsh(iel))) = elementary_values(indices);
      ncounter = ncounter + spu.nsh(iel)*spv.nsh(iel);
    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_curlu_curlv_3d: singular map in element number %d', iel)
    end
  end

  if (nargout == 1)
    varargout{1} = sparse (rows(1:ncounter), cols(1:ncounter), ...
                           values(1:ncounter), spv.ndof, spu.ndof);
  elseif (nargout == 3 || nargout == 0)
    varargout{1} = rows(1:ncounter);
    varargout{2} = cols(1:ncounter);
    varargout{3} = values(1:ncounter);
  else
    error ('op_curlu_curlv_3d: wrong number of output arguments')
  end

end

%% COPY OF THE FIRST VERSION OF THE FUNCTION (MORE UNDERSTANDABLE)
%
% function mat = op_curlu_curlv_3d (spu, spv, msh, coeff)
%
%   mat = spalloc(spv.ndof, spu.ndof, 1);
%   for iel = 1:msh.nel
%     if (all (msh.jacdet(:,iel)))
%       mat_loc = zeros (spv.nsh(iel), spu.nsh(iel));
%       for idof = 1:spv.nsh(iel)
%         ishc = reshape(spv.shape_function_curls(:,:,idof,iel), 3, []);
%         for jdof = 1:spu.nsh(iel)
%           jshc = reshape (spu.shape_function_curls(:,:,jdof,iel), 3, []);
% % The cycle on the quadrature points is vectorized
%           %for inode = 1:msh.nqn
%             mat_loc(idof, jdof) = mat_loc(idof, jdof) + ...
%               sum (msh.jacdet(:,iel) .* msh.quad_weights(:, iel) .*...
%               sum (ishc .* jshc, 1).' .* coeff (:, iel));
%           %end
%         end
%       end
%       mat(spv.connectivity(:,iel), spu.connectivity(:,iel)) = ...
%         mat(spv.connectivity(:,iel), spu.connectivity(:,iel)) + mat_loc;
%     else
%       warning ('geopdes:jacdet_zero_at_quad_node', 'op_curlu_curlv_3d: singular map in element number %d', iel)
%     end
%   end
%
% end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [state] = tapas_sampler_mc3g(data, model, inference, state)
%% Sampler for population mcmc (mc3) with generalized temperature.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

ollh = state.llh{1};
olpp = state.llh{2};

% Get the joint
%for i = 2:nn - 1
%    lpp = lpp + sum(state.llh{i}, 1);
%end

T = model.graph{1}.htheta.T;
[nv, nc] = size(T);
v = reshape([1:nv * nc], nv, nc);

% Doesn't need to do anything
if nc == 1
    return
end

c = 0;
for l = 1 : inference.mc3it
    s1 = ceil(rand() * nv);
    s2 = ceil(rand() * (nc - 1));
    % Realinate
    [s1, s2i] = ind2sub([nv, nc], v(s1, s2));
    [s1, s2j] = ind2sub([nv, nc], v(s1, s2 + 1));

    % Now we need to adjust the probability
    plpp1 = model.graph{2}.llh_sn(state.graph{2}, ...
        state.graph{3}, model.graph{2}.htheta, s1, s2, s2j);
    plpp2 = model.graph{2}.llh_sn(state.graph{2}, ...
        state.graph{3}, model.graph{2}.htheta, s1, s2 + 1, s2i);

    % Compute the ratio
    p = exp((ollh(s1, s2i) .* T(s1, s2 + 1) + ...
        ollh(s1, s2j) .* T(s1, s2) + ...
        plpp1 + plpp2) - ...
        (ollh(s1, s2i) .* T(s1, s2) + ...
        ollh(s1, s2j) .* T(s1, s2 + 1) + ...
        olpp(s1, s2) + olpp(s1, s2 + 1) ));
    if rand() < p
        c = c + 1;
        v(s1, [s2, s2 + 1]) = v(s1, [s2 + 1, s2]);
        olpp(s1, s2) = plpp1;
        olpp(s1, s2 + 1) = plpp2;
    end
end
%:fprintf(1, 'Swap rate %0.5f\n', c/inference.mc3it);
state.llh{1} = state.llh{1}(v);
state.llh{2} = olpp;
state.graph{2}.y = state.graph{2}.y(v);

end

author: translationalneuromodeling

repo: tapas


text:
%************************************************************************
%                       Bezier_circle                                   *
%                                                                       *
%           This program for practising CAO - EMMC11                    *
%           Professor :     Pierre BECKERS                              *
%           (Aerospace Laboratory, University of Liege, Belgium)        *
%           Student   :     NGUYEN QUOC DUAN                            *
%                          ( ngqduan@yahoo.com )                        *
%                                                                       *
%************************************************************************

%*** Problem ************************************************************
% Given 6 points (pi = 3.1416) :                                        *
%   Pi = [rcos(-pi/3+i*pi/3)    rsin(-pi/3+i*pi/3)] for  i = 0,1,...5   *
% Defined the Bezier curve starting in P0, ending in P5 and passing     *
% through points P1, P2, P3,P4. Explain the method and show all         *
% hypotheses.                                                           *
% ***********************************************************************

clear all; close all; clc;

% INPUT DATA

r=1;            % the length of the radius
pi=4*atan(1);   % the value of pi
P=[];           % coordinates matrix of 6 initial points
for i=1:6
    P(i,:)=[r*cos(-pi/3+(i-1)*pi/3) r*sin(-pi/3+(i-1)*pi/3)];
end

% COMPUTE THE COEFFICIENTS MATRICES

n=5;                % the polynomial degree of the Bezier curve
coefx=zeros(3,3);
coefy=zeros(3,3);
t=[0:0.2:1];        % matrix of the values of parameter at 6 given points
for k=1:3
    for j=1:3
        i=j-1;
        coefx(k,j)=C(i,n)*(B_para(i,n,t(k))-B_para(n-i,n,t(k)));
        coefy(k,j)=C(i,n)*(B_para(i,n,t(k))+B_para(n-i,n,t(k)));
    end
end

% COMPUTE THE COORDINATES OF BEZIER POINTS

% Initial values
xB=zeros(6,1);          % x-coordinates of 6 Bezier points
yB=zeros(6,1);          % y-coordinates of 6 Bezier points
% Solve the two equations systems
    xB(1:3)=coefx\P(1:3,1);  % x-coordinates of the first 3 Bezier points
    yB(1:3)=coefy\P(1:3,2);  % y-coordinates of the first 3 Bezier points
% Apply the symmetry of 6 given points
for i=4:6
    xB(i)=-xB(6-i+1);
    yB(i)=yB(6-i+1);
end

% OUTPUT THE SOLUTIONS
number=[0:5]';
disp('6 Bezier points coordinates : ' )
[number xB yB]

% PLOT FIGURE
figure(1);
h=gcf;
set(h,'name','Bezier control points & Bezier curve');
set(h,'NumberTitle','off');
axis equal;
title('Bezier control points & Bezier curve');

% Conecting 6 Bezier points by straight segments
plot(xB,yB,'r-')
hold on;

% Locate 4 Bezier control points
text(xB(2),yB(2),'Q1');
text(xB(3),yB(3),'Q2');
text(xB(4),yB(4),'Q3');
text(xB(5),yB(5),'Q4');
hold on;
% plot Bezier curve
precision=0.01;
for k=0:1/precision
    tp=k*precision;
    xp(k+1)=0;
    yp(k+1)=0;
    for j=1:6
        i=j-1;
        xpj=xB(j)*C(i,n)*B_para(i,n,tp);
        ypj=yB(j)*C(i,n)*B_para(i,n,tp);
        xp(k+1)=xpj+xp(k+1);
        yp(k+1)=ypj+yp(k+1);
    end
end
plot(xp,yp,'g-')
hold on;

% Locate 6 given points
text(P(1,1),P(1,2),'Po');
text(P(2,1),P(2,2),'P1');
text(P(3,1),P(3,2),'P2');
text(P(4,1),P(4,2),'P3');
text(P(5,1),P(5,2),'P4');
text(P(6,1),P(6,2),'P5');

author: Sable

repo: mcbench-benchmarks


text:
function [U,D] = unique_images(dir_name)
  % UNIQUE_IMAGES Generate a list of unique and duplicate images in a
  % directory.
  %
  % [U,D] = unique_images(dir_name)
  %
  % Inputs:
  %   dir_name  path to directory
  % Outputs:
  %   U list of unique images/image paths
  %   D list of duplicate iamges, so that [U,D] are all images
  %

  assert(isdir(dir_name));
  files = imdir(dir_name);
  U = {};
  ucount = 0;
  D = {};
  dcount = 0;
  valid = true(size(files));

  for ii = 1:numel(files)
    fii = [dir_name '/' files(ii).name];
    imii = imread(fii);
    for jj = (ii+1):numel(files)
      imjj = imread([dir_name '/' files(jj).name]);
      if all(size(imii) == size(imjj)) && all(imii(:) == imjj(:))
        valid(jj) = false;
      end
    end
  end
  U = files(valid);
  D = files(~valid);

end

author: alecjacobson

repo: gptoolbox


text:
% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% A quadrotor is a simple flying vehicle with 4 propellers with their thrust
% axes pointing upwards.  This demo shows an animation of a flying quadrotor
% with a nested control system.

sl_quadrotor

% Running the simulation shows the vehicle takes off and flying in a cirle.
sim('sl_quadrotor');

author: Allopart

repo: rbpf-gmapping


text:
function [confmatrix] = cfmatrix2(actual, predict, classlist, per, printout)
% CFMATRIX2 calculates the confusion matrix for any prediction
% algorithm ( prediction algorithm generates a list of classes to which
% each test feature vector is assigned );
%
% Outputs: confusion matrix
%
%                 Actual Classes
%                   p       n
%              ___|_____|______|
%    Predicted  p'|     |      |
%      Classes  n'|     |      |
%
%           Also the TP, FP, FN and TN are output for each class based on
%           http://en.wikipedia.org/wiki/Confusion_matrix
%           The Precision, Sensitivity and Specificity for each class have
%           also been added in this update along with the overall accuracy
%           of the model ( ModelAccuracy ).
%
%
% Further description of the outputs:
%
% True Postive [TP] = Condition Present + Positive result
% False Positive [FP] = Condition absent + Positive result [Type I error]
% False (invalid) Negative [FN] = Condition present + Negative result [Type II error]
% True (accurate) Negative [TN] = Condition absent + Negative result
% Precision(class) = TP(class) / ( TP(class) + FP(class) )
% Sensitivity(class) = Recall(class) = TruePositiveRate(class)
% = TP(class) / ( TP(class) + FN(class) )
% Specificity ( mostly used in 2 class problems )=
% TrueNegativeRate(class)
% = TN(class) / ( TN(class) + FP(class) )
%
% Inputs:
%
% 1. actual / 2. predict
% The inputs provided are the 'actual' classes vector
% and the 'predict'ed classes vector. The actual classes are the classes
% to which the input feature vectors belong. The predicted classes are the
% class to which the input feature vectors are predicted to belong to,
% based on a prediction algorithm.
% The length of actual class vector and the predicted class vector need to
% be the same. If they are not the same, an error message is displayed.
% 3. classlist
% The third input provides the list of all the classes {p,n,...} for which
% the classification is being done. All classes are numbers.
% 4. per = 1/0 (default = 0)
% This parameter when set to 1 provides the values in the confusion matrix
% as percentages. The default provides the values in numbers.
% 5. printout = 1/0 ( default = 1 )
% This parameter when set to 1 provides output on the matlab terminal and
% can be used to suppress output by setting to 0. ( default = 1 ). Assuming
% 'printout' of output use case would be more common and at the same time
% provided option to suppress output when the number of classes can be very
% large.
%
% Example:
% >> a = [ 1 2 3 1 2 3 1 1 2 3 2 1 1 2 3];
% >> b = [ 1 2 3 1 2 3 1 1 1 2 2 1 2 1 3];
% >> Cf = cfmatrix2(a, b, [1 2 3], 0, 1);
% is equivalent to
% >> Cf = cfmatrix2(a, b);
% The values of classlist(unique from actual), per(0), printout(1) are set
% to the respective defaults.
%
%
% [Avinash Uppuluri: avinash_uv@yahoo.com: Last modified: 03/28/2012]
%
% Changes added for 03/28/2012 upload
% a. Pre-initialize confmatrix
% b. Simplified logic making the code more readable and faster;
%    (based on comments from an interviewer who reviewed the code)
% c. Provide input variable 'printout' as an option to suppress output to
%    screen ( output to display is still the default (printout = 1)
%    assuming that will be the more common use case ).
% d. Added Precision(class), Sensitivity(class), Specificity(class) and
%    the overall accuracy of model calculations.

% If classlist not entered: make classlist equal to all
% unique elements of actual
if (nargin < 2)
   error('Not enough input arguments. Need atleast two vectors as input');
elseif (nargin == 2)
    classlist = unique(actual); % default values from actual
    per = 0;
    printout = 1;
elseif (nargin == 3)
    per = 0; % default is numbers and input 1 or higher for percentage
    printout = 1;
elseif (nargin == 4)
    printout = 1; % default is silent output ( 0 ); one or higher printsout
elseif (nargin > 5)
   error('Too many input arguments.');
end

if (length(actual) ~= length(predict))
    error('First two inputs need to be vectors with equal size.');
elseif ((size(actual,1) ~= 1) && (size(actual,2) ~= 1))
    error('First input needs to be a vector and not a matrix');
elseif ((size(predict,1) ~= 1) && (size(predict,2) ~= 1))
    error('Second input needs to be a vector and not a matrix');
end
format short g;
n_class = length(classlist);
confmatrix = zeros(n_class);
line_two = '----------';
line_three = '_________|';

for i = 1:n_class
    for j = 1:n_class
        m = (predict == classlist(i) ...
           & actual  == classlist(j));
        confmatrix(i,j) = sum(m);
    end
    line_two = strcat(line_two,'---',num2str(classlist(i)),'-----');
    line_three = strcat(line_three,'__________');
end

% True Postive [TP] = Condition Present + Positive result
% False Positive [FP] = Condition absent + Positive result [Type I error]
% False (invalid) Negative [FN] = Condition present + Negative result [Type II error]
% True (accurate) Negative [TN] = Condition absent + Negative result
% Precision(class) = TP(class) / ( TP(class) + FP(class) )
% Sensitivity(class) = Recall(class) = TruePositiveRate(class)
% = TP(class) / ( TP(class) + FN(class) )
% Specificity ( mostly used in 2 class problems )=
% TrueNegativeRate(class)
% = TN(class) / ( TN(class) + FP(class) )

TPFPFNTN    = zeros(4, n_class);
Precision   = zeros(1, n_class);
Sensitivity = zeros(1, n_class);
Specificity = zeros(1, n_class);

temps1 = sprintf('    TP  ');
temps2 = sprintf('    FP  ');
temps3 = sprintf('    FN  ');
temps4 = sprintf('    TN  ');
temps5 = sprintf('Preci.  ');
temps6 = sprintf('Sensi.  ');
temps7 = sprintf('Speci.  ');

for i = 1:n_class
    TPFPFNTN(1, i) = confmatrix(i,i); % TP
    temps1 = strcat(temps1,sprintf(' |   %2.2f    ',TPFPFNTN(1, i)));
    TPFPFNTN(2, i) = sum(confmatrix(i,:))-confmatrix(i,i); % FP
    temps2 = strcat(temps2,sprintf(' |   %2.2f    ',TPFPFNTN(2, i) ));
    TPFPFNTN(3, i) = sum(confmatrix(:,i))-confmatrix(i,i); % FN
    temps3 = strcat(temps3,sprintf(' |   %2.2f    ',TPFPFNTN(3, i) ));
    TPFPFNTN(4, i) = sum(confmatrix(:)) - sum(confmatrix(i,:)) -...
        sum(confmatrix(:,i)) + confmatrix(i,i); % TN
    temps4 = strcat(temps4,sprintf(' |   %2.2f    ',TPFPFNTN(4, i) ));
    % Precision(class) = TP(class) / ( TP(class) + FP(class) )
    Precision(i)   = TPFPFNTN(1, i) / sum(confmatrix(i,:));
    temps5 = strcat(temps5,sprintf(' |   %1.2f    ',Precision(i) ));
    % Sensitivity(class) = Recall(class) = TruePositiveRate(class)
    % = TP(class) / ( TP(class) + FN(class) )
    Sensitivity(i) = TPFPFNTN(1, i) / sum(confmatrix(:,i));
    temps6 = strcat(temps6,sprintf(' |   %1.2f    ',Sensitivity(i) ));
    % Specificity ( mostly used in 2 class problems )=
    % TrueNegativeRate(class)
    % = TN(class) / ( TN(class) + FP(class) )
    Specificity(i) = TPFPFNTN(4, i) / ( TPFPFNTN(4, i) + TPFPFNTN(2, i) );
    temps7 = strcat(temps7,sprintf(' |   %1.2f    ',Specificity(i) ));
end

ModelAccuracy = sum(diag(confmatrix))/sum(confmatrix(:));
temps8 = sprintf('Model Accuracy is %1.2f ',ModelAccuracy);

if (per > 0) % ( if > 0 implies true; < 0 implies false )
    confmatrix = (confmatrix ./ length(actual)).*100;
end

if ( printout > 0 ) % ( if > 0 printout; < 0 no printout )
    disp('------------------------------------------');
    disp('             Actual Classes');
    disp(line_two);
    disp('Predicted|                     ');
    disp('  Classes|                     ');
    disp(line_three);

    for i = 1:n_class
        temps = sprintf('       %d             ',i);
        for j = 1:n_class
            temps = strcat(temps,sprintf(' |    %2.1f    ',confmatrix(i,j)));
        end
        disp(temps);
        clear temps
    end
    disp('------------------------------------------');

    disp('------------------------------------------');
    disp('             Actual Classes');
    disp(line_two);
    disp(temps1); disp(temps2); disp(temps3); disp(temps4);
    disp(temps5); disp(temps6); disp(temps7);
    disp('------------------------------------------');
    disp(temps8);
    disp('------------------------------------------');
end
clear temps1 temps2 temps3 temps4 temps5 temps6 temps7 temps8


text:
function hhof = plotHofFeatures(H,hax,maxv,colors,binedges,psize,offset)

if nargin<7,
  offset = [0 0];
end

hhof = [];
for xi = 1:size(H,2),
  cx = (psize/2 + 1 + (xi-1)*psize);
  for yi = 1:size(H,1),
    cy = (psize/2 + 1 + (yi-1)*psize);

    for bini = 1:size(H,3),
      tmp = linspace(binedges(bini),binedges(bini+1),20);
      xcurr = cx + [0,psize/2*cos(tmp),0];
      ycurr = cy + [0,psize/2*sin(tmp),0];
      hhof(yi,xi,bini) = patch(xcurr+offset(1),ycurr+offset(2),colors(bini,:),...
        'LineStyle','none','FaceAlpha',min(1,H(yi,xi,bini)/maxv),'Parent',hax);
    end

  end
end

author: kristinbranson

repo: JAABA


text:
function writeSmoothAbs()
%
% FUNCTION:
%   This function is used to write the code for smoothAbs
%

MaxOrder = 8;  %Should be at least 1
N = MaxOrder + 1;

%% Set up the constraints:

% The value and slope should be 1 at the right transition point, and zero
% for any higher derivatives.
xRight = ones(1,N);
yRight = zeros(1,N);
yRight(1) = 1; yRight(2) = 1;

% The value and slope should be 1 at the right transition point, and zero
% for any higher derivatives.
xLeft = -ones(1,N);
yLeft = zeros(1,N);
yLeft(1) = 1; yLeft(2) = -1;

% All Constraints are expressed across all derivatives
dLevel = 0:N;

%% Begin writing the file:

fid = fopen('smoothAbs.m','w');

fprintf(fid,'function y = smoothAbs(x,alpha,order)\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% y = smoothAbs(x,alpha,order)\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% FUNCTION:\n');
fprintf(fid,'%%   This function is a smooth version of abs(). Smoothing is done using\n');
fprintf(fid,'%%   a piecewise polynomial, which is Nth order continuous at the\n');
fprintf(fid,'%%   transitions, where N is a value between 0 and 8, and input as ''order''\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% INPUTS:\n');
fprintf(fid,'%%   x = a vector of matrix of real numbers\n');
fprintf(fid,'%%   alpha = smoothing parameter: -alpha<x<alpha will be smoothed\n');
fprintf(fid,'%%   order = the order of the smoothing (0=no smoothing)\n');
fprintf(fid,['%%           order must be in the set {0,1,..,' num2str(MaxOrder) '}\n']);
fprintf(fid,'%%\n');
fprintf(fid,'%% OUTPUTS:\n');
fprintf(fid,'%%   y = a smoothed version of abs(x)\n%%\n\n');
fprintf(fid,'switch order \n');
fprintf(fid,'    case 0 %% No Smoothing \n');
fprintf(fid,'        y=abs(x); \n');

for i=2:N
    x = [xLeft(1:i),xRight(1:i)];
    y = [yLeft(1:i),yRight(1:i)];
    d = [dLevel(1:i),dLevel(1:i)];
    [~,pSym] = SolveCoeff(x,y,d);
    str = 'p = [';
    for j=2:length(pSym)  %Leading coefficient is always zero
        if j==length(pSym)
            str = [str char(pSym(j)) ']; \n']; %#ok<*AGROW>
        else
            str = [str char(pSym(j)) ','];
        end
    end
    fprintf(fid,['    case ' num2str(i-1) ' \n']);
    fprintf(fid,'        y=x; \n');
    fprintf(fid,'        c1 = x < -alpha; \n');
    fprintf(fid,'        c2 = x > alpha; \n');
    fprintf(fid,'        y(c1)=-x(c1); \n');
    fprintf(fid,['        ' str]);
    fprintf(fid,'        idx = ~c1&~c2;\n');
    fprintf(fid,'        y(idx)=alpha*polyval(p,x(idx)/alpha); \n');
end

fprintf(fid,'    otherwise \n');
fprintf(fid,'        error(''Order not supported'');\n');
fprintf(fid,'end \n');

fprintf(fid,'end \n');

fclose(fid);
end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [beta_draws,sigma_draws,IV_draws,C_draws,D,gamma,storage1,storage2]=irfIV6(betahat,n,m,p,k,T,names,startdate,enddate,X,Y,endo,data_endo,data_exo,const,pref,strctident,IRFt,IRFperiods,It,Bu)

% inputs:  - matrix 'betahat': OLS estimate for beta
%          - matrix 'sigmahat': OLS estimate for sigma
%          - integer 'IRFperiods': number of periods for IRFs
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
%          - cell 'signrestable': table recording the sign restriction input from the user
%          - cell 'signresperiods': table containing the periods corresponding to each restriction
%          - string 'ShockwithInstrument' Name of the shock where the instrument belongs to
% outputs: - cell 'struct_irf_record': record of the draws for the orthogonalised IRFs
%          - matrix 'D_record': record of the draws for the structural matrix D
%          - matrix 'gamma_record': record of the gibbs sampler draws for the structural disturbances variance-covariance matrix gamma
%          - matrix 'hd_record': record of historical decompositions
%          - matrix 'ETA_record': record of structural Shocks
%          - vector 'IVcorrelation': vector of correlation with proxy proxy
%          - vector 'PofIVcorrelation': corresponding P-value

%% Phase 1:% %% IV F-Test
% Load IV and make it comparable with the reduced form errors
[EPSIV,IVcut,EPSt,sigmahatIV,B,EPS,df,sigma_hat,inv_sigma_hat,beginInstrument,EndInstrument,IV,txt,OverlapIVinY,cut1,cut2,cut3,cut4]=bear.loadIV(betahat,k,n,Y,X,T,p,names,startdate,enddate,strctident,pref);

%%  Check strength of instrument
F_stat_test_unc   = NaN*ones(1,n);
signif_test_unc   = NaN*ones(1,n);
Betas_test_unc    = NaN*ones(1,n);
Rsquared_test_unc = NaN*ones(1,n);

for i = 1:n
    yyy      = EPSt(i,:)'; %%get residuals
    xxx_rest = ones(size(EPSt,2),1); %%vector of constants
    xxx      = [ones(size(EPSt,2),1), IVcut]; %%constants plus instrument

    % restricted model
    coeff_rest = (xxx_rest'*xxx_rest)^(-1)*xxx_rest'*yyy;
    resid_rest = yyy - xxx_rest*coeff_rest;
    RSS_rest   = resid_rest'*resid_rest;

    % unrestricted model
    coeff = (xxx'*xxx)^(-1)*xxx'*yyy;
    resid = yyy - xxx*coeff;
    RSS   = resid'*resid;

    % compute statistic
    F_stat   = (RSS_rest-RSS)*(size(yyy,1)-size(xxx,2))/(RSS*(size(xxx,2)-size(xxx_rest,2)));
    F_pvalue = 1-fcdf(F_stat,size(xxx,2)-size(xxx_rest,2),size(yyy,1)-size(xxx,2));

    F_stat_test_unc(i) = F_stat;
    significance = 3*(F_pvalue < 0.01);
    significance = significance + 2*(F_pvalue < 0.05 & F_pvalue > 0.01);
    significance = significance + 1*(F_pvalue < 0.1 & F_pvalue > 0.05);
    signif_test_unc(i) = significance;

    % regression statistics
    coeff  = (xxx'*xxx)^(-1)*xxx'*yyy;
    fitted = xxx*coeff;
    resid  = yyy - fitted;
    yyy_demean           = yyy - mean(yyy);
    Rsquared_test_unc(i) = 1 - resid'*resid/(yyy_demean'*yyy_demean);

    Betas_test_unc(i) = coeff(2);
end

%% Print F Statistic

table_validitytest_rows = char('coefficient', 'significance level', 'number of obs', 'F statistic', 'R squared');
table_validitytest_rows = cellstr(table_validitytest_rows);
table_validitytest_values  = [Betas_test_unc;  signif_test_unc;  length(IVcut)*ones(1,n); F_stat_test_unc; Rsquared_test_unc ];
table_validitytest_columns = cellstr(endo)';

table = array2table(table_validitytest_values);
for kk=1:length(table_validitytest_columns)
    table.Properties.VariableNames(kk) = table_validitytest_columns(1,kk);
end
table_rows = cell2table(table_validitytest_rows);
table=[table_rows,table];

for kk=1
    table.Properties.VariableNames(kk) = cellstr(char('Validitytest'));
end

% print proxy var statistics in the command window and in the results file
filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'at');
%two empty lines
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
%print Proxy VAR statistics title
fprintf('%s\n','Proxy VAR statistics');
fprintf(fid,'%s\n','Proxy VAR statistics');
fprintf('%25s %15s %15s %15s %15s\n','','coefficient','significance','number of obs','F statistic');
fprintf(fid,'%25s %15s %15s %15s %15s\n','','coefficient','significance','number of obs','F statistic');
% handle the endogenous
   for jj=1:n
      values=[Betas_test_unc(1,jj);  signif_test_unc(1,jj);  length(IVcut); F_stat_test_unc(1,jj)];
      fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1}),values);
      fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1}),values);
   end
%empty line
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fclose(fid);

%% first stage regression

%step 2: Regress the first reduced form shock on the instrument
Shock = EPSIV(:,1);
[nobs,~] = size(IVcut);
XX = [ones(nobs,1) IVcut];
[~,nvar] = size(XX);
%get OLS estimate
XpXi = (XX'*XX)\eye(nvar);
betaIV=XpXi*(XX'*Shock);
%get predicted value
IVpred = XX*betaIV;

%% second stage regression
%step 3: Regress the other reduced form shocks on the predicted value
ImpactIRFIV = zeros(n,1);
ImpactIRFIV(1,1) = 1;

for hh=2:n
Shock = EPSIV(:,hh);
[nobs,~] = size(IVpred);
IVpredtemp = [ones(nobs,1) IVpred];
[~,nvar] = size(IVpredtemp);
IVpIVi = (IVpredtemp'*IVpredtemp)\eye(nvar);
betaIV2=IVpIVi*(IVpredtemp'*Shock);
ImpactIRFIV(hh,1) = betaIV2(2,1); %should be equal to Mu from 2:end
end

% step 5: Create the structural matrix and only fill the first column as
% this is the only one identified
%Dols=eye(n,n);
%%another way to retrieve b11 (the scalar that scales the IRF to be a 1sdt Shock) is simply
C=chol(bear.nspd(sigmahatIV),'lower');

if IRFt==5 % D & gamma output
    % step 5: Create the structural matrix and only fill the first column as
% this is the only one identified
D=eye(n,n);
% Step 6: Replace the first Column in the Cholesky Decomposition by
%the structural impact matrix computed above
% D(1:end,1) = ImpactIRFIV*b11;

gamma=eye(n);
%%another way to retrieve b11 (the scalar that scales the IRF to be a 1sdt Shock) is simply
b=ImpactIRFIV;
%%Recover the vector q that maps the first column of C into b such that Cq=b;
q = C\b;
%%b11 is the euclidian length of q
b11q = 1/norm(q);
if strctident.stdeviation ==1
D(1:end,1) = ImpactIRFIV*b11q;
else
D(1:end,1) = ImpactIRFIV;
end
end

%% Preliminiaries for sign restrictions
Acc=It-Bu;  %Minimum accepted draws for the rotation matrix

%% Storage cells and other preliminiaries for the bootstrap part
endo_artificial = zeros(T+p,n); %initialize artifical data matrix
%preliminaries for IV block bootstrap
BlockSize = p; %this can be made frequency specific %only used for bootstraptype2
nBlock = ceil((T)/BlockSize); %only used for bootstraptype2
startofcommonsample= find(OverlapIVinY==1, 1, 'first');
endofcommonsample= find(OverlapIVinY==1, 1, 'last');
%EPShatcut = EPS(OverlapIVinY,:);
beta_draws = nan(k*n,Acc*10);
sigma_draws=nan(n^2,Acc*10);
IV_draws=nan(n,Acc*10);
C_draws=nan(n^2,Acc*10);
beta = betahat;
B=reshape(beta,k,n);

%%%%while AA<=BB with BB=1000;
%%%%%%%%%%%%%%%%%%%%Bootstrap to approximate distribution of reduced form and IV coefficients%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hbar=bear.parfor_progressbar(Acc*10,'IV progress');
for AA=1:Acc*10 %parfor loop possible here? %%%% are 100000 draws necessary? parfor AA=1:Acc*100
    stationary=0;
while stationary==0
        %%Get reduced form coefficients and IV column
        %% STEP 1: generate the artificial data
        if strctident.bootstraptype==1
            % Initialize wild bootstrap by drawing residuals and flipping the sign at
            % random get reduced form residuals
            EPS  = Y-X*B;
            rotationvector = 1-2*(rand(T,1)>0.5); %the rotation vector is of the same size as EPS. We need to determine where EPS starts in terms of IV and where it ends
            EPSrotate = EPS.*(rotationvector*ones(1,n));
            %create IV specific rotation vector that is 1 if the samples dont
            %overlapp and the EPS rotation vector otherwhise
            %         IVrotationvector = rotationvector(OverlapIVinY);
            %         IVrotate = IV.*IVrotationvector;

            IVrotationvector = ones(length(IV),1);
            for kk=startofcommonsample:(endofcommonsample-startofcommonsample)
                if OverlapIVinY(kk,1)==1
                    IVrotationvector(kk,1) = rotationvector(kk,1);
                else
                    IVrotationvector(kk,1) = 1;
                end
            end
            IVrotate = IV.*IVrotationvector;

        else
            %%Moving block bootstrap bootstrap
            EPS  = Y-X*B;
            IVblocks = zeros(BlockSize,size(IV,2),length(IV)-BlockSize+1); %create moving blocks for the instrument
            for yy = 1:length(IV)-BlockSize+1
                IVblocks(:,:,yy) = IV(yy:BlockSize+yy-1,:);
            end

            EPSblocks = zeros(BlockSize,n,T-BlockSize+1); %create moving blocks for the reduced form errors
            for yy = 1:T-BlockSize+1
                EPSblocks(:,:,yy) = EPS(yy:BlockSize+yy-1,:);
            end

            %the first block that can be accomodated by the instrument (i.e. where reduced form errors and instrument overlap
            %is given by EPSblocks(:,:,startofcommonsaple);

            %the last block of reduced form errors that can be accomodated by
            %the instrument is EPSblocks(:,:,endofcommonsample-BlockSize);

            %preliminaries for centering of the bootstrapped VAR errors
            centering = zeros(BlockSize,n);
            for yy = 1:BlockSize
                centering(yy,:) = mean(EPS(yy:T-BlockSize+yy,:,1),1); %running mean
            end

            centering = repmat(centering,[nBlock,1]);
            centering = centering(1:T,:);

            %center the bootstrapped proxy variables
            IVcentering = zeros(BlockSize,size(IV,2));
            for yy = 1:BlockSize %moving average of the instrument
                subIV = IV(yy:size(IV,1)-BlockSize+yy,:);
                IVcentering(yy,:) = mean(subIV((subIV(:,1) ~= 0),1),1); %running mean
            end
            IVcentering = repmat(IVcentering,[nBlock,1]);
            IVcentering = IVcentering(1:size(IV,1),:);

            pickblock = randi([startofcommonsample endofcommonsample-BlockSize+1],nBlock,1);
            pickblockIV = pickblock-startofcommonsample+1;
            %now check if we can replace some blocks in the beginning of the sample
            %with other residuals, that are not accomodated by the instrument
            missing_beginning = startofcommonsample;
            %missing_end = T-endofcommonsample;

            blocks_beginning = floor(missing_beginning/BlockSize); %number of blocks that we can possibly replace in the beginning
            blocks_end = floor((T-endofcommonsample)/BlockSize); %number of blocks that we can possibly replace in the end

            if blocks_beginning > 0 %if the instrument and the sample dont coincide for at least the length of the first block
                pickblock_beginning = randi([1 startofcommonsample],blocks_beginning,1);
                pickblock(1:blocks_beginning,:)=pickblock_beginning(1:blocks_beginning);
            end

            if blocks_end > 0 %if we can resample a block from the blocks that are not accomodated by the instrument do so
                pickblock_end = randi([endofcommonsample-BlockSize+2 T-BlockSize+1],blocks_end,1);
                howmanyblock = length((nBlock-blocks_end) : nBlock) ;
                if howmanyblock > length(pickblock_end)
                pickblock((nBlock-blocks_end)+1 : nBlock)=pickblock_end(1:blocks_end);
                else
                pickblock((nBlock-blocks_end) : nBlock)=pickblock_end(1:blocks_end);
                end
            end

            EPSrotate = zeros(nBlock*BlockSize,n);
            IVrotate = zeros(nBlock*BlockSize,size(IV,2));
            for yy = 1:nBlock
                EPSrotate(1+BlockSize*(yy-1):BlockSize*yy,:) = EPSblocks(:,:,pickblock(yy,1));
                IVrotate(1+BlockSize*(yy-1):BlockSize*yy,:) = IVblocks(:,:,pickblockIV(yy,1));
            end

            EPSrotate = EPSrotate(1:T,:); %drop last pseudo observations
            IVrotate = IVrotate(1:T,:); %drop last pseudo observations
%             IVrotate = IVrotate(1:length(IV),:);
            IVrotate = IVrotate(OverlapIVinY,:); %drop observations where instrument is not available
            %EPScutrotate = EPSrotate(OverlapIVinY,:);
            %center the bootstrapped residuals and proxies by subtracting
            %the mean
            EPSrotate = EPSrotate - centering;
            for yy = 1:size(IVrotate,2)
                IVrotate((IVrotate(:,yy)~=0),yy) = IVrotate((IVrotate(:,yy)~=0),yy) - IVcentering((IVrotate(:,yy)~=0),yy);
            end

        end

        %% STEP 2: Generate artificial data
        [endo_artificial_gen]=bear.gen_fake_data(endo_artificial,data_endo,B,EPSrotate,const,p,n,T);

        %% STEP 3: estimate reduced form VAR on artificial data.
        [~,betadraw,sigmadraw,Xdraw,~,Ydraw]=bear.olsvar(endo_artificial_gen,data_exo,const,p);
        %%%%% do we need this test for irft6 only? or is it also useful for
        %%%%% IRFt5
        if IRFt==5
        stationary=1;
        elseif IRFt==6
        [stationary]=bear.checkstable(betadraw,n,p,k); %switches stationary to 0, if the draw is not stationary
        end
end
        if IRFt==5
            %% STEP 4: identify the model
    [Ddraw]=bear.irfiv_ols_for_bootstrap_GK(txt,beginInstrument,EndInstrument,names,IVrotate,betadraw,sigmadraw,m,n,Xdraw,Ydraw,k,p,enddate,startdate,endo,IRFperiods,pref,cut1,cut2,cut3,cut4,strctident);
    %% Step 5: Calculate impulse responses and store them
    [~,ortirfmatrix]=bear.irfsim(betadraw,Ddraw,n,m,p,k,IRFperiods);
       for jj=1:IRFperiods
       storage1{AA,1}(:,:,jj)=ortirfmatrix(:,:,jj);
       end
       storage2{AA,1}=Ddraw;
%        beta_gibbs(:,AA)=betadraw;
%        sigma_gibbs(:,AA)=bear.vec(sigmadraw);
        elseif IRFt==6
        %% Step 4: Identifity the proxy VAR equation
        [~,~,q,C]=irfiv_ols_for_ivres(txt,names,IVrotate, betadraw,sigmadraw,m,n,Xdraw,Ydraw,k,p,enddate,startdate, endo, cut1, cut2,cut3,cut4);
beta_draws(:,AA)=betadraw;
sigma_draws(:,AA)=bear.vec(sigmadraw);
%D_draws(:,AA)=bear.vec(D);
IV_draws(:,AA)=q/norm(q);
C_draws(:,AA)=bear.vec(C);
%IVrotate_draws(:,AA)=IVrotate;
        end
%%%AA = AA+1 ;
hbar.iterate(1); % update progress by one iteration
end
close(hbar);   %close progress bar

%%%%statistics here?

author: european-central-bank

repo: BEAR-toolbox


text:
function varargout = plotcoeffs(f, varargin)
%PLOTCOEFFS   Display coefficients graphically.
%   PLOTCOEFFS(F) calls PLOTCOEFFS(F.SMOOTHPART).
%
%   This method is mainly aimed at Chebfun developers.
%
% See also PLOT.

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

% Deal with an empty input:
if ( isempty(f) )
    if ( nargout == 1 )
        [varargout{1:nargout}] = plot([]);
    end
    return
end

% Call PLOTCOEFFS() on the SMOOTHPART of F.
[varargout{1:nargout}] = plotcoeffs(f.smoothPart, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% input: uv_c1 N x 2 point array in c1 frame,
% output: uv_c2 N by 2 matrix in C2 image coordinates
% R_c2c1 the rotation from C1 frame to C2 frame, t_c2c1 the translation
% from c1 frame to c2 frame, it is the position of the camera c1 origin in
% c2 frame, n is the average normal direction of the feature point in the
% C1 frame, d is the distance from the feature to the c1 optical center
% along the minus n direction
% for normalize() input is 2 x N matrix, output is 2 by N matrix,
% project_points2() input 3 by N matrix, output 2 by N matrix
function uv_c2=rotate_with_dist_fc_c2c1(cam, uv_c1, R_c2c1, t_c2c1, n, d)

xn=normalize(uv_c1',cam.fc,cam.cc,cam.kc,cam.alpha);
xn2=((R_c2c1-(t_c2c1*n'/d)))*[xn;ones(1,size(xn,2))];
uv_c2=project_points2(xn2,zeros(3,1),zeros(3,1),cam.fc,cam.cc,cam.kc,cam.alpha);
uv_c2=uv_c2';

repo: ekfmonoslam


text:
%  qtkp_shift_solve -- Solves the problem (alpha@T1@T2@...@Tp-lambda*I)*y=c
%  where @ stands for the kronecker product for quasi-triangular matrices Ti
%
%  ::
%
% 	y = qtkp_shift_solve(T,n,c,lambda)
% 	y = qtkp_shift_solve(T,n,c,lambda,alpha)
% 	y = qtkp_shift_solve(T,[],...)
%
%  Args:
%
%     - **T** [cell array]: each element Ti, i=1,2,...,p is a
%       quasi-triangular (square) matrix
%
%     - **n** [vector|{[]}]: number of rows (or columns) for Ti,
%       i=1,2,...,p.
%
%     - **c** [vector]: right-hand side of the equality
%
%     - **lambda** [scalar]: See formula above
%
%     - **alpha** [scalar|2 x 2 matrix|{1}]: First element of the kronecker
%       product
%
%  Returns:
%     :
%
%     - **y** [N x 1 vector]: where N=n1 x n2 x ... x np
%
%  See also : utils.kronecker.kp_shift_solve
%

author: jmaih


text:
function sqi = stdsqi(signal)
%stdSQI Standard Deviation SQI
%
% Returns the normalized standard deviation of a signal on the upper
% 5% of the signals amplitude
%
%
% Input:
%   signal:         single channel (F)ECG [1xN double]
%
% Output:
%   sqi:            resulting sSQI for segment
%
% --
% fecgsyn toolbox, version 1.2, March 2017
% Released under the GNU General Public License
%
% Copyright (C) 2017  Joachim Behar & Fernando Andreotti
% Department of Engineering Science, University of Oxford
% joachim.behar@oxfordalumni.org, fernando.andreotti@eng.ox.ac.uk
%
%
% For more information visit: http://www.fecgsyn.com
%
% Referencing this work
% (SQI indices)
% Andreotti, F., Gräßer, F., Malberg, H., & Zaunseder, S. (2017). Non-Invasive Fetal ECG Signal Quality Assessment for
% Multichannel Heart Rate Estimation. IEEE Trans. Biomed. Eng., (in press).
%
% (FECGSYN Toolbox)
% Behar, J., Andreotti, F., Zaunseder, S., Li, Q., Oster, J., & Clifford, G. D. (2014). An ECG Model for Simulating
% Maternal-Foetal Activity Mixtures on Abdominal ECG Recordings. Physiol. Meas., 35(8), 1537–1550.
%
%
% Last updated : 15-03-2017
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

y = sort(signal,'descend');
nfac = median(y(1:round(0.05*length(y))));
sqi=std(signal)/nfac;
sqi(sqi>1) = 1; % limit to one, in the exceptional case it is larger than unit

end

author: fernandoandreotti


text:
function obj = removeContact(obj, contact)
    % Remove contact constraints defined for the system
    %
    % Parameters:
    % contact: the name of the contact @type cellstr

    % validate the input argument
    if ischar(contact)
        contact_name = {contact};
    elseif isa(contact, 'ContactFrame')
        contact_name = {contact.Name};
    elseif iscell(contact)
        contact_name = cell(numel(contact),1);
        for i=1:numel(contact)
            c = contact{i};
            if ischar(c)
                contact_name{i} = c;
            elseif isa(c, 'ContactFrame')
                contact_name{i} = c.Name;
            else
                error('You must provide the name of the contact or the ContactFrame object.')
            end
        end
    else
        error('You must provide the name of the contact or the ContactFrame object.')
    end

    for i=1:length(contact_name)
        name = contact_name{i};

        % remove the associated holonomic constraints
        obj.removeHolonomicConstraint(name);

        % remove the associated unilateral constraints
        % 1. friction cone
        fc_cstr = ['fc', name];
        if isfield(obj.UnilateralConstraints, fc_cstr)
            obj.removeUnilateralConstraint(fc_cstr);
        end
        % 2. ZMP constraints
        zmp_cstr = ['zmp', name];
        if isfield(obj.UnilateralConstraints, zmp_cstr)
            obj.removeUnilateralConstraint(zmp_cstr);
        end
    end
end

author: ayonga

repo: frost-dev

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


text:
 function [xs, info] = pl_iot(x, Ab, data, R, varargin)
%function [xs, info] = pl_iot(x, Ab, data, R, [options])
%|
%| Generic penalized-likelihood minimization,
%| for arbitrary negative log-likelihood with convex non-quadratic penalty,
%| via incremental optimization transfer using separable quadratic surrogates.
%|
%| cost(x) = sum_i h(data_i; [Ax]_i) + R(x),
%|
%| in
%|	x	[np 1]		initial estimate
%|	Ab	[nd np]		system matrix (see Gblock)
%|	data	{cell}		whatever data is needed for the likelihood
%|				data{1} must have size "[nb,na]", so for 3D
%|				one must use reshaper(yi, '2d').
%|	R			penalty object (see Robject.m)
%|
%| options
%|	dercurv	{function_handle} function returning derivatives and curvatures
%|				of negative log-likeihood via the call:
%|		[deriv curv] = dercurv(data, Ab{m}*x, curvtype, iblock, nblock)
%|				or choose 'trl' or 'wls' (default)
%|	niter	#		# total iterations (default: 1+0, pure OS!)
%|	os	#		how many "warmup" OS iterations (default: 1)
%|	riter	#		# of penalty subiterations (default: 3).
%|	curvtype ''		curvature type:
%|					'pc' : precomputed (default)
%|					'oc' : optimal (ensures monotonic)
%|	pixmax	(value)		upper bound on pixel values (default: infinity)
%|	pixmin	(value)		lower bound on pixel values (default: 0)
%|	isave	[]		list of iterations to archive
%|				(default: [] 'last)
%|	gi	[nd]		precomputed Ab * 1 factors.
%|	chat	(value)		verbosity (default: 0)
%|	userfun			user defined function handle (see default below)
%|
%| out
%|	xs	[np nsave]	estimates each (saved) iteration
%|	info	[niter+1 ?]	userfun output.  default is cpu times
%|
%| Copyright 2005-3-8, Jeff Fessler, University of Michigan

if nargin == 1 && streq(x, 'test')
	run_mfile_local('pwls_example')
return
end
if nargin < 4, ir_usage, end

cpu etic

Ab = block_ob(Ab, 'ensure'); % make it a block object (if not already)
nblock = block_ob(Ab, 'n');
starts = subset_start(nblock);

% defaults
arg.dercurv = 'wls';
arg.curvtype = 'pc';
arg.riter = 3;
arg.os = 1;
arg.niter = arg.os + 0; % default is pure OS, no IOT
arg.isave = [];
arg.pixmin = 0;
arg.pixmax = inf;
arg.gi = [];
arg.update_even_if_denom_0 = true;
arg.chat = false;
arg.userfun = @userfun_default;

% options
arg = vararg_pair(arg, varargin);
arg.isave = iter_saver(arg.isave, arg.niter);

if streq(arg.dercurv, 'wls')
	arg.dercurv = @wls_dercurv;
elseif streq(arg.dercurv, 'trl')
	arg.dercurv = @trl_dercurv;
end

if ~isa(arg.dercurv, 'function_handle')
	error 'dercurv not function handle?'
end

ticker(mfilename, 1, arg.niter)

[nb na] = size(data{1});

% g_i = sum_j g_ij.  caution: this requires real g_ij and g_ij >= 0
if isempty(arg.gi)
	arg.gi = sum(Ab');
	arg.gi = max(reale(arg.gi), 0); % trick: just in case...
	arg.gi = reshape(arg.gi, [nb na]);
end

% precompute likelihood-term denominator if needed
if streq(arg.curvtype, 'pc')
	[dummy curvi] = feval(arg.dercurv, data, 0, 'pc');
	ldenom = Ab' * col(arg.gi .* curvi); % one denominator shared by all subsets
else
	ldenoms = zeros(np, nblock);
end

np = length(x);
x = min(x, arg.pixmax);
x = max(x, arg.pixmin);
xs = zeros(np, length(arg.isave));
if any(arg.isave == 0)
	xs(:, arg.isave == 0) = x;
end

%info = zeros(arg.niter, ?); % trick: do not initialize since size may change

%
% precompute gradient-related state vectors, usually by OS-SPS warmup iterations
%
vvm = zeros(np,nblock);

if arg.os == 0
	for iblock=1:nblock
%		vvm(:,iblock) = Ab{ii}
		error 'not done'
	end
end

%
% SPS-OS iterations to "warm up" (often 1 or 2 iterations suffices)
%
for iter = 1:arg.os
	for iset = 1:nblock
		ticker([mfilename ': os'], [iter iset], [arg.os nblock])

		iblock = starts(iset);
		ia = iblock:nblock:na;

		li = Ab{iblock} * x;	% l=A*x "line integrals"
		li = reshape(li, nb, length(ia));
		[dhi curvi] = feval(arg.dercurv, data, li, ...
			arg.curvtype, iblock, nblock);

		if ~streq(arg.curvtype, 'pc') % on-the-fly curvatures
			ldenoms(:,iblock) = ...
				Ab{iblock}' * col(arg.gi(:,ia) .* curvi);
			ldenom = nblock * ldenoms(:,iblock);
		end

		lnum = ldenom .* x - nblock * (Ab{iblock}' * dhi(:));

		if iter == arg.os % save last gradient-related state vectors
			vvm(:,iblock) = lnum / nblock;
		end

		x = inner_update(x, lnum, ldenom, R, arg);
	end
	if any(arg.isave == iter)
		xs(:, arg.isave == iter) = x;
	end
	info(1+iter,:) = feval(arg.userfun);
	if arg.chat, printf('Range %g %g', min(x), max(x)), end
end

%
% At this point we have initialized vv_m and ldenoms
% and initialized x based on the last subset.
% It is often logical to update x using all subsets now.
% This is "almost free" since no new likelihood gradients are used.
%
vv = sum(vvm,2);
if streq(arg.curvtype, 'pc')
	ldenoms = ldenom / nblock;
else
	ldenom = sum(ldenoms,2);
end
if 1
	x = inner_update(x, vv, ldenom, R, arg);
	if any(arg.isave == -1)
		xs(:, arg.isave == -1) = x;
	end
end

%
% IOT iterations
%
for iter = (1+arg.os):arg.niter
	for iset = 1:nblock
		ticker([mfilename ': iot'], [iter iset], [arg.niter nblock])

		iblock = starts(iset);
		ia = iblock:nblock:na;

		li = Ab{iblock} * x;	% l=A*x "line integrals"
		li = reshape(li, nb, length(ia));
		[dhi curvi] = feval(arg.dercurv, data, li, ...
			arg.curvtype, iblock, nblock);

		if ~streq(arg.curvtype, 'pc')
			ldenom = ldenom - ldenoms(:,iblock);
			ldenoms(:,iblock) = ...
				Ab{iblock}' * col(arg.gi(:,ia) .* curvi);
			ldenom = ldenom + ldenoms(:,iblock);
		end

		vv = vv - vvm(:,iblock);
		if streq(arg.curvtype, 'pc')
			vvm(:,iblock) = ldenoms .* x - (Ab{iblock}' * dhi(:));
		else
			vvm(:,iblock) = ldenoms(:,iblock) .* x - (Ab{iblock}' * dhi(:));
		end
		vv = vv + vvm(:,iblock);

		x = inner_update(x, vv, ldenom, R, arg);
	end

	if any(arg.isave == iter)
		xs(:, arg.isave == iter) = x;
	end
	info(1+iter,:) = feval(arg.userfun);
	if arg.chat, printf('Range %g %g', min(x), max(x)), end
end

%
% inner_update()
% given SPS for likelihood, update x using regularization subiterations
%
function x = inner_update(x, lnum, ldenom, R, arg)

for ii=1:arg.riter
	if isempty(R)
		num = lnum;
		rdenom = 0;
	else
		rdenom = R.denom(R, x);
		num = lnum - R.cgrad(R, x) + rdenom .* x;
	end

	% update
	if arg.update_even_if_denom_0
		x = num ./ (ldenom + rdenom);
	else
		old = x(ldenom == 0);
		x = num ./ (ldenom + rdenom);
		x(ldenom == 0) = old;
	end
	x = max(x, arg.pixmin);	% lower bound
	x = min(x, arg.pixmax);	% upper bound
end

% default user function.
% using this evalin('caller', ...) trick, one can compute anything of interest
function out = userfun_default
%x = evalin('caller', 'x');
out = [cpu('etoc')];

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%% Filename and Format string
% This code loads in the contents of the data file with textscan a block at a
% time and converts the cellarray from textscan into a double array

filename='waferdata.csv';
FormatString=['%*f%*f%*f%*f' repmat('%f',1,9) repmat('%*f',1,7)];

%% Parameters
Headers=3;
NumLines=1e6; % Total number of lines to read
BlockLines=10000; % Size of block
NumBlocks=NumLines/BlockLines; % Number of blocks

%% Open file
fid = fopen(filename); % Open file

%% Preallocate data
data=zeros(NumLines,9); % Pre-allocate space for data

%% Get Headers
cellchunk=textscan(fid,'%s',3,'delimiter','\n'); % Get first 3 lines

%% Read in blocks
for Block=1:NumBlocks
    data((Block-1)*BlockLines+(1:BlockLines),:) = ...
	cell2mat(textscan(fid,FormatString,BlockLines,'delimiter',','));
    disp(['Block ' num2str(Block) ]); % Display current block to show progress
end

%% Close file
fclose(fid);


text:
%ROTZ Rotation about Z axis
%
% R = ROTZ(THETA) is an SO(3) rotation matrix (3x3) representing a rotation of THETA
% radians about the z-axis.
%
% R = ROTZ(THETA, 'deg') as above but THETA is in degrees.
%
% See also ROTX, ROTY, ANGVEC2R, ROT2.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB 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.
%
% RTB 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 RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function R = rotz(t, deg)
    if nargin > 1 && strcmp(deg, 'deg')
        t = t *pi/180;
    end

    ct = cos(t);
    st = sin(t);
    R = [
        ct  -st  0
        st   ct  0
        0    0   1
        ];


text:
function x = p19_start ( i, x )

%*****************************************************************************80
%
%% P19_START returns a starting point for problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the starting point.
%
%    Output, real X, the starting point.
%
  if ( i == 1 )
    x = 0.0;
  elseif ( i == 2 )
    x = 1.0;
  elseif ( i == 3 )
    x = 0.5;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P19_START - Fatal error!\n' );
    fprintf ( 1, '  Illegal starting point index = %d\n', i );
    error ( 'P19_START - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = log(f)
%LOG   Natural logarithm of a CHEBFUN3.
%   LOG(F) is the natural logarithm of F. This function gives an error
%   message if the function passes through or becomes numerically close to
%   zero.
%
% See also CHEBFUN3/COMPOSE.

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

% Empty check:
if ( isempty(f) )
    return
end

if ( isreal(f) )
    % Positive/negative test.
    ss = singleSignTest(f);  % Returns TRUE if there is no sign change.
    if ( ~ss )
        error('CHEBFUN:CHEBFUN3:log:notSmooth', ...
            'Sign change detected. Unable to represent the result.');
    end
end

f = compose(f, @log);

end


text:
%% Constrained Total Variation Projection with ADDM
% Test for ADMM algorithm convergence on a randomized example.

addpath('../');
addpath('../toolbox/');

%%
% We want to solve a problem of the form

%%
%   |min_{x in C, norm(K*x,1)<=1} -<x,w> + epsilon/2*norm(x)^2

%%
% Here |K| is some sort of "gradient" operator (here we use a random
% matrix) and C

%%
% This can be written as

%%
%   |min_x F(K*x) + G(x)|

%%
% where |F(u) = i_{norm(u,1)<=1}| and
%  |G(x) = i_{C}(x) - <f,w> + epsilon/2*norm(f)^2|.

%%
% Inner product shortcut.

dotp = @(u,v)sum(u(:).*v(:));

%%
% Dimension of the problem.

n = 300;

%%
% Number of computed "gradient".

p = 500;

%%
% Number of affine constraint

r = 10;

%%
% Regularization.

epsilon = .1;

%%
% Gradient operator.

K = randn(p,n);

%%
% Linear function to optimize.

w = randn(n,1);

%%
% Constraint operator.

A = randn(r,n);
y = randn(r,1)*0;

%%
% Projector on A*f=y.

pA = A'*(A*A')^(-1);
ProjC = @(f)f + pA*(y-A*f);

%%
% Projection on L1 ball.

ProxF = @(u,rho)perform_l1ball_projection(u,1);
ProxFS = compute_dual_prox(ProxF);

%%
% Proximal operator of G.

ProxG = @(x,tau)ProjC( (x+tau*w)/(1+tau*epsilon) );

%%
% Callback to record information during the iterations.

F = @(x)-dotp(x,w) + epsilon/2*norm(x(:))^2;
Constr = @(x)norm(K*x,1);
options.report = @(x)struct('F', F(x), 'Constr', Constr(x));

%%
% Run the algorihtm.

options.niter = 5000;
[f,R] = perform_admm(zeros(n,1), K, K', ProxFS, ProxG, options);

%%
% Retrieve the F and constraint function values.

f = s2v(R,'F');
constr = s2v(R,'Constr');

%%
% Display.

clf;
subplot(2,1,1);
plot(f);
title('Energy');
subplot(2,1,2);
plot(log10(abs(constr-1)));
title('Constraint');

author: gpeyre

repo: matlab-toolboxes


text:
function VAR = MoveCluster(VAR,i_fish,sourceIDs,destIDs)
nClus = length(VAR(i_fish).ClusGroup{sourceIDs{1}});
clus_list = intersect(1:nClus,sourceIDs{2});

Source = VAR(i_fish).ClusGroup{sourceIDs{1}}(clus_list);

dest_clus_list = destIDs{2}(1:length(clus_list));
VAR(i_fish).ClusGroup{destIDs{1}}(dest_clus_list) = Source;

VAR(i_fish).ClusGroup{sourceIDs{1}}(clus_list) = [];
end

author: xiuyechen

repo: FishExplorer


text:
function [MSE, Entropy, SNRdB, sdV] = QuantSNR (Yq, Xq, FPDF, ScaleF)
% Calculate the SNR (dB) for a quantizer defined by a table for a
% signal with a given probability density function.
%
% Yq - Nlev quantizer output levels (in ascending order)
% Xq - Nlev-1 quantizer decision levels (in ascending order)
% FPDF - Cell array of function handles {Farea, Fmean, Fvar}
% Farea - Pointer to a function to calculate the integral of a probability
%   density function in an interval.
%                 b
%   Farea(a,b) = Int p(x) dx,
%                 a
%   where (a,b) is the interval and p(x) is the probability density
%   function.
% Fmean - Pointer to a function to calculate the mean of a probability
%   density function in an interval.
%                 b
%   Fmean(a,b) = Int x p(x) dx,
%                 a
%   where (a,b) is the interval and p(x) is the probability density
%   function.
% Fvar - Pointer to a function to calculate the second moment of a
%   probability density function in an interval.
%                b
%   Fvar(a,b) = Int x^2 p(x) dx,
%                a
%   where (a,b) is the interval and p(x) is the probability density
%   function.
% ScaleF - Scaling factor for the signal (optional, default 1)

% MSE - Resulting mean square error
% Entropy - Resulting entropy (bits)
% SNRdB - Resulting signal-to-noise ratio dB
% sdV - Peak to average ratio for the quantizer. The peak level (overload)
%   V of the quantizer is defined here as the the distance from the mean to
%   the largest output level plus a half step size. Then
%     sdV = sd / V,
%   where sd is the standard deviation of the signal.

if (nargin <= 3)
  ScaleF = 1;
end

% The scale factor represents a scaling of the signal
% For convenience we scale the quantizer by 1/ScaleF and later multipy
% the computed MSE by ScaleF^2 achieve the same effect as scaling the
% signal

% Quantizers are really a combination of two operations
%  - determine an interval k, such that Xq(k-1) <= x < Xq(k)
%  - output a value Yq(k)
% Consider a scaled signal g*x, then k can be determined as the value
% which satisfies
%   Xq(k-1)/g <= x < Xq(k)/g
% The quantization error for the scaled signal is
%   e = g*x - Yq(k)
%     = g * (x - Yq(k)/g)
% The bracketed term is the error for a signal x (not scaled) applied
% to a quantizer with break points Xq(.)/g and output levels Yq(.)/g.

Nlev = length(Yq);
if (ScaleF ~=1)
  if (Nlev > 1)
    Xq = Xq / ScaleF;
  end
  Yq = Yq / ScaleF;
end

% Calculate the mean square error
MSE = ScaleF^2 * QuantMSE (Yq, Xq, FPDF);

% Calculate the entropy
if (nargout >= 2)
  Entropy = QuantEntropy (Xq, FPDF{1});
end

% Calculate the SNR in dB
if (nargout >= 3)
  Fmean = FPDF{2};
  Fvar = FPDF{3};
  Xmean = feval(Fmean, -Inf, Inf);
  Svar = feval(Fvar, -Inf, Inf) - Xmean^2;
  SNRdB = 10 * log10(ScaleF^2 * Svar / MSE);
end

% Calculate the peak to average ratio for the quantizer
if (nargout >= 4)
  if (Nlev > 1)
    VL = (Yq(1)-Xq(1)) + Yq(1);
    VH = (Yq(end)-Xq(end)) + Yq(end);
    V = max(abs(VL-Xmean), abs(VH-Xmean));
  else
    V = Inf;
  end
  sdV = sqrt(Svar) / V;
end

return