text:
function submit(partId, webSubmit)
%SUBMIT Submit your code and output to the pgm-class servers
%   SUBMIT() will connect to the pgm-class server and submit your solution
%   There is no penalty for submitting, so go ahead and try this!
%
%   If this function does not work for you, use the web-submission mechanism.
%   Call the submitWeb function, which will produce a file for the part you
%   wish to submit. Then, submit the file to the class servers using the
%   "Web Submission" button on the Programming Assignments page on the course
%   website.
%
%   webSubmit is a boolean variable that specifies whether to prepare
%   a file for web-submission (if webSubmit = 1) or to submit directly
%   to the server (if webSubmit = 0, or is not specified); you should call
%   submitWeb if you want to do a web-submission.

  warning off all;

  fprintf('==\n== [pgm-class] Submitting Solutions | Programming Assignment %s\n==\n', ...
          homework_id());
  if ~exist('partId', 'var') || isempty(partId)
    partId = promptPart();
  end

  if ~exist('webSubmit', 'var') || isempty(webSubmit)
    webSubmit = 0; % submit directly by default
  end

  % Check valid partId
  partNames = validParts();
  if ~isValidPartId(partId)
    fprintf('!! Invalid assignment part selected.\n');
    fprintf('!! Expected an integer from 1 to %d.\n', numel(partNames) + 1);
    fprintf('!! Submission Cancelled\n');
    return
  end

  if ~exist('228_login_data.mat','file')
    [login password] = loginPrompt();
    save('228_login_data.mat','login','password');
  else
    load('228_login_data.mat');
    [login password] = quickLogin(login,password);
    save('228_login_data.mat','login','password');
  end
  if isempty(login)
    fprintf('!! Submission Cancelled\n');
    return
  end

  fprintf('\n== Connecting to pgm-class ... ');
  if exist('OCTAVE_VERSION')
    fflush(stdout);
  end

  % Setup submit list
  if partId == numel(partNames) + 1
    submitParts = 1:numel(partNames);
  else
    submitParts = [partId];
  end

  for s = 1:numel(submitParts)
    % Submit this part
    partId = submitParts(s);

    samplePartId = 2 * (partId - 1) + 1;
    testPartId = samplePartId + 1;

    for thisPartId = [samplePartId testPartId]
        if (~webSubmit) % submit directly to server

            [login, ch, signature, auxstring] = getChallenge(login, thisPartId);
            if isempty(login) || isempty(ch) || isempty(signature)
                % Some error occured, error string in first return element.
                fprintf('\n!! Error: %s\n\n', login);
                return
            end

            % Attempt Submission with Challenge
            ch_resp = challengeResponse(login, password, ch);

            [result, str] = submitSolution(login, ch_resp, thisPartId, ...
                output(thisPartId, auxstring), source(partId), signature);

            if (~isTest(thisPartId))
                partName = partNames{partId};
            else
                partName = [partNames{partId} ' (test)'];
            end

            fprintf('\n== [pgm-class] Submitted Assignment %s - Part %d - %s\n', ...
                homework_id(), partId, partName);
            fprintf('== %s\n', strtrim(str));

            if exist('OCTAVE_VERSION')
                fflush(stdout);
            end

        else % make web submission files

            [result] = submitSolutionWeb(login, thisPartId, output(thisPartId), ...
                source(partId));
            result = base64encode(result);

            if (~isTest(thisPartId))
                partType = 'sample';
            else
                partType = 'test';
            end

            fprintf('\nSave as submission file [submit_pa%s_part%d_%s.txt (enter to accept default)]:', ...
                homework_id(), partId, partType);
            saveAsFile = input('', 's');
            if (isempty(saveAsFile))
                saveAsFile = sprintf('submit_pa%s_part%d_%s.txt', homework_id(), partId, partType);
            end

            fid = fopen(saveAsFile, 'w');
            if (fid)
                fwrite(fid, result);
                fclose(fid);
                fprintf('\nSaved your solutions to %s.\n\n', saveAsFile);
                fprintf(['You can now submit your solutions using the ' ...
                    'Web Submission button\non the Assignments page\n']);
            else
                fprintf('Unable to save to %s\n\n', saveAsFile);
                fprintf(['You can create a submission file by saving the \n' ...
                    'following text in a file: (press enter to continue)\n\n']);
                pause;
                fprintf(result);
            end

        end
    end

  end

end

% ================== CONFIGURABLES FOR EACH HOMEWORK ==================

function id = homework_id()
  id = '2';
end

function [partNames] = validParts()
  partNames = { 'phenotypeGivenGenotypeMendelianFactor', ...
                'phenotypeGivenGenotypeFactor', ...
                'genotypeGivenAlleleFreqsFactor', ...
                'genotypeGivenParentsGenotypesFactor', ...
                'constructGeneticNetwork', ...
                'phenotypeGivenCopiesFactor', ...
                'constructDecoupledGeneticNetwork', ...
                'constructSigmoidPhenotypeFactor'
              };
end

function srcs = sources()
  % Separated by part
  srcs = { { 'phenotypeGivenGenotypeMendelianFactor.m' }, ...
           { 'phenotypeGivenGenotypeFactor.m' }, ...
           { 'genotypeGivenAlleleFreqsFactor.m' }, ...
           { 'genotypeGivenParentsGenotypesFactor.m' }, ...
           { 'constructGeneticNetwork.m' }, ...
           { 'phenotypeGivenCopiesFactor.m' }, ...
           { 'constructDecoupledGeneticNetwork.m' }, ...
           { 'constructSigmoidPhenotypeFactor.m' }
         };
end

% specifies which parts are test parts
function result = isTest(partId)
  if (bitand(partId, 0))
    result = true;
  else
    result = false;
  end
end

function out = output(partId, auxstring)

  if partId == 1
    isDominant = 1;
    genotypeVar = 1;
    phenotypeVar = 3;
    pggmf(1) = phenotypeGivenGenotypeMendelianFactor(isDominant, genotypeVar, phenotypeVar);
    pggmf(1) = SortFactorLaterVars(pggmf(1));
    out = SerializeFactorsFg(pggmf);
  elseif partId == 2
    isDominant = 0;
    genotypeVar = 1;
    phenotypeVar = 3;
    pggmf(1) = phenotypeGivenGenotypeMendelianFactor(isDominant, genotypeVar, phenotypeVar);
    pggmf(1) = SortFactorLaterVars(pggmf(1));
    out = SerializeFactorsFg(pggmf);
  elseif partId == 3
    alphaList = [0.8; 0.6; 0.1];
    genotypeVar = 1;
    phenotypeVar = 3;
    pggf(1) = phenotypeGivenGenotypeFactor(alphaList, genotypeVar, phenotypeVar);
    pggf(1) = SortFactorLaterVars(pggf(1));
    out = SerializeFactorsFg(pggf);
  elseif partId == 4
    alphaList = [0.2; 0.5; 0.9];
    genotypeVar = 1;
    phenotypeVar = 3;
    pggf(1) = phenotypeGivenGenotypeFactor(alphaList, genotypeVar, phenotypeVar);
    pggf(1) = SortFactorLaterVars(pggf(1));
    out = SerializeFactorsFg(pggf);
  elseif partId == 5
    alleleFreqs = [0.1; 0.9];
    genotypeVar = 1;
    ggaff(1) = genotypeGivenAlleleFreqsFactor(alleleFreqs, genotypeVar);
    out = SerializeFactorsFg(ggaff);
  elseif partId == 6
    alleleFreqs = [0.98; 0.02];
    genotypeVar = 1;
    ggaff(1) = genotypeGivenAlleleFreqsFactor(alleleFreqs, genotypeVar);
    out = SerializeFactorsFg(ggaff);
  elseif partId == 7
    alleleList = {'F', 'f'};
    genotypeVarChild = 3;
    genotypeVarParentOne = 1;
    genotypeVarParentTwo = 2;
    ggpgf(1) = genotypeGivenParentsGenotypesFactor(length(alleleList), genotypeVarChild, genotypeVarParentOne, genotypeVarParentTwo);
    ggpgf(1) = SortFactorLaterVars(ggpgf(1));
    out = SerializeFactorsFg(ggpgf);
  elseif partId == 8
    alleleList = {'A', 'a', 'n'};
    genotypeVarChild = 3;
    genotypeVarParentOne = 1;
    genotypeVarParentTwo = 2;
    ggpgf(1) = genotypeGivenParentsGenotypesFactor(length(alleleList), genotypeVarChild, genotypeVarParentOne, genotypeVarParentTwo);
    ggpgf(1) = SortFactorLaterVars(ggpgf(1));
    out = SerializeFactorsFg(ggpgf);
  elseif partId == 9
    alleleFreqs = [0.1; 0.9];
    alphaList = [0.8; 0.6; 0.1];
    pedigree = struct('parents', [0,0;1,3;0,0;1,3;2,6;0,0;2,6;4,9;0,0]);
    pedigree.names = {'Ira','James','Robin','Eva','Jason','Rene','Benjamin','Sandra','Aaron'};
    cgn = constructGeneticNetwork(pedigree, alleleFreqs, alphaList);
    for i = 1:length(cgn)
        cgn(i) = SortFactorLaterVars(cgn(i));
    end
    out = SerializeFactorsFg(sortStruct(cgn));
  elseif partId == 10
    alleleFreqs = [0.1; 0.9];
    alphaList = [0.8; 0.6; 0.1];
    pedigree = struct('parents', [0,0;0,0;2,1;0,0;2,1;0,0;0,0;3,4;5,7;5,6]);
	pedigree.names = {'Alan','Vivian','Alice','Larry','Beth','Henry','Leon','Frank','Amy', 'Martin'};
    cgn = constructGeneticNetwork(pedigree, alleleFreqs, alphaList);
    for i = 1:length(cgn)
        cgn(i) = SortFactorLaterVars(cgn(i));
    end
    out = SerializeFactorsFg(sortStruct(cgn));
  elseif partId == 11
    alphaListThree = [0.8; 0.6; 0.1; 0.5; 0.05; 0.01];
    phenotypeVar = 3;
    genotypeVarMotherCopy = 1;
    genotypeVarFatherCopy = 2;
    numAlleles = 3;
    pgpaf(1) = phenotypeGivenCopiesFactor(alphaListThree, numAlleles, genotypeVarMotherCopy, genotypeVarFatherCopy, phenotypeVar);
    pgpaf(1) = SortFactorLaterVars(pgpaf(1));
    out = SerializeFactorsFg(pgpaf);
  elseif partId == 12
    alphaListThree = [0.001; 0.009; 0.3; 0.2; 0.75; 0.95];
    phenotypeVar = 3;
    genotypeVarMotherCopy = 1;
    genotypeVarFatherCopy = 2;
    numAlleles = 3;
    pgpaf(1) = phenotypeGivenCopiesFactor(alphaListThree, numAlleles, genotypeVarMotherCopy, genotypeVarFatherCopy, phenotypeVar);
    pgpaf(1) = SortFactorLaterVars(pgpaf(1));
    out = SerializeFactorsFg(pgpaf);
  elseif partId == 13
    pedigree = struct('parents', [0,0;1,3;0,0;1,3;2,6;0,0;2,6;4,9;0,0]);
    pedigree.names = {'Ira','James','Robin','Eva','Jason','Rene','Benjamin','Sandra','Aaron'};
    alphaListThree = [0.8; 0.6; 0.1; 0.5; 0.05; 0.01];
    alleleFreqsThree = [0.1; 0.7; 0.2];
    cdgn = constructDecoupledGeneticNetwork(pedigree, alleleFreqsThree, alphaListThree);
    for i = 1:length(cdgn)
        cdgn(i) = SortFactorLaterVars(cdgn(i));
    end
    out = SerializeFactorsFg(sortStruct(cdgn));
  elseif partId == 14
    pedigree = struct('parents', [0,0;0,0;2,1;0,0;2,1;0,0;0,0;3,4;5,7;5,6]);
	pedigree.names = {'Alan','Vivian','Alice','Larry','Beth','Henry','Leon','Frank','Amy', 'Martin'};
    alphaListThree = [0.8; 0.6; 0.1; 0.5; 0.05; 0.01];
    alleleFreqsThree = [0.1; 0.7; 0.2];
    cdgn = constructDecoupledGeneticNetwork(pedigree, alleleFreqsThree, alphaListThree);
    for i = 1:length(cdgn)
        cdgn(i) = SortFactorLaterVars(cdgn(i));
    end
    out = SerializeFactorsFg(sortStruct(cdgn));
  elseif partId == 15
    alleleWeights = {[3, -3], [0.9, -0.8]};
    phenotypeVar = 3;
    geneCopyVarParentOneList = [1; 2];
    geneCopyVarParentTwoList = [4; 5];
    cspf(1) = constructSigmoidPhenotypeFactor(alleleWeights, geneCopyVarParentOneList, geneCopyVarParentTwoList, phenotypeVar);
    cspf(1) = SortFactorLaterVars(cspf(1));
    out = SerializeFactorsFg(cspf);
  elseif partId == 16
    alleleWeights = {[0.01, -.2], [1, -.5]};
    phenotypeVar = 3;
    geneCopyVarParentOneList = [1; 2];
    geneCopyVarParentTwoList = [4; 5];
    cspf(1) = constructSigmoidPhenotypeFactor(alleleWeights, geneCopyVarParentOneList, geneCopyVarParentTwoList, phenotypeVar);
    cspf(1) = SortFactorLaterVars(cspf(1));
    out = SerializeFactorsFg(cspf);
  end
end

function out = SerializeFactorsFg(F)
% Serializes a factor struct array into the .fg format for libDAI
% http://cs.ru.nl/~jorism/libDAI/doc/fileformats.html
%
% To avoid incompatibilities with EOL markers, make sure you write the
% string to a file using the appropriate file type ('wt' for windows, 'w'
% for unix)

lines = cell(5*numel(F) + 1, 1);

lines{1} = sprintf('%d\n', numel(F));
lineIdx = 2;
for i = 1:numel(F)
  lines{lineIdx} = sprintf('\n%d\n', numel(F(i).var));
  lineIdx = lineIdx + 1;

  lines{lineIdx} = sprintf('%s\n', num2str(F(i).var(:)')); % ensure that we put in a row vector
  lineIdx = lineIdx + 1;

  lines{lineIdx} = sprintf('%s\n', num2str(F(i).card(:)')); % ensure that we put in a row vector
  lineIdx = lineIdx + 1;

  lines{lineIdx} = sprintf('%d\n', numel(F(i).val));
  lineIdx = lineIdx + 1;

  % Internal storage of factor vals is already in the same indexing order
  % as what libDAI expects, so we don't need to convert the indices.
  vals = [0:(numel(F(i).val) - 1); F(i).val(:)'];
  lines{lineIdx} = sprintf('%d %0.8g\n', vals);
  lineIdx = lineIdx + 1;
end

out = sprintf('%s', lines{:});

end

function out = SerializePedigree(pedigree)
% SerializePedigree Serializes a pedigree into a text string.

numPeople = numel(pedigree.names);

lines = cell(2*numPeople + 1, 1);

lines{1} = sprintf('%d\n', numPeople);
lineIdx = 2;
for i = 1:numPeople
  lines{lineIdx} = sprintf('%s\n', pedigree.names{i});
  lineIdx = lineIdx + 1;
end

for i = 1:numPeople
  lines{lineIdx} = sprintf('%d %d\n', pedigree.parents(i, 1), pedigree.parents(i, 2));
  lineIdx = lineIdx + 1;
end

out = sprintf('%s', lines{:});
end

function G = SortFactorLaterVars(F)
% SortFactorLaterVars Returns the input factor with a sorted variable
% ordering but leaves the first variable alone.
%
% SortFactorLaterVars(F) returns a new factor that contains exactly the same
% variables and values but in which the variables are ordered differently
% (the ordering is as specified in the .var field). Specifically, the
% returned factor will have its first variable  ordered in ascending order.
%
% For instance, if F.var = [3 1 2], and G = SortFactorVars(F), tnen
% G.var = [1 2 3], and G.card, G.val will have been appropriately reordered
% (from F.card and F.val) with respect to this sorted variable ordering.
%

if (numel(F.var) > 1)
  laterVars = F.var(2:end);
  [slVars, lorder] = sort(laterVars);
  G.var = [F.var(1) slVars];

  order = [1 (lorder+1)];

  G.card = F.card(order);
  G.val = zeros(numel(F.val), 1);

  assignmentsInF = IndexToAssignment(1:numel(F.val), F.card);
  assignmentsInG = assignmentsInF(:, order);
  G.val(AssignmentToIndex(assignmentsInG, G.card)) = F.val;

else % nothing to do
  G = F;
end

end

function sListSorted = sortStruct(sList)

for i = 1 : length(sList)
  S{i} = sprintf('%d%d%f', sList(i).var, sList(i).card, sList(i).val);
end

[S I] = sort(S);
sListSorted = sList(I);

end

%%% Remove -staging when you deploy!
function url = challenge_url()
  url = 'http://www.coursera.org/pgm/assignment/challenge';
end

%%% Remove -staging when you deploy!
function url = submit_url()
  url = 'http://www.coursera.org/pgm/assignment/submit';
end

% ========================= CHALLENGE HELPERS =========================

function src = source(partId)
  src = '';
  src_files = sources();
  if partId <= numel(src_files)
      flist = src_files{partId};
      for i = 1:numel(flist)
          fid = fopen(flist{i});
          if (fid == -1)
            error('Error opening %s (is it missing?)', flist{i});
          end
          line = fgets(fid);
          while ischar(line)
            src = [src line];
            line = fgets(fid);
          end
          fclose(fid);
          src = [src '||||||||'];
      end
  end
end

function ret = isValidPartId(partId)
  partNames = validParts();
  ret = (~isempty(partId)) && (partId >= 1) && (partId <= numel(partNames) + 1);
end

function partId = promptPart()
  fprintf('== Select which part(s) to submit:\n', ...
          homework_id());
  partNames = validParts();
  srcFiles = sources();
  for i = 1:numel(partNames)
    fprintf('==   %d) %s [', i, partNames{i});
    fprintf(' %s ', srcFiles{i}{:});
    fprintf(']\n');
  end
  fprintf('==   %d) All of the above \n==\nEnter your choice [1-%d]: ', ...
          numel(partNames) + 1, numel(partNames) + 1);
  selPart = input('', 's');
  partId = str2num(selPart);
  if ~isValidPartId(partId)
    partId = -1;
  end
end

function [email,ch,signature,auxstring] = getChallenge(email, part)
  str = urlread(challenge_url(), 'post', {'email_address', email, 'assignment_part_sid', [homework_id() '-' num2str(part)], 'response_encoding', 'delim'});

  str = strtrim(str);
  r = struct;
  while(numel(str) > 0)
    [f, str] = strtok (str, '|');
    [v, str] = strtok (str, '|');
    r = setfield(r, f, v);
  end

  email = getfield(r, 'email_address');
  ch = getfield(r, 'challenge_key');
  signature = getfield(r, 'state');
  auxstring = getfield(r, 'challenge_aux_data');
end

function [result, str] = submitSolutionWeb(email, part, output, source)

  result = ['{"assignment_part_sid":"' base64encode([homework_id() '-' num2str(part)], '') '",' ...
            '"email_address":"' base64encode(email, '') '",' ...
            '"submission":"' base64encode(output, '') '",' ...
            '"submission_aux":"' base64encode(source, '') '"' ...
            '}'];
  str = 'Web-submission';
end

function [result, str] = submitSolution(email, ch_resp, part, output, ...
                                        source, signature)

  params = {'assignment_part_sid', [homework_id() '-' num2str(part)], ...
            'email_address', email, ...
            'submission', base64encode(output, ''), ...
            'submission_aux', base64encode(source, ''), ...
            'challenge_response', ch_resp, ...
            'state', signature};

  str = urlread(submit_url(), 'post', params);

  % Parse str to read for success / failure
  result = 0;

end

% =========================== LOGIN HELPERS ===========================

function [login password] = loginPrompt()
  % Prompt for password
  [login password] = basicPrompt();

  if isempty(login) || isempty(password)
    login = []; password = [];
  end
end

function [login password] = basicPrompt()
  login = input('Login (Email address): ', 's');
  password = input('Submission Password (from Assignments page): ', 's');
end

function [login password] = quickLogin(login,password)
  cont_token = input(['You are currently logged in as ' login '. Is this you? (y/n - type n to reenter password)'],'s');
  if(cont_token(1)=='Y'||cont_token(1)=='y')
    return;
  else
    [login password] = loginPrompt();
  end
end

function [str] = challengeResponse(email, passwd, challenge)
  str = sha1([challenge passwd]);
end

% =============================== SHA-1 ================================

function hash = sha1(str)

  % Initialize variables
  h0 = uint32(1732584193);
  h1 = uint32(4023233417);
  h2 = uint32(2562383102);
  h3 = uint32(271733878);
  h4 = uint32(3285377520);

  % Convert to word array
  strlen = numel(str);

  % Break string into chars and append the bit 1 to the message
  mC = [double(str) 128];
  mC = [mC zeros(1, 4-mod(numel(mC), 4), 'uint8')];

  numB = strlen * 8;
  if exist('idivide')
    numC = idivide(uint32(numB + 65), 512, 'ceil');
  else
    numC = ceil(double(numB + 65)/512);
  end
  numW = numC * 16;
  mW = zeros(numW, 1, 'uint32');

  idx = 1;
  for i = 1:4:strlen + 1
    mW(idx) = bitor(bitor(bitor( ...
                  bitshift(uint32(mC(i)), 24), ...
                  bitshift(uint32(mC(i+1)), 16)), ...
                  bitshift(uint32(mC(i+2)), 8)), ...
                  uint32(mC(i+3)));
    idx = idx + 1;
  end

  % Append length of message
  mW(numW - 1) = uint32(bitshift(uint64(numB), -32));
  mW(numW) = uint32(bitshift(bitshift(uint64(numB), 32), -32));

  % Process the message in successive 512-bit chs
  for cId = 1 : double(numC)
    cSt = (cId - 1) * 16 + 1;
    cEnd = cId * 16;
    ch = mW(cSt : cEnd);

    % Extend the sixteen 32-bit words into eighty 32-bit words
    for j = 17 : 80
      ch(j) = ch(j - 3);
      ch(j) = bitxor(ch(j), ch(j - 8));
      ch(j) = bitxor(ch(j), ch(j - 14));
      ch(j) = bitxor(ch(j), ch(j - 16));
      ch(j) = bitrotate(ch(j), 1);
    end

    % Initialize hash value for this ch
    a = h0;
    b = h1;
    c = h2;
    d = h3;
    e = h4;

    % Main loop
    for i = 1 : 80
      if(i >= 1 && i <= 20)
        f = bitor(bitand(b, c), bitand(bitcmp(b), d));
        k = uint32(1518500249);
      elseif(i >= 21 && i <= 40)
        f = bitxor(bitxor(b, c), d);
        k = uint32(1859775393);
      elseif(i >= 41 && i <= 60)
        f = bitor(bitor(bitand(b, c), bitand(b, d)), bitand(c, d));
        k = uint32(2400959708);
      elseif(i >= 61 && i <= 80)
        f = bitxor(bitxor(b, c), d);
        k = uint32(3395469782);
      end

      t = bitrotate(a, 5);
      t = bitadd(t, f);
      t = bitadd(t, e);
      t = bitadd(t, k);
      t = bitadd(t, ch(i));
      e = d;
      d = c;
      c = bitrotate(b, 30);
      b = a;
      a = t;

    end
    h0 = bitadd(h0, a);
    h1 = bitadd(h1, b);
    h2 = bitadd(h2, c);
    h3 = bitadd(h3, d);
    h4 = bitadd(h4, e);

  end

  hash = reshape(dec2hex(double([h0 h1 h2 h3 h4]), 8)', [1 40]);

  hash = lower(hash);

end

function ret = bitadd(iA, iB)
  ret = double(iA) + double(iB);
  ret = bitset(ret, 33, 0);
  ret = uint32(ret);
end

function ret = bitrotate(iA, places)
  t = bitshift(iA, places - 32);
  ret = bitshift(iA, places);
  ret = bitor(ret, t);
end

% =========================== Base64 Encoder ============================
% Thanks to Peter John Acklam
%

function y = base64encode(x, eol)
%BASE64ENCODE Perform base64 encoding on a string.
%
%   BASE64ENCODE(STR, EOL) encode the given string STR.  EOL is the line ending
%   sequence to use; it is optional and defaults to '\n' (ASCII decimal 10).
%   The returned encoded string is broken into lines of no more than 76
%   characters each, and each line will end with EOL unless it is empty.  Let
%   EOL be empty if you do not want the encoded string broken into lines.
%
%   STR and EOL don't have to be strings (i.e., char arrays).  The only
%   requirement is that they are vectors containing values in the range 0-255.
%
%   This function may be used to encode strings into the Base64 encoding
%   specified in RFC 2045 - MIME (Multipurpose Internet Mail Extensions).  The
%   Base64 encoding is designed to represent arbitrary sequences of octets in a
%   form that need not be humanly readable.  A 65-character subset
%   ([A-Za-z0-9+/=]) of US-ASCII is used, enabling 6 bits to be represented per
%   printable character.
%
%   Examples
%   --------
%
%   If you want to encode a large file, you should encode it in chunks that are
%   a multiple of 57 bytes.  This ensures that the base64 lines line up and
%   that you do not end up with padding in the middle.  57 bytes of data fills
%   one complete base64 line (76 == 57*4/3):
%
%   If ifid and ofid are two file identifiers opened for reading and writing,
%   respectively, then you can base64 encode the data with
%
%      while ~feof(ifid)
%         fwrite(ofid, base64encode(fread(ifid, 60*57)));
%      end
%
%   or, if you have enough memory,
%
%      fwrite(ofid, base64encode(fread(ifid)));
%
%   See also BASE64DECODE.

%   Author:      Peter John Acklam
%   Time-stamp:  2004-02-03 21:36:56 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   if isnumeric(x)
      x = num2str(x);
   end

   % make sure we have the EOL value
   if nargin < 2
      eol = sprintf('\n');
   else
      if sum(size(eol) > 1) > 1
         error('EOL must be a vector.');
      end
      if any(eol(:) > 255)
         error('EOL can not contain values larger than 255.');
      end
   end

   if sum(size(x) > 1) > 1
      error('STR must be a vector.');
   end

   x   = uint8(x);
   eol = uint8(eol);

   ndbytes = length(x);                 % number of decoded bytes
   nchunks = ceil(ndbytes / 3);         % number of chunks/groups
   nebytes = 4 * nchunks;               % number of encoded bytes

   % add padding if necessary, to make the length of x a multiple of 3
   if rem(ndbytes, 3)
      x(end+1 : 3*nchunks) = 0;
   end

   x = reshape(x, [3, nchunks]);        % reshape the data
   y = repmat(uint8(0), 4, nchunks);    % for the encoded data

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Split up every 3 bytes into 4 pieces
   %
   %    aaaaaabb bbbbcccc ccdddddd
   %
   % to form
   %
   %    00aaaaaa 00bbbbbb 00cccccc 00dddddd
   %
   y(1,:) = bitshift(x(1,:), -2);                  % 6 highest bits of x(1,:)

   y(2,:) = bitshift(bitand(x(1,:), 3), 4);        % 2 lowest bits of x(1,:)
   y(2,:) = bitor(y(2,:), bitshift(x(2,:), -4));   % 4 highest bits of x(2,:)

   y(3,:) = bitshift(bitand(x(2,:), 15), 2);       % 4 lowest bits of x(2,:)
   y(3,:) = bitor(y(3,:), bitshift(x(3,:), -6));   % 2 highest bits of x(3,:)

   y(4,:) = bitand(x(3,:), 63);                    % 6 lowest bits of x(3,:)

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Now perform the following mapping
   %
   %   0  - 25  ->  A-Z
   %   26 - 51  ->  a-z
   %   52 - 61  ->  0-9
   %   62       ->  +
   %   63       ->  /
   %
   % We could use a mapping vector like
   %
   %   ['A':'Z', 'a':'z', '0':'9', '+/']
   %
   % but that would require an index vector of class double.
   %
   z = repmat(uint8(0), size(y));
   i =           y <= 25;  z(i) = 'A'      + double(y(i));
   i = 26 <= y & y <= 51;  z(i) = 'a' - 26 + double(y(i));
   i = 52 <= y & y <= 61;  z(i) = '0' - 52 + double(y(i));
   i =           y == 62;  z(i) = '+';
   i =           y == 63;  z(i) = '/';
   y = z;

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Add padding if necessary.
   %
   npbytes = 3 * nchunks - ndbytes;     % number of padding bytes
   if npbytes
      y(end-npbytes+1 : end) = '=';     % '=' is used for padding
   end

   if isempty(eol)

      % reshape to a row vector
      y = reshape(y, [1, nebytes]);

   else

      nlines = ceil(nebytes / 76);      % number of lines
      neolbytes = length(eol);          % number of bytes in eol string

      % pad data so it becomes a multiple of 76 elements
      y = [y(:) ; zeros(76 * nlines - numel(y), 1)];
      y(nebytes + 1 : 76 * nlines) = 0;
      y = reshape(y, 76, nlines);

      % insert eol strings
      eol = eol(:);
      y(end + 1 : end + neolbytes, :) = eol(:, ones(1, nlines));

      % remove padding, but keep the last eol string
      m = nebytes + neolbytes * (nlines - 1);
      n = (76+neolbytes)*nlines - neolbytes;
      y(m+1 : n) = '';

      % extract and reshape to row vector
      y = reshape(y, 1, m+neolbytes);

   end

   % output is a character array
   y = char(y);

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [L,S] = low_rank_matrix_completion(params)
% Given partial observation of a low rank matrix, attempts to complete it.
%
% function [L,S] = low_rank_matrix_completion(A)
%
% This example demonstrates how to use the geometry factory for the
% embedded submanifold of fixed-rank matrices, fixedrankembeddedfactory.
% This geometry is described in the paper
% "Low-rank matrix completion by Riemannian optimization"
% Bart Vandereycken - SIAM Journal on Optimization, 2013.
%
% This can be a starting point for many optimization problems of the form:
%
% minimize f(X) such that rank(X) = k, size(X) = [m, n].
%
% Note that the code is long because it showcases quite a few features of
% Manopt: most of the code is optional.
%
% Input:  None. This example file generate random data.
%
% Output: None.

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Nicolas Boumal, July 15, 2014
% Contributors:
%
% Change log:
%
global P;
global PA;
global problem;

k = 1; % rank
A = params.A;
P = params.Omega;
[m,n] = size(A);
%P = sparse(randi([0 1],m,n));

% Random data generation. First, choose the size of the problem.
% We will complete a matrix of size mxn of rank k:
%m = 200;
%n = 500;
%k = 10;
% Generate a random mxn matrix A of rank k
%L = randn(m, k);
%R = randn(n, k);
%A = L*R';
% Generate a random mask for observed entries: P(i, j) = 1 if the entry
% (i, j) of A is observed, and 0 otherwise.
%fraction = 4 * k*(m+n-k)/(m*n);
%P = sparse(rand(m, n) <= fraction);
% Hence, we know the nonzero entries in PA:
PA = P.*A;

% Pick the manifold of matrices of size mxn of fixed rank k.
problem.M = fixedrankembeddedfactory(m, n, k);

% Define the problem cost function. The input X is a structure with
% fields U, S, V representing a rank k matrix as U*S*V'.
% f(X) = 1/2 * || P.*(X-A) ||^2
problem.cost = @cost;

% Define the Euclidean gradient of the cost function, that is, the
% gradient of f(X) seen as a standard function of X.
% nabla f(X) = P.*(X-A)
problem.egrad = @egrad;

% This is optional, but it's nice if you have it.
% Define the Euclidean Hessian of the cost at X, along H, where H is
% represented as a tangent vector: a structure with fields Up, Vp, M.
% This is the directional derivative of nabla f(X) at X along Xdot:
% nabla^2 f(X)[Xdot] = P.*Xdot
problem.ehess = @euclidean_hessian;

% Check consistency of the gradient and the Hessian. Useful if you
% adapt this example for a new cost function and you would like to make
% sure there is no mistake.
% warning('off', 'manopt:fixedrankembeddedfactory:exp');
% checkgradient(problem); pause;
% checkhessian(problem); pause;

% Compute an initial guess. Points on the manifold are represented as
% structures with three fields: U, S and V. U and V need to be
% orthonormal, S needs to be diagonal.
[U, S, V] = svds(PA, k);
X0.U = U;
X0.S = S;
X0.V = V;

% Minimize the cost function using Riemannian trust-regions, starting
% from the initial guess X0.
X = trustregions(problem, X0);

% The reconstructed matrix is X, represented as a structure with fields
% U, S and V.
L = X.U*X.S*X.V'; % low-rank
S = A - L; % sparse

fprintf('||X-A||_F = %g\n', norm(L - A, 'fro'));

% Alternatively, we could decide to use a solver such as
% steepestdescent or conjugategradient. These solvers need to solve a
% line-search problem at each iteration. Standard line searches in
% Manopt have generic purpose systems to do this. But for the problem
% at hand, it so happens that we can rather accurately guess how far
% the line-search should look, and it would be a waste to not use that.
% Look up the paper referenced above for the mathematical explanation
% of the code below.
%
% To tell Manopt about this special information, we specify the
% linesearch hint function in the problem structure. Notice that this
% is not the same thing as specifying a linesearch function in the
% options structure.
%
% Both the SD and the CG solvers will detect that we
% specify the hint function below, and they will use an appropriate
% linesearch algorithm by default, as a result. Typically, they will
% try the step t*H first, then if it does not satisfy an Armijo
% criterion, they will decrease t geometrically until satisfaction or
% failure.
%
% Just like the cost, egrad and ehess functions, the linesearch
% function could use a store structure if you like. The present code
% does not use the store structure, which means quite a bit of the
% computations are made redundantly, and as a result a better method
% could appear slower. See the Manopt tutorial about caching when you
% are ready to switch from a proof-of-concept code to an efficient
% code.
%
% The inputs are X (a point on the manifold) and H, a tangent vector at
% X that is assumed to be a descent direction. That is, there exists a
% positive t such that f(Retraction_X(tH)) < f(X). The function below
% is supposed to output a "t" that it is a good "guess" at such a t.
% problem.linesearch = @linesearch_helper;

% Notice that for this solver, the Hessian is not needed.
% [Xcg, xcost, info, options] = conjugategradient(problem, X0);
%
% fprintf('Take a look at the options that CG used:\n');
% disp(options);
% fprintf('And see how many trials were made at each line search call:\n');
% info_ls = [info.linesearch];
% disp([info_ls.costevals]);
%
% fprintf('Try it again without the linesearch helper.\n');

% Remove the linesearch helper from the problem structure.
% problem = rmfield(problem, 'linesearch');
%
% [Xcg, xcost, info, options] = conjugategradient(problem, X0);
%
% fprintf('Take a look at the options that CG used:\n');
% disp(options);
% fprintf('And see how many trials were made at each line search call:\n');
% info_ls = [info.linesearch];
% disp([info_ls.costevals]);

% If the problem has a small enough dimension, we may (for analysis
% purposes) compute the spectrum of the Hessian at a point X. This may
% help in studying the conditioning of a problem. If you don't provide
% the Hessian, Manopt will approximate the Hessian with finite
% differences of the gradient and try to estimate its "spectrum" (it's
% not a proper linear operator). This can give some intuition, but
% should not be relied upon.
% if problem.M.dim() < 100
%   fprintf('Computing the spectrum of the Hessian...');
%   s = hessianspectrum(problem, X);
%   hist(s);
% end

end

author: andrewssobral

repo: lrslibrary


text:
function help_mex

d = dir('*.mexsg64');		% ".dll'

for i = 1:length(d)
	n = d(i).name;
	[path,name,ext] = fileparts(n);
	myeval(name, 'disp(lasterr)');
end

function myeval(s1, s2)

disp(s1)
eval(s1, s2)
disp('')

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [ft3] = mi32ft3(mi3)
% Convert volume from cubic miles to cubic feet.
% Chad Greene 2012
ft3 = mi3*147197952000;

author: Sable

repo: mcbench-benchmarks


text:
function [TorF, vstr, rdate] = have_feature_ktrlink()
%HAVE_FEATURE_KTRLINK  Detect availability/version info for KTRLINK
%
%   Feature detection function implementing 'ktrlink' tag for HAVE_FEATURE
%   to detect availability/version of Artelys Knitro prior to version 9.0.0,
%   which required the MATLAB Optimization Toolbox.
%
%   See also HAVE_FEATURE, HAVE_FEATURE_KNITRO, KTRLINK.

%   MP-Opt-Model
%   Copyright (c) 2004-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%% ktrlink for pre-Knitro 9.0, requires Optim Toolbox
TorF = exist('ktrlink', 'file');
vstr = '';
rdate = '';
if TorF
    try
        str = evalc(['[x fval] = ktrlink(@(x)1,1);']);
    end
    TorF = exist('fval', 'var') && fval == 1;
    if TorF
        pat = 'KNITRO ([^\s]+)\n|Knitro ([^\s]+)\n';
        [s,e,tE,m,t] = regexp(str, pat);
        if ~isempty(t)
            vstr = t{1}{1};
        end
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function f = cost(X)
  global P;
  global PA;
  % Note that it is very much inefficient to explicitly construct the
  % matrix X in this way. Seen as we only need to know the entries
  % of Xmat corresponding to the mask P, it would be far more
  % efficient to compute those only.
  Xmat = X.U*X.S*X.V';
  f = .5*norm( P.*Xmat - PA , 'fro')^2;
end


text:
function [Tx,Ty,Tz,nr,A] = surfcv(x,y,z,f,const);
%SURFCV 3-D surface of constant value plot.
%   SURFCV(X,Y,Z,F,CONST) draws a surface of constant value
%   for a function of three variables:  F(X,Y,Z) = CONST .
%   The arrays X,Y,Z define the coordinates for F and must
%   be monotonic and 3-D plaid (as if produced by meshgrid).
%   The surface is calculated by linear interpolation of the
%   values of F and is drawn with triangular patches using
%   FILL3.  The color of each patch is determined by its
%   orientation with respect to a specified direction.
%   F must be an M-by-N-by-P volume array.
%
%   SURFCV(F,CONST) assumes X=1:N, Y=1:M, Z=1:P.
%
%   H = SURFCV(...) returns a column vector of handles
%   to PATCH objects, one handle per patch.
%
%   [TX,TY,TZ] = SURFCV(X,Y,Z,F,CONST); returns 3-by-K
%   arrays specifying x,y,z coordinates of the vertices
%   of the triangular patches, where K is the number
%   of patches.  FILL3(TX,TY,TZ,C) can be used to draw
%   the surface.
%
%   [TX,TY,TZ,NR] = SURFCV(...) also returns 3-by-K
%   arrays of normal vectors for each patch which
%   can be used to compute lighting conditions of the
%   surface.  Vectors are normalized to unity and
%   point in the direction of increasing F.
%
%   [TX,TY,TZ,NR,A] = SURFCV(...) returns the area
%   of each patch.
%
%   Example 1: To visualize the surface
%              1./(X.^2+Y.^2) +  1./(Y.^2+Z.^2) + 1./(Z.^2+X.^2) = 2
%              over the range -3 < X < 3, -3 < Y < 3, -3 < Z < 3:
%
%      [x,y,z] = meshgrid(-3:.4:3, -3:.4:3, -3:.4:3);
%      f = 1./(x.^2+y.^2) +  1./(y.^2+z.^2) + 1./(z.^2+x.^2);
%      const = 2;
%      surfcv(x,y,z,f,const);
%
%   Example 2: The routine does not require the grid to be Cartesian
%              and can handle a more general problem of the form:
%                 F(U(X,Y,Z),V(X,Y,Z),W(X,Y,Z)) = CONST
%              This allows to change surface appearance not only
%              by modifying the function F, but also by deforming
%              the grid:
%
%      [x,y,z] = meshgrid(-3:.4:3, -3:.4:3, -3:.4:3);
%      f = 1./(x.^2+y.^2) +  1./(y.^2+z.^2) + 1./(z.^2+x.^2);
%      const = 2;
%      z = z + 0.2.*(x.^2 - y.^2);      %  Grid deformation
%      surfcv(x,y,z,f,const);
%

%  Copyright (c) 1997 Ruslan L. Davidchack
%                     University of Kansas, Lawrence, KS
%                     e-mail: ruslan@ukans.edu
%                     URL:    http://www.ukans.edu/home/ruslan
%      created:  Sep 29, 1997;
%      modified: Dec 23, 1998;

%Short description of the algorithm.
%   Function F is defined on the M-by-N-by-P 3-D grid which
%   divides the space into (M-1)*(N-1)*(P-1) cubes.
%   F specifies values of the function at the vertices.
%   The surface of constant value cuts a patch through a
%   cube if F < CONST for some vertices of the cube and
%   F > CONST for others.
%     Note: F is slightly modified on the input to assure
%           that it does not exactly equal CONST on a vertex,
%           and thus a vertex of a surface element belongs to
%           only one edge (see below).
%   Every vertex of the patch lays on the edge of the
%   cube and its coordinates can be determined by linear
%   interpolation between the values of the function on the
%   ends of the edge. Two vertices of a patch are connected
%   if they lay on the edges that belong to the same face of
%   a cube. If a patch has more than three vertices, it is
%   subdivided into triangular elements by additional
%   connections between vertices that do not belong to the
%   same face.
%   The idea of the algorithm is to group the edges which
%   contain vertices of the surface patches into pairs according
%   to their affiliation with the faces (there can be only 2 or 4
%   patch vertices on a face).  Then the edge pairs are grouped
%   according to their affiliation with the cubes.  The "cube-edge"
%   array, thus created, is then decomposed into triplets of patch
%   vertices.  There can be more than one patch per cube and
%   the algorithm handles this case as well.
%Naming convensions
%   ve - vertex, ed - edge, fa - face, cu - cube
%   nv - number of vertices, ne - number of edges, etc.
%   ce - coordinates of edges, cf - coordinates of faces, etc.
%   ie - indices of edges, ia - indices of faces, etc.

%Simple case to follow the algorithm
%      [x,y,z] = meshgrid(1:4, 1:3, 1:2); const = 15;
%      f = x.^2 + y.^2 + z.^2;
%Different cool surfaces:
% 1.   [x,y,z] = meshgrid(-3:.2:3, -3:.2:3, -3:.2:1.6); const = 7;
%      f = 1./(x.^2+y.^2)+1./(y.^2+z.^2)+1./(z.^2+x.^2)+x.^2+y.^2+z.^2;
%
% 2.   [x,y,z] = meshgrid(-2:.2:2, -2:.2:2, -2:.2:2);  const = 0;
%      f = (x.^2 + y.^2 + z.^2).^2 - 4*(x.^2 - y.^2 - z.^2);
%
% 3.   [x,y,z] = meshgrid(-2:.2:2, -2:.2:2, -2:.2:2);  const = 0;
%      f = x.^2 + y.^2 - z.^2.*(2 - z)./(2 + z);

  if nargin == 2,
    f = x; const = y;
    [x,y,z] = meshgrid(1:size(f,2),1:size(f,1),1:size(f,3));
  elseif nargin ~=5
    disp('Wrong number of input arguments'); return;
  end,

%Make sure F == CONST is never true on the grid
  f1 = f; rngf = max(max(max(f))) - min(min(min(f)));
  if ~isfinite(rngf), rngf = 1; end,
  f1(find(f == const)) = const + 1e-12*rngf;

%Number of vertices in each direction:  nv = [M N P]
%   Total number of vertices = prod(nv)
  nv = size(f1);  iv = reshape(1:prod(nv),nv);
  [ix,iy,iz] = meshgrid(1:nv(2),1:nv(1),1:nv(3));

%Number of edges of each orientation:
%			ne = [(M-1)*N*P  M*(N-1)*P  M*N*(P-1)]
%   Total number of edges = sum(ne)
  ne = prod(nv) - prod(nv([2 1 1;3 3 2]));

%Determine indices of vertices that belong to each edge
%   Indexing of vertices is the same as in f(:)
  ed = zeros(sum(ne),2);
%   Edges parallel to y-axis
  ed((1:ne(1)),1) = reshape(iv(1:nv(1)-1,:,:),ne(1),1);
  ed((1:ne(1)),2) = ed((1:ne(1)),1) + 1;
%   Edges parallel to x-axis
  ie = (1:ne(2)) + ne(1);
  ed(ie,1) = reshape(iv(:,1:nv(2)-1,:),ne(2),1);
  ed(ie,2) = ed(ie,1) + nv(1);
%   Edges parallel to z-axis
  ie = (1:ne(3)) + ne(1) + ne(2);
  ed(ie,1) = (1:ne(3))';
  ed(ie,2) = ed(ie,1) + nv(1)*nv(2);

%Select edges crossed by the surface
  se = find(prod(f1(ed)' - const) < 0)';
  if isempty(se),
    disp('There is no surface in the specified range');
    if nargout > 0, Tx = []; end;
    if nargout > 2, Ty = []; Tz = []; end;
    if nargout > 3, nr = []; end;
    if nargout > 4, A = []; end;
    return; end;
  ed = ed(se,:);

%Restore coordinates of the selected edges
  ce = zeros(length(se),4);
  ce(:,1:3) = [ix(ed(:,1)) iy(ed(:,1)) iz(ed(:,1))];
  ce(:,4) = (se > ne(1)) + (se > ne(1)+ne(2)) + 1;

%Determine coordinates of the patch vertices
  dd = diff(f1(ed)')';  i1 = find(dd < 0);
  dd = (const - f1(ed(:,1)))./dd;
  np = zeros(3,length(se));  np(1,:) = diff(x(ed)');
  np(2,:) = diff(y(ed)');    np(3,:) = diff(z(ed)');
  xv = np(1,:)'.*dd + x(ed(:,1));
  yv = np(2,:)'.*dd + y(ed(:,1));
  zv = np(3,:)'.*dd + z(ed(:,1));
  np(:,i1) = -np(:,i1);

  ie1 = find(ce(:,4) == 1); ie2 = find(ce(:,4) == 2);
  ie3 = find(ce(:,4) == 3);
  ie = cumsum([length(ie1); length(ie2); length(ie3)]);

%Find faces containing the selected edges (maximum of 4 faces per edge)
  cf = repmat([(1:length(se))' ce],1,4);
%Edge-to-face coordinate transformation matrix
  etof = [0 0 0 0 0   0 -1  0  0 0   0 0 0 0  2   0  0  0 -1  2;...
	  0 0 0 0 0   0  0  0 -1 0   0 0 0 0 -1   0  0 -1  0 -1;...
          0 0 0 0 0   0  0 -1  0 0   0 0 0 0 -1   0 -1  0  0 -1];
  ie = diff([0;ie]);
  cf = cf + [repmat(etof(1,:),ie(1),1);repmat(etof(2,:),ie(2),1);...
             repmat(etof(3,:),ie(3),1)];
  cf = reshape(cf',5,4*length(se));

%Discard faces that are outside the range of x,y,z
  cf(:,find(prod(cf) == 0)) = [];  cf = cf';
  i1 = find(cf(:,5) ~= 2);  cf(i1(find(cf(i1,3) > nv(1)-1)),:) = [];
  i2 = find(cf(:,5) < 3);   cf(i2(find(cf(i2,2) > nv(2)-1)),:) = [];
  i3 = find(cf(:,5) > 1);   cf(i3(find(cf(i3,4) > nv(3)-1)),:) = [];

%Find indices of the selected faces
  i1 = find(cf(:,5)==1); i2 = find(cf(:,5)==2); i3 = find(cf(:,5)==3);
%Number of faces of each orientation:
%                ne = [(M-1)*(N-1)*P  M*(N-1)*(P-1)  (M-1)*N*(P-1)]
  nf = [nv(3)*(nv(2)-1)*(nv(1)-1) nv(1)*(nv(2)-1)*(nv(3)-1)...
        nv(2)*(nv(3)-1)*(nv(1)-1)];
  fa = zeros(size(cf,1),1);
  fa(i1) = cf(i1,3)+((cf(i1,2)-1)+(cf(i1,4)-1)*(nv(2)-1))*(nv(1)-1);
  fa(i2) = cf(i2,3)+((cf(i2,2)-1)+(cf(i2,4)-1)*(nv(2)-1))*nv(1)+nf(1);
  fa(i3) = cf(i3,3)+((cf(i3,2)-1)+(cf(i3,4)-1)*nv(2))*(nv(1)-1)+...
	   nf(1) + nf(2);
  [fa ia] = sort(fa);

%!!!!!!!! Check whether fac values come in pairs.
%   Each face can have only 2 or 4 selected edges.
%   The check below can be removed after debugging process
%   is completed
if mod(length(fa),2) ~= 0,
  disp('Odd number!'); return,
else
  if sum(diff(reshape(fa,2,length(fa)/2))) > 0,
  disp('Odd number!'); return, end, end;
%!!!!!!!!

%Combine face coord. with indices of edges that belong to the faces
  cf = [fa(1:2:end) reshape(cf(ia,1),2,length(fa)/2)' ...
        cf(ia(1:2:end),2:5)];

%Find cubes containing selected faces
  i1 = find(cf(:,7)==1); i2 = find(cf(:,7)==2); i3 = find(cf(:,7)==3);
  cc = repmat(cf(:,2:6),1,2);
  ftoc = [zeros(3,7) [0 0 -1;0 -1 0;-1 0 0]];
  cc = cc + [repmat(ftoc(1,:),length(i1),1);...
   repmat(ftoc(2,:),length(i2),1);repmat(ftoc(3,:),length(i3),1)];
  cc = reshape(cc',5,size(cf,1)*2);
  cc(:,find(prod(cc) == 0 | cc(4,:) > nv(1)-1 | cc(3,:) > nv(2)-1 |...
            cc(5,:) > nv(3)-1)) = [];  cc = cc';

  cu = cc(:,4) + ((cc(:,3)-1) + (cc(:,5)-1)*(nv(2)-1))*(nv(1)-1);
  [cu ic] = sort(cu);
%Combine indices of selected cubes with indices of edges that belong
%   to them into a cube-edge index array
  cc = [cu cc(ic,1:2)];

%Combine edges into triplets composing surface patches
  ip = [];   cc = sortrows(cc);
  while 1,
%Indices of patches
    ic = find(diff([0;cc(:,1);0]) ~= 0);
%Indices of triangular patches
    ic3 = ic(find(diff(ic) == 3));
%Indices of patches with more than 3 vertices
    icm = ic(find(diff(ic) > 3));

%Select triangular patches and take them out of the cube-edge array
    ic(end) = [];   ip = [ip [cc(ic,2:3)';cc(ic+1,3)']];
    if ~isempty(ic3),
      irm = repmat(ic3,1,3) + repmat(0:2,length(ic3),1);
    else irm = []; end,
    irm = [irm(:);icm];   cc(icm+1,2) = cc(icm,3);   cc(irm,:) = [];

%Sort new cube-edge index array or quit if it is empty
    if ~isempty(cc),  cc = sortrows(cc);  else,
      disp(['Number of patches = ' num2str(size(ip,2))]);
      break; end

%Find cases when there is more than one patch per cube
    irm = find(sum(abs(diff(cc)')) == 0);
    if ~isempty(irm),
      %disp(['More than one patch in ' int2str(length(irm)) ' cubes']);
      irm = repmat(irm,2,1)+repmat([0;1],1,length(irm)); irm = irm(:);
      cc(irm,:) = []; end,

  end;

  Tx = xv(ip);  Ty = yv(ip);  Tz = zv(ip);

%This is to compute area of each patch and the normal vector
  x2 = Tx(2,:)-Tx(1,:); y2 = Ty(2,:)-Ty(1,:); z2 = Tz(2,:)-Tz(1,:);
  x3 = Tx(3,:)-Tx(1,:); y3 = Ty(3,:)-Ty(1,:); z3 = Tz(3,:)-Tz(1,:);
  if nargout == 5
     A = sqrt((x2.*y3-x3.*y2).^2 + (y2.*z3-y3.*z2).^2 +...
              (z2.*x3-z3.*x2).^2)./2; end
  nr = [y2.*z3-y3.*z2; z2.*x3-z3.*x2; x2.*y3-x3.*y2];
  nr = nr./repmat(sqrt(sum(nr.^2)),3,1);
  i1 = find(sum(nr.*np(:,ip(1,:))) < 0);
  nr(:,i1) = -nr(:,i1);
  i2 = ip(2,i1); ip(2,i1) = ip(3,i1); ip(3,i1) = i2;
  Tx = xv(ip);  Ty = yv(ip);  Tz = zv(ip);
  if nargout < 3,
    va = [-2 3 1];
% Colored according to the patch orientation
    C = va*nr;
% Colored according to the patch distance from the origin
%    C = sqrt(Tx.^2 + Ty.^2 + Tz.^2);
    clf;  H = fill3(Tx,Ty,Tz,C);  axis image;
    set(gca,'xlim',[min(min(Tx)) max(max(Tx))],...
            'ylim',[min(min(Ty)) max(max(Ty))],...
            'zlim',[min(min(Tz)) max(max(Tz))],'box','on');
    view(va);
    Tx = H;
  end
  if nargout == 0, clear Tx, end
  if nargout < 4, clear nr, end


text:
%% computeKappas.m
% This scripts computes normalization coefficient for each image in database.
%
% 02-08-11 Michal Uricar

clc;
close all; clearvars;

%% Timestamp

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

%% load GT - TRN

gt_fname = './MAT/GT_nImages-all.mat';
load(gt_fname);

%% compute kappas

kappa = inf(1, numel(gt));

for i = 1 : numel(gt)
    GT = gt{i};
    C = ( (GT(:, 2)+GT(:, 6))/2 + (GT(:, 3)+GT(:, 7))/2 )/2;
    kappa(i) = 1/sqrt(sum(( C - (GT(:, 4)+GT(:, 5))/2 ).^2));
end;

%% save kappas

save('./MAT/kappas.mat', 'kappa');

%% load GT - VAL

gt_fname = './MAT/GTVal_nImages-all.mat';
load(gt_fname);

%% compute kappas

kappa = inf(1, numel(gt));

for i = 1 : numel(gt)
    GT = gt{i};
    C = ( (GT(:, 2)+GT(:, 6))/2 + (GT(:, 3)+GT(:, 7))/2 )/2;
    kappa(i) = 1/sqrt(sum(( C - (GT(:, 4)+GT(:, 5))/2 ).^2));
end;

%% save kappas

save('./MAT/kappasVAL.mat', 'kappa');

%% load GT - TST

gt_fname = './MAT/GTTst_nImages-all.mat';
load(gt_fname);

%% compute kappas

kappa = inf(1, numel(gt));

for i = 1 : numel(gt)
    GT = gt{i};
    C = ( (GT(:, 2)+GT(:, 6))/2 + (GT(:, 3)+GT(:, 7))/2 )/2;
    kappa(i) = 1/sqrt(sum(( C - (GT(:, 4)+GT(:, 5))/2 ).^2));
end;

%% save kappas

save('./MAT/kappasTST.mat', 'kappa');

%% Timestamp

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

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function nmt_peaksearch_helper
% simplifies Callback for peak search functionality in GUI
global st

cfg=[];
cfg.searchradius = [str2num(get(st.nmt.gui.searchradius1,'String')) str2num(get(st.nmt.gui.searchradius2,'String'))];

cfg.peaktype=get(st.nmt.gui.peaktype,'string');
cfg.peaktype=cfg.peaktype{get(st.nmt.gui.peaktype,'Value')};

peakdomain=get(st.nmt.gui.peakdomain,'string');
peakdomain=peakdomain{get(st.nmt.gui.peakdomain,'Value')};

switch(peakdomain)
    case 'spatial'
        cfg.time = st.nmt.cfg.time_idx;
    case 'temporal'
        cfg.vox = st.nmt.cfg.vox_idx;
    case 'spatiotemporal'
        % nothing to do, this is default behavior
    otherwise
        error('well this is unexpected...')
end

[v,t]=nmt_peaksearch(cfg);

nmt_repos(v,t);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [MHz] = GHz2MHz(GHz)
% Convert frequency from gigahertz to megahertz.
% Chad A. Greene 2012
MHz = GHz*1e+3;


text:
function k = findcols(A, b)
%FINDCOLS Find indices of a given column within a matrix.
%
%   FINDCOLS(A, B) returns a row vector with the indices of the columns
%   in the matrix A that are identical to the column vector B.  If no
%   columns in A are identical to B, an empty vector is returned.
%
%   The methods uses a for-loop, but it uses less memory and is in many
%   cases a lot faster than the vectorized methods
%
%      find( all( A == repmat(b, 1, size(A, 2)), 1 ) )
%      find( all( A == b(:,ones(size(A, 2), 1)), 1 ) )
%
%   See also FIND, FINDROWS.

%   Author:      Peter J. Acklam
%   Time-stamp:  2002-03-03 13:51:19 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   k = find( A(1,:) == b(1) );
   for j = 2:size(A, 1)
      k = k( A(j,k) == b(j) );
      if isempty(k)
         return
      end
   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
% this code is revised based on ILSVRC 2013 (http://www.image-net.org/challenges/LSVRC/2013/)

function  zeroShot = zeroShot_top_recall_Phrase(Nre, tuple_confs_cell, tuple_labels_cell, sub_bboxes_cell, obj_bboxes_cell)

load('gt.mat','gt_tuple_label','gt_obj_bboxes','gt_sub_bboxes');
load('zeroShot.mat','zeroShot');

for ii = 1 : 1000
    det = find(zeroShot{ii} == 0);
    gt_tuple_label{ii}(det,:) = [];
    gt_obj_bboxes{ii}(det,:) = [];
    gt_sub_bboxes{ii}(det,:) = [];
end

%num_imgs = length(gt_tuple_label);
num_imgs = 1000;
for i=1:num_imgs
    [tuple_confs_cell{i}, ind] = sort(tuple_confs_cell{i},'descend');
    if length(ind) >= Nre
        tuple_confs_cell{i} = tuple_confs_cell{i}(1:Nre);
        tuple_labels_cell{i} = tuple_labels_cell{i}(ind(1:Nre),:);
        obj_bboxes_cell{i} = obj_bboxes_cell{i}(ind(1:Nre),:);
        sub_bboxes_cell{i} = sub_bboxes_cell{i}(ind(1:Nre),:);
    else
        tuple_labels_cell{i} = tuple_labels_cell{i}(ind,:);
        obj_bboxes_cell{i} = obj_bboxes_cell{i}(ind,:);
        sub_bboxes_cell{i} = sub_bboxes_cell{i}(ind,:);
    end
end

num_pos_tuple = 0;
for ii = 1 : num_imgs
    num_pos_tuple = num_pos_tuple + size(gt_tuple_label{ii},1);
end

tp_cell = cell(1,num_imgs);
fp_cell = cell(1,num_imgs);

gt_thr = 0.5;
% iterate over images
for i=1:num_imgs

    gt_tupLabel = gt_tuple_label{i};
    if ~isempty(gt_obj_bboxes{i})
        gt_box_entity = [min(gt_obj_bboxes{i}(:,1:2),gt_sub_bboxes{i}(:,1:2)),max(gt_obj_bboxes{i}(:,3:4),gt_sub_bboxes{i}(:,3:4))];
    else
        gt_box_entity = [];
    end

    num_gt_tuple = size(gt_tupLabel,1);
    gt_detected = zeros(1,num_gt_tuple);

    labels = tuple_labels_cell{i};
    boxObj = obj_bboxes_cell{i};
    boxSub = sub_bboxes_cell{i};
    if ~isempty(boxObj)
        box_entity_our  = [min(boxObj(:,1:2), boxSub(:,1:2)), max(boxObj(:,3:4), boxSub(:,3:4))];
    else
        box_entity_our  = [];
    end

    num_obj = size(labels,1);
    tp = zeros(1,num_obj);
    fp = zeros(1,num_obj);
    for j=1:num_obj

        bbO = box_entity_our(j,:);
        ovmax = -inf;
        kmax = -1;

        for k=1:num_gt_tuple
            if norm(labels(j,:) - gt_tupLabel(k,:),2) ~= 0
                continue;
            end
            if gt_detected(k) > 0
                continue;
            end

            bbgtO = gt_box_entity(k,:);

            biO=[max(bbO(1),bbgtO(1)) ; max(bbO(2),bbgtO(2)) ; min(bbO(3),bbgtO(3)) ; min(bbO(4),bbgtO(4))];
            iwO=biO(3)-biO(1)+1;
            ihO=biO(4)-biO(2)+1;

            if iwO>0 & ihO>0
                % compute overlap as area of intersection / area of union
                uaO=(bbO(3)-bbO(1)+1)*(bbO(4)-bbO(2)+1)+...
                   (bbgtO(3)-bbgtO(1)+1)*(bbgtO(4)-bbgtO(2)+1)-...
                   iwO*ihO;
                ov =iwO*ihO/uaO;

                % makes sure that this object is detected according
                % to its individual threshold
                if ov >= gt_thr && ov > ovmax
                    ovmax=ov;
                    kmax=k;
                end
            end
        end

        if kmax > 0
            tp(j) = 1;
            gt_detected(kmax) = 1;
        else
            fp(j) = 1;
        end
    end

    % put back into global vector
    tp_cell{i} = tp;
    fp_cell{i} = fp;

end

t = tic;
tp_all = [];
fp_all = [];
confs = [];
for ii = 1 : num_imgs
tp_all = [tp_all; tp_cell{ii}(:) ];
fp_all = [fp_all; fp_cell{ii}(:) ];
confs = [confs; tuple_confs_cell{ii}(:)];
end

[confs, ind] = sort(confs,'descend');
tp_all = tp_all(ind);
fp_all = fp_all(ind);

tp = cumsum(tp_all );
fp = cumsum(fp_all );
recall =(tp/num_pos_tuple);
zeroShot = recall(end);

end

author: Prof-Lu-Cewu

repo: Visual-Relationship-Detection


text:
% DROSPLOTCHIPDISTANCES Plot the accuracy figures appearing in the paper
% FORMAT
% DESC Plot the accuracy figures appearing in the paper
%
% COPYRIGHT : Antti Honkela, 2009

% SHEFFIELDML

tfnames = drosTF.names;
FONTSIZE = 8;
styles = {'bo-', 'rd--', 'm*--', 'gs--', 'k--'};
t = [20, 100, 250];

rankings = {};
for k=1:length(tfnames),
  tf = tfnames{k};
  nonmutarankings{k} = {indrank.(tf), [], corrrank.(tf)};
  if isfield(mutarank, tf),
    rankings{k} = {indrank.(tf), mutarank.(tf), corrrank.(tf)};
  else
    rankings{k} = nonmutarankings{k};
  end
end

clear accs;
clear pvals;
figure(1);
for k=1:length(tfnames),
  subplot(3, length(tfnames), k); %3*(k-1)+2);
  tf = tfnames{k};

  [accs(:, :, k), pvals(:, :, k)] = drosPlotAccuracyBars(...
      rankings{k}, ...
      chip_validation.(tf), t, styles, [], drosexp, drosinsitu);
  set(gca, 'FontSize', FONTSIZE);
  title(sprintf('Global ChIP: %s', tfnames{k}));
  xlabel('Top N to consider');
  if k==1,
    ylabel('Relative enrichment (%)');
  end
end

clear accs;
clear pvals;
for k=1:length(tfnames),
  subplot(3, length(tfnames), length(tfnames)+k);
  tf = tfnames{k};

  if isfield(mutant_validation, tf),
    [accs(:, :, k), pvals(:, :, k)] = drosPlotAccuracyBars(...
	nonmutarankings{k}, ...
	mutant_validation.(tf), t, styles, [], drosexp, drosinsitu);
    set(gca, 'FontSize', FONTSIZE);
    title(sprintf('Global knock-outs: %s', tfnames{k}));
    xlabel('Top N to consider');
    if k==1,
      ylabel('Relative enrichment (%)');
    end
  end
end

subplot(3, length(tfnames), 3*length(tfnames)-[1, 0]);
if exist('INCLUDE_TSNI') && INCLUDE_TSNI,
  bar(rand(5));
else
  bar(rand(length(rankings{1})));
end
hold on
plot(1:2, t(1:2), 'k--');

axis([-10 -9 -10 -9]);
set(gca, 'FontSize', FONTSIZE);
axis off
if exist('INCLUDE_TSNI') && INCLUDE_TSNI,
  legend('Single-target models', 'Multiple-target models', 'TSNI', ...
	 'Knock-outs', 'Correlation', 'Random', 'Location', 'North');
else
  legend('Single-target models', ... %'Multiple-target models', ...
	 'Knock-outs', 'Correlation', 'Random', 'Location', 'North');
end

set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [20 20])
set(gcf, 'PaperPosition', [0 0 17.4 18])
hold off

clear accs;
clear pvals;
figure(2);
for k=1:length(tfnames),
  subplot(3, length(tfnames), k); %3*(k-1)+2);
  tf = tfnames{k};

  [accs(:, :, k), pvals(:, :, k)] = drosPlotAccuracyBars(...
      rankings{k}, ...
      chip_validation.(tf), t, styles, 1, drosexp, drosinsitu);
  set(gca, 'FontSize', FONTSIZE);
  title(sprintf('Focused ChIP: %s', tfnames{k}));
  xlabel('Top N to consider');
  if k==1,
    ylabel('Relative enrichment (%)');
  end
end

clear accs;
clear pvals;
for k=1:length(tfnames),
  subplot(3, length(tfnames), length(tfnames)+k);
  tf = tfnames{k};

  if isfield(mutant_validation, tf),
    [accs(:, :, k), pvals(:, :, k)] = drosPlotAccuracyBars(...
	nonmutarankings{k}, ...
	mutant_validation.(tf), t, styles, 1, drosexp, drosinsitu);
    set(gca, 'FontSize', FONTSIZE);
    title(sprintf('Focused knock-outs: %s', tfnames{k}));
    xlabel('Top N to consider');
    if k==1,
      ylabel('Relative enrichment (%)');
    end
  end
end

subplot(3, length(tfnames), 3*length(tfnames)-[1, 0]);
if exist('INCLUDE_TSNI') && INCLUDE_TSNI,
  bar(rand(5));
else
  bar(rand(length(rankings{1})));
end
hold on
plot(1:2, 1:2, 'k-.');
plot(1:2, 1:2, 'k--');

axis([-10 -9 -10 -9]);
set(gca, 'FontSize', FONTSIZE);
axis off
if exist('INCLUDE_TSNI') && INCLUDE_TSNI,
  legend('Single-target models', 'Multiple-target models', 'TSNI', ...
	 'Knock-outs', 'Correlation', 'Random', 'Location', 'North');
else
  legend('Single-target models', ... %'Multiple-target models', ...
	 'Knock-outs', 'Correlation', 'Random', 'Location', 'North');
end

set(gcf, 'PaperUnits', 'centimeters');
set(gcf, 'PaperSize', [20 20])
set(gcf, 'PaperPosition', [0 0 17.4 18])
hold off

% figure(3);
% tf = 'twi';

% clear accs;
% clear pvals;

% rankings = {indrank.(tf), ssrank.(tf)};

% subplot(2, 3, 1);
% [accs(:, :, k), pvals(:, :, k)] = drosPlotAccuracyBars(rankings, ...
%     chip_validation.(tf), t, styles, [], drosexp, drosinsitu);
% set(gca, 'FontSize', FONTSIZE);
% title(sprintf('Global ChIP: %s', tf));
% xlabel('Top N to consider');
% ylabel('Relative enrichment (%)');

% subplot(2, 3, 2);
% [accs(:, :, k), pvals(:, :, k)] = drosPlotAccuracyBars(rankings, ...
%     mutant_validation.(tf), t, styles, [], drosexp, drosinsitu);
% set(gca, 'FontSize', FONTSIZE);
% title(sprintf('Global knock-outs: %s', tf));
% xlabel('Top N to consider');
% %ylabel('Relative enrichment (%)');

% subplot(2, 3, 4);
% [accs(:, :, k), pvals(:, :, k)] = drosPlotAccuracyBars(rankings, ...
%     chip_validation.(tf), t, styles, 1, drosexp, drosinsitu);
% set(gca, 'FontSize', FONTSIZE);
% title(sprintf('Focused ChIP: %s', tf));
% xlabel('Top N to consider');
% ylabel('Relative enrichment (%)');

% subplot(2, 3, 5);
% [accs(:, :, k), pvals(:, :, k)] = drosPlotAccuracyBars(rankings, ...
%     mutant_validation.(tf), t, styles, 1, drosexp, drosinsitu);
% set(gca, 'FontSize', FONTSIZE);
% title(sprintf('Focused knock-outs: %s', tf));
% xlabel('Top N to consider');
% %ylabel('Relative enrichment (%)');

% subplot(2, 3, [3,6]);
% bar(rand(2));
% hold on
% plot(1:2, t(1:2), 'k--');
% plot(1:2, 1:2, 'k-.');

% axis([-10 -9 -10 -9]);
% set(gca, 'FontSize', FONTSIZE);
% axis off
% legend('n=12', 'n=7', ...
%        'Random', 'Filtered', 'Location', 'West');

% set(gcf, 'PaperUnits', 'centimeters');
% set(gcf, 'PaperSize', [20 20])
% set(gcf, 'PaperPosition', [0 0 9.65 9.65])
% %set(gcf, 'PaperSize', [20 20])
% %set(gcf, 'PaperPosition', [0 0 8.7 8.7])
% hold off

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [SPP] = point_to_hae_newton(R_TGT_COA, Rdot_TGT_COA, ARP_COA, VARP_COA, SCP, HAE0, delta_MAX, NLIM)
% POINT_TO_HAE_NEWTON transforms pixel row, col to a constant height
% above the ellipsoed via algorithm in SICD Image Projections.
%
% GPP = point_to_hae_newton(R_TGT_COA, Rdot_TGT_COA, ARP_COA, VARP_COA, SCP, HAE0)
%
% Inputs:
%    R_TGT_COA    - range to the ARP at COA
%    Rdot_TGT_COA - range rate relative to the ARP at COA
%    ARP_COA      - aperture reference position at tCOA
%    VARP_COA     - velocity at tCOA
%    SCP          - scene center point (ECF meters)
%    HAE0         - Surface height (m) above the WGS-84 reference ellipsoid
%                   for projection point SPP
%    delta_MAX    - Threshold for convergence of iterative projection
%                   sequence.
%    NLIM         - Maximum number of iterations allowed.
%
% Outputs:
%    SPP          - [3xN] Surface Projection Point position on the HAE0
%                   surface and along the R/Rdot contour
%
% Authors: Rocco Corsetti, NGA/IB
%          Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if ~exist('delta_MAX','var')
    delta_MAX = 1e-3;
end
if ~exist('NLIM','var')
    NLIM = 10;
end
SPP = zeros(3,numel(R_TGT_COA));

%% 9. Precise R/Rdot To Constant HAE Surface Projection

% This way of looping through points really adds no efficiency in MATLAB.
% Really we shouled vectorize the matrix operations and not loop.  However,
% this is an easy way to make its behaviour match the other projection
% functions in the toolbox (point_to_ground_plane, point_to_hae).
for points_i = 1:numel(R_TGT_COA)

    % Rename variables from SICD projection document notation to "Spotlight
    % Synthetic Aperture Radar (SAR) Sensor Model" notation
    R = SCP(:);
    h = HAE0;

    Range_obs = R_TGT_COA(points_i);
    Doppler_obs = -Rdot_TGT_COA(points_i);

    R_ARP_curr = ARP_COA(:,(points_i));
    V_ARP_curr = VARP_COA(:,(points_i));

    R0 = R;
    a = 6378137.0; % WGS84 semimajor axis (m)
    b = 6356752.314245; % WGS84 semiminor axis (m)

    delta = Inf;
    iters = 0;
    while (max(abs(delta)) > delta_MAX) && (iters < NLIM)

        F1 = -9999*ones(3,3);
        F2 = -9999*ones(3,3);
        F3 = -9999*ones(3,3);

        for i = 1:3
            for j = 1:3
                switch j
                    case(1)
                        eps = -0.5;
                    case(2)
                        eps = 0;
                    case(3)
                        eps = +0.5;
                end

                R = R0;
                R(i) = R(i) + eps;
                X = R(1); Y = R(2); Z = R(3);

                r_est = R - R_ARP_curr;

                Range_est = norm(r_est);
                Doppler_est = dot(V_ARP_curr,r_est)/norm(r_est);

                F1(i,j) = Range_obs - Range_est;
                F2(i,j) = Doppler_obs - Doppler_est;
                F3(i,j) = ((X^2 + Y^2)/(a+h)^2) + ((Z^2)/(b+h)^2) - 1;
            end
        end

        dF1dX = (F1(1,3) - F1(1,1)) / (2*0.5);
        dF1dY = (F1(2,3) - F1(2,1)) / (2*0.5);
        dF1dZ = (F1(3,3) - F1(3,1)) / (2*0.5);

        dF2dX = (F2(1,3) - F2(1,1)) / (2*0.5);
        dF2dY = (F2(2,3) - F2(2,1)) / (2*0.5);
        dF2dZ = (F2(3,3) - F2(3,1)) / (2*0.5);

        dF3dX = (F3(1,3) - F3(1,1)) / (2*0.5);
        dF3dY = (F3(2,3) - F3(2,1)) / (2*0.5);
        dF3dZ = (F3(3,3) - F3(3,1)) / (2*0.5);

        f = -[F1(1,2)
            F2(1,2)
            F3(1,2)];

        B = [dF1dX dF1dY dF1dZ
            dF2dX dF2dY dF2dZ
            dF3dX dF3dY dF3dZ];

        delta = B\f;

        R0 = R0 + delta;

        iters = iters + 1;
    end

    SPP(:,points_i) = R0; % Convert back to notation from SICD projection document
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
clear;

%Sensors
%%Acc set1
A1=diag([0.999987092709754 1]);
B1=diag([3.5596949329112e-006 1e-7]);
C1=[1 1];
sR1=(8.55478077790315e-007)^0.5;
sP1=diag([(4.9087e-007)^0.5 1e-3]);
[A1, Q1, R1]=dc2dc_v000(A1,B1*B1',sR1*sR1',1/100,2,1);
B1=Q1^0.5;
sR1=R1^0.5;

A2=A1;A3=A1;A4=A1;A5=A1;
B2=B1;B3=B1;B4=B1;B5=B1;
C2=C1;C3=C1;C4=C1;C5=C1;
sR2=sR1;sR3=sR1;sR4=sR1;sR5=sR1;
sP2=sP1;sP3=sP1;sP4=sP1;sP5=sP1;

%%Acc set2
A11=1;
B11=1e-7;%5.35244311625171e-008;
C11=1;
sR11=(8.55478077790315e-007)^0.5;
sP11=(6.8307e-007*20)^0.5;
[A11, Q11, R11]=dc2dc_v000(A11,B11*B11',sR11*sR11',1/100,2,1);
B11=Q11^0.5;
sR11=R11^0.5;

A12=A11;A13=A11;A14=A11;A15=A11;
B12=B11;B13=B11;B14=B11;B15=B11;
C12=C11;C13=C11;C14=C11;C15=C11;
sR12=sR11;sR13=sR11;sR14=sR11;sR15=sR11;
sP12=sP11;sP13=sP11;sP14=sP11;sP15=sP11;

%gyros
A6=0.999954464990211;
B6=7.50002159071526e-007*10;
C6=1;
sR6=(4.40958675010914e-008)^0.5;
vr_a=markov1st_v000(A6 ,B6,[], 2,100);
sP6=(vr_a(3))^0.5;
[A6, Q6, R6]=dc2dc_v000(A6,B6*B6',sR6*sR6',1/100,2,1);
B6=Q6^0.5;
sR6=R6^0.5;

A7=A6;B7=B6;C7=C6;sR7=sR6;sP7=sP6;

% A7=0.999977286299473;
% B7=2.06907835784507e-006;
% C7=1;
% sR7=(3.74167215805145e-007)^0.5;
% sP7=(9.4241e-008)^0.5;
% [A7, Q7, R7]=dc2dc_v000(A7,B7*B7',sR7*sR7',1/100,2,1);
% B7=Q7^0.5;
% sR7=R7^0.5;

MA1=[cos(0) sin(0) 0;
    cos(0) sin(0) 0;
    cos(-pi/5) sin(-pi/5) 0;
    cos(-pi/2) sin(-pi/2) 0;
    cos(-pi/2) sin(-pi/2) 0];

MG1=[0 0 1;0 0 1];

MA2=[cos(pi/10) sin(pi/10) 0;
    cos(pi/8) sin(pi/8) 0;
    cos(pi/6) sin(pi/6) 0;
    cos(pi/4) sin(pi/4) 0;
    cos(pi/2) sin(pi/2) 0];

%%% Main model
%M=[MA1;MG];
M=[MA1;MG1;MA2];

%imu observation structure
Rimu=sR1*sR1';
Rimu=diagmat_v000(sR2*sR2',Rimu,1);
Rimu=diagmat_v000(sR3*sR3',Rimu,1);
Rimu=diagmat_v000(sR4*sR4',Rimu,1);
Rimu=diagmat_v000(sR5*sR5',Rimu,1);
Rimu=diagmat_v000(sR6*sR6',Rimu,1);
Rimu=diagmat_v000(sR7*sR7',Rimu,1);
Rimu=diagmat_v000(sR11*sR11',Rimu,1);
Rimu=diagmat_v000(sR12*sR12',Rimu,1);
Rimu=diagmat_v000(sR13*sR13',Rimu,1);
Rimu=diagmat_v000(sR14*sR14',Rimu,1);
Rimu=diagmat_v000(sR15*sR15',Rimu,1);

Cimu=C1;
Cimu=diagmat_v000(C2,Cimu,1);
Cimu=diagmat_v000(C3,Cimu,1);
Cimu=diagmat_v000(C4,Cimu,1);
Cimu=diagmat_v000(C5,Cimu,1);
Cimu=diagmat_v000(C6,Cimu,1);
Cimu=diagmat_v000(C7,Cimu,1);
Cimu=diagmat_v000(C11,Cimu,1);
Cimu=diagmat_v000(C12,Cimu,1);
Cimu=diagmat_v000(C13,Cimu,1);
Cimu=diagmat_v000(C14,Cimu,1);
Cimu=diagmat_v000(C15,Cimu,1);

[T Mls]=cp_Tparam_v000(M,Rimu);
Cimu_obs=[T*Cimu];
Rimu_obs=T*Rimu*T';
Rimu_inp=Mls*Rimu*Mls';

A=A1;
A=diagmat_v000(A2,A,1);
A=diagmat_v000(A3,A,1);
A=diagmat_v000(A4,A,1);
A=diagmat_v000(A5,A,1);
A=diagmat_v000(A6,A,1);
A=diagmat_v000(A7,A,1);
A=diagmat_v000(A11,A,1);
A=diagmat_v000(A12,A,1);
A=diagmat_v000(A13,A,1);
A=diagmat_v000(A14,A,1);
A=diagmat_v000(A15,A,1);

Q=B1*B1';
Q=diagmat_v000(B2*B2',Q,1);
Q=diagmat_v000(B3*B3',Q,1);
Q=diagmat_v000(B4*B4',Q,1);
Q=diagmat_v000(B5*B5',Q,1);
Q=diagmat_v000(B6*B6',Q,1);
Q=diagmat_v000(B7*B7',Q,1);
Q=diagmat_v000(B11*B11',Q,1);
Q=diagmat_v000(B12*B12',Q,1);
Q=diagmat_v000(B13*B13',Q,1);
Q=diagmat_v000(B14*B14',Q,1);
Q=diagmat_v000(B15*B15',Q,1);

P=sP1*sP1';
P=diagmat_v000(sP2*sP2',P,1);
P=diagmat_v000(sP3*sP3',P,1);
P=diagmat_v000(sP4*sP4',P,1);
P=diagmat_v000(sP5*sP5',P,1);
P=diagmat_v000(sP6*sP6',P,1);
P=diagmat_v000(sP7*sP7',P,1);
P=diagmat_v000(sP11*sP11',P,1);
P=diagmat_v000(sP12*sP12',P,1);
P=diagmat_v000(sP13*sP13',P,1);
P=diagmat_v000(sP14*sP14',P,1);
P=diagmat_v000(sP15*sP15',P,1);

Pwls=P;

% %%%model for the equivalent systems
%Equivalent models for identical systems
[Mls_e1, Ae1, Be1, Ce1, sRe1, sPe1]=cp_redsys_v000(MA1(:,1:2), A1, B1, C1, sR1, sP1);
[Mls_e2, Ae2, Be2, Ce2, sRe2, sPe2]=cp_redsys_v000(MA2(:,1:2), A11, B11, C11, sR11, sP11);
[Mls_e3, Ae3, Be3, Ce3, sRe3, sPe3]=cp_redsys_v000(MG1(:,3), A6, B6, C6, sR6, sP6);

%imu observation structure
Rimue=sRe1*sRe1';
Rimue=diagmat_v000(sR6*sR6',Rimue,1);
Rimue=diagmat_v000(sR7*sR7',Rimue,1);
Rimue=diagmat_v000(sRe2*sRe2',Rimue,1);

Cimue=Ce1;
Cimue=diagmat_v000(C6,Cimue,1);
Cimue=diagmat_v000(C7,Cimue,1);
Cimue=diagmat_v000(Ce2,Cimue,1);

Me=[1 0 0;0 1 0;MG1;1 0 0;0 1 0];
[Te Mlse]=cp_Tparam_v000(Me,Rimue);
Cimu_obse=[Te*Cimue];
Rimu_obse=Te*Rimue*Te';
Rimu_inpe=Mlse*Rimue*Mlse';

%single system model
Ae=Ae1;
Ae=diagmat_v000(A6,Ae,1);
Ae=diagmat_v000(A7,Ae,1);
Ae=diagmat_v000(Ae2,Ae,1);

Qe=Be1*Be1';
Qe=diagmat_v000(B6*B6',Qe,1);
Qe=diagmat_v000(B7*B7',Qe,1);
Qe=diagmat_v000(Be2*Be2',Qe,1);

Pe=sPe1*sPe1';
Pe=diagmat_v000(sP6*sP6',Pe,1);
Pe=diagmat_v000(sP7*sP7',Pe,1);
Pe=diagmat_v000(sPe2*sPe2',Pe,1);

%%%%%%%%%%%%%%%%%%
nit=60*100;
nst=size(A,1);
nste=size(Ae,1);

x=(P^0.5)*randn(nst,1);
xest=zeros(size(A,1),1);
xeste=zeros(size(Ae,1),1);

sRimu=Rimu^0.5;
nsen=size(sRimu,1);
sQ=Q^0.5;

%debug variables
debctr=0;
debx=zeros(nst,nit);
debp=zeros(nst,nit);
debxe=zeros(nste,nit);
debpe=zeros(nste,nit);
debuopt=zeros(3,nit);
debuopte=zeros(3,nit);
debuwls=zeros(3,nit);
debpopt=zeros(3,nit);
debpopte=zeros(3,nit);
debpwls=zeros(3,nit);
debeq=zeros(5,nit);

% debctr=debctr+1;
% debx(:,debctr)=xest;
% debp(:,debctr)=diag(P).^0.5;
% debxe(:,debctr)=xeste;
% debpe(:,debctr)=diag(Pe).^0.5;

%Start main loop
for in=1:nit
    %observations
    u=[0;0;0];

    y=M*u+Cimu*x+sRimu*randn(nsen,1);
    imu_inp=Mls*y;
    imu_obs=T*y;

    ye=[Mls_e1*y(1:5);y(6:7);Mls_e2*y(8:end)];
    imu_inpe=Mlse*ye;
    imu_obse=Te*ye;

    %Kalman
    K=P*Cimu_obs'*inv(Cimu_obs*P*Cimu_obs'+Rimu_obs);
    P=P-K*Cimu_obs*P;
    dz=K*(imu_obs-Cimu_obs*xest);
    xest=xest+dz;

    Ke=Pe*Cimu_obse'*inv(Cimu_obse*Pe*Cimu_obse'+Rimu_obse);
    Pe=Pe-Ke*Cimu_obse*Pe;
    dze=Ke*(imu_obse-Cimu_obse*xeste);
    xeste=xeste+dze;

    mx_a=Mls*Cimu*P*Cimu'*Mls';
    mx_b=Mlse*Cimue*Pe*Cimue'*Mlse';
    mx_c=Mls*Cimu*Pwls*Cimu'*Mls';

    debctr=debctr+1;
    debx(:,debctr)=xest;
    debp(:,debctr)=diag(P).^0.5;
    debxe(:,debctr)=xeste;
    debpe(:,debctr)=diag(Pe).^0.5;
    debuopt(:,debctr)=Mls*y-Mls*Cimu*xest;
    debuopte(:,debctr)=Mlse*ye-Mlse*Cimue*xeste;
    debuwls(:,debctr)=Mls*(y);
    debpopt(:,debctr)=diag(mx_a+Rimu_inp).^0.5;
    debpopte(:,debctr)=diag(mx_b+Rimu_inpe).^0.5;
    debpwls(:,debctr)=diag(mx_c+Rimu_inp).^0.5;
    debeq(1:2,debctr)=Mls_e1*y(1:5);
    debeq(3,debctr)=Mls_e3*y(6:7);
    debeq(4:5,debctr)=Mls_e2*y(8:end);

    %Compute the states of next cycle
    x=A*x+sQ*randn(nst,1);

    %prediction for the next states
    xest=A*xest;
    P=A*P*A'+Q;
    Pwls=A*Pwls*A'+Q;

    xeste=Ae*xeste;
    Pe=Ae*Pe*Ae'+Qe;

end

figure(10);
yy=1;
tt=[1:10:debctr]/60;
plot(tt,debuopte(yy,1:10:debctr));
hold on;
plot(tt,debuopt(yy,1:10:debctr));
plot(tt,debuwls(yy,1:10:debctr),'k');
plot(tt,debeq(0+yy,1:10:debctr),'r');
plot(tt,debeq(3+yy,1:10:debctr),'m');

plot(tt,debpopte(yy,1:10:debctr),'--');
plot(tt,debpopt(yy,1:10:debctr),'--');
plot(tt,debpwls(yy,1:10:debctr),'--k');
grid;
xlabel('Time (min)');
ylabel('X-Acceleration Error(m/sec^2)');
%legend('Opt. S. (u^{opt})', 'Red. Opt. S. (u^{red})','WLS S. (u^{wls})','Set1 Opt. (u^{A1})','Set2 Opt. (u^{A2})','Opt. SD.', 'Red. Opt. SD.','WLS SD');
legend('Opt. S. (u^{opt})', 'Red. Opt. S. (u^{red})','WLS S. (u^{wls})','Set1 Opt. (u^{A1})','Set2 Opt. (u^{A2})');

figure(11);
yy=2;
tt=[1:10:debctr]/60;
plot(tt,debuopte(yy,1:10:debctr));
hold on;
plot(tt,debuopt(yy,1:10:debctr));
plot(tt,debuwls(yy,1:10:debctr),'k');
plot(tt,debeq(0+yy,1:10:debctr),'r');
plot(tt,debeq(3+yy,1:10:debctr),'m');

plot(tt,debpopte(yy,1:10:debctr),'--');
plot(tt,debpopt(yy,1:10:debctr),'--');
plot(tt,debpwls(yy,1:10:debctr),'--k');
grid;
xlabel('Time (min)');
ylabel('Y-Acceleration Error(m/sec^2)');
%legend('Opt. S. (u^{opt})', 'Red. Opt. S. (u^{red})','WLS S. (u^{wls})','Set1 Opt. (u^{A1})','Set2 Opt. (u^{A2})','Opt. SD.', 'Red. Opt. SD.','WLS SD');
legend('Opt. S. (u^{opt})', 'Red. Opt. S. (u^{red})','WLS S. (u^{wls})','Set1 Opt. (u^{A1})','Set2 Opt. (u^{A2})');

return;

repo: ekfmonoslam


text:
% SOLVE_LINEAR_ELASTICITY: Solve a linear elasticity problem on a NURBS domain.
%                          For a planar domain it is the plane strain model.
%
% The function solves the linear elasticity problem
%
%      - div (sigma(u)) = f    in Omega = F((0,1)^n)
%      sigma(u) \cdot n = g    on Gamma_N
%                     u = h    on Gamma_D
%
% with   sigma(u) = mu*(grad(u) + grad(u)^t) + lambda*div(u)*I.
%
%   u:          displacement vector
%   sigma:      Cauchy stress tensor
%   lambda, mu: Lame' parameters
%   I:          identity tensor
%
% USAGE:
%
%  [geometry, msh, space, u] = solve_linear_elasticity (problem_data, method_data)
%
% INPUT:
%
%  problem_data: a structure with data of the problem. It contains the fields:
%    - geo_name:     name of the file containing the geometry
%    - nmnn_sides:   sides with Neumann boundary condition (may be empty)
%    - drchlt_sides: sides with Dirichlet boundary condition
%    - press_sides:  sides with pressure boundary condition (may be empty)
%    - symm_sides:   sides with symmetry boundary condition (may be empty)
%    - lambda_lame:  first Lame' parameter
%    - mu_lame:      second Lame' parameter
%    - f:            source term
%    - h:            function for Dirichlet boundary condition
%    - g:            function for Neumann condition (if nmnn_sides is not empty)
%
%  method_data : a structure with discretization data. Its fields are:
%    - degree:     degree of the spline functions.
%    - regularity: continuity of the spline functions.
%    - nsub:       number of subelements with respect to the geometry mesh
%                   (nsub=1 leaves the mesh unchanged)
%    - nquad:      number of points for Gaussian quadrature rule
%
% OUTPUT:
%
%  geometry: geometry structure (see geo_load)
%  msh:      mesh object that defines the quadrature rule (see msh_cartesian)
%  space:    space object that defines the discrete basis functions (see sp_vector)
%  u:        the computed degrees of freedom
%
% See also EX_LIN_ELAST_HORSESHOE for an example.
%
% Copyright (C) 2010 Carlo de Falco
% Copyright (C) 2011, 2015 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 [geometry, msh, sp, u] = ...
              solve_linear_elasticity (problem_data, method_data)

% Extract the fields from the data structures into local variables
data_names = fieldnames (problem_data);
for iopt  = 1:numel (data_names)
  eval ([data_names{iopt} '= problem_data.(data_names{iopt});']);
end
data_names = fieldnames (method_data);
for iopt  = 1:numel (data_names)
  eval ([data_names{iopt} '= method_data.(data_names{iopt});']);
end

% Construct geometry structure
geometry = geo_load (geo_name);
degelev  = max (degree - (geometry.nurbs.order-1), 0);
nurbs    = nrbdegelev (geometry.nurbs, degelev);
[rknots, zeta, nknots] = kntrefine (nurbs.knots, nsub-1, nurbs.order-1, regularity);

nurbs    = nrbkntins (nurbs, nknots);
geometry = geo_load (nurbs);

% Construct msh structure
rule     = msh_gauss_nodes (nquad);
[qn, qw] = msh_set_quad_nodes (geometry.nurbs.knots, rule);
msh      = msh_cartesian (geometry.nurbs.knots, qn, qw, geometry);

% Construct space structure
space_scalar = sp_nurbs (nurbs, msh);
scalar_spaces = repmat ({space_scalar}, 1, msh.rdim);
sp = sp_vector (scalar_spaces, msh);
clear space_scalar scalar_spaces

% Assemble the matrices
mat    = op_su_ev_tp (sp, sp, msh, lambda_lame, mu_lame);
rhs    = op_f_v_tp (sp, msh, f);

% Apply Neumann boundary conditions
for iside = nmnn_sides
% Restrict the function handle to the specified side, in any dimension, gside = @(x,y) g(x,y,iside)
  gside = @(varargin) g(varargin{:},iside);
  dofs = sp.boundary(iside).dofs;
  rhs(dofs) = rhs(dofs) + op_f_v_tp (sp.boundary(iside), msh.boundary(iside), gside);
end

% Apply pressure conditions
for iside = press_sides
  msh_side = msh_eval_boundary_side (msh, iside);
  sp_side  = sp_eval_boundary_side (sp, msh_side);

  x = cell (msh_side.rdim, 1);
  for idim = 1:msh_side.rdim
    x{idim} = reshape (msh_side.geo_map(idim,:,:), msh_side.nqn, msh_side.nel);
  end
  pval = reshape (p (x{:}, iside), msh_side.nqn, msh_side.nel);

  rhs(sp_side.dofs) = rhs(sp_side.dofs) - op_pn_v (sp_side, msh_side, pval);
end

% Apply symmetry conditions
symm_dofs = [];
for iside = symm_sides
  if (~strcmpi (sp.transform, 'grad-preserving'))
    error ('The symmetry condition is only implemented for spaces with grad-preserving transform')
  end
  msh_side = msh_eval_boundary_side (msh, iside);
  for idim = 1:msh.rdim
    normal_comp(idim,:) = reshape (msh_side.normal(idim,:,:), 1, msh_side.nqn*msh_side.nel);
  end

  parallel_to_axes = false;
  for ind = 1:msh.rdim
    ind2 = setdiff (1:msh.rdim, ind);
    if (all (all (abs (normal_comp(ind2,:)) < 1e-10)))
      symm_dofs = union (symm_dofs, sp.boundary(iside).dofs(sp.boundary(iside).comp_dofs{ind}));
      parallel_to_axes = true;
      break
    end
  end
  if (~parallel_to_axes)
    error ('solve_linear_elasticity: We have only implemented the symmetry condition for boundaries parallel to the axes')
  end

end

% Apply Dirichlet boundary conditions
u = zeros (sp.ndof, 1);
[u_drchlt, drchlt_dofs] = sp_drchlt_l2_proj (sp, msh, h, drchlt_sides);
u(drchlt_dofs) = u_drchlt;

int_dofs = setdiff (1:sp.ndof, [drchlt_dofs, symm_dofs]);
rhs(int_dofs) = rhs(int_dofs) - mat (int_dofs, drchlt_dofs) * u_drchlt;

% Solve the linear system
u(int_dofs) = mat(int_dofs, int_dofs) \ rhs(int_dofs);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


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

%%
clear; close all; clc;

%% Syntax
% |[A,B,C]=plane_fit(x,y,z);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function rollPitchYaw = rollPitchYawFromRotation(R)

    % ROLLPITCHYAWFROMROTATION converts a rotation matrix into Euler angles.
    %                          The Euler angles convention follows the one
    %                          of iDyntree and is such that the rotation
    %                          matrix is:  R = Rz(yaw)*Ry(pitch)*Rx(roll).
    %
    % FORMAT: rollPitchYaw = rollPitchYawFromRotation(R)
    %
    % INPUT:  - R = [3 * 3] rotation matrix
    %
    % OUTPUT: - rollPitchYaw = [3 * 1] vector of Euler angles [rad]
    %
    % Authors: Daniele Pucci, Marie Charbonneau, Gabriele Nava
    %
    %          all authors are with the Italian Istitute of Technology (IIT)
    %          email: name.surname@iit.it
    %
    % Genoa, Dec 2017
    %

    %% --- Initialization ---

    % For documentation, see also:
    %
    % http://wiki.icub.org/codyco/dox/html/idyntree/html/classiDynTree_1_1Rotation.html#a600352007d9250f7f227f21db85611f2
    %
    % http://www.geometrictools.com/Documentation/EulerAngles.pdf
    %
    rollPitchYaw = zeros(3,1);

    if (R(3,1) < +1)

        if (R(3,1) > -1)
            rollPitchYaw(2) = asin(-R(3,1));
            rollPitchYaw(3) = atan2(R(2,1),R(1,1));
            rollPitchYaw(1) = atan2(R(3,2), R(3,3));
        else
            % Not a unique solution : roll − yaw = atan2(−R23,R22)
            rollPitchYaw(2) = pi/2;
            rollPitchYaw(3) =-atan2(-R(2,3),R(2,2));
            rollPitchYaw(1) = 0;
        end
    else
        % Not a unique solution : roll − yaw = atan2(−R23,R22)
        rollPitchYaw(2) = -pi/2;
        rollPitchYaw(3) = atan2(-R(2,3),R(2,2));
        rollPitchYaw(1) = 0;
    end
end

author: robotology

repo: whole-body-controllers


text:
function transmat = ex_transmat(nt)
%EX_TRANSMAT  Example transition probability matrix definition.

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

transmat = cell(1, nt);
T = [ 0.158655253931457; 0.682689492137086; 0.158655253931457 ];
[transmat{:}] = deal(T * ones(1,3));
transmat{1} = T;


text:
 function y = reale(x, arg2)
%function y = reale(x)
%function y = reale(x, tol)
%function y = reale(x, 'warn')
%function y = reale(x, 'error')
%function y = reale(x, 'report')
%function y = reale(x, 'disp')
% return real part of complex data (with error checking).
% checks that imaginary part is negligible (or warning etc. if not)
% Copyright Jeff Fessler, The University of Michigan

com = 'error';
tol = 1e-13;
if nargin > 1
	if ischar(arg2)
		com = arg2;
	elseif isnumeric(arg2)
		tol = arg2;
	end
end

if streq(com, 'disp')
	;
elseif streq(com, 'warn')
	onlywarn = 1;
elseif streq(com, 'error')
	onlywarn = 0;
elseif streq(com, 'report')
	;
else
	error(sprintf('bad argument %s', com))
end

if max(abs(x(:))) == 0
	if any(imag(x(:)) ~= 0)
		error 'max real 0, but imaginary!'
	else
		y = real(x);
		return
end
end

frac = max(abs(imag(x(:)))) / max(abs(x(:)));
if streq(com, 'report')
	printf('imaginary part %g%%', frac * 100)
	return
end

if frac > tol
	t = sprintf('%s: imaginary fraction %g', mfilename, frac);
	if streq(com, 'disp')
		disp(t)
	elseif onlywarn
		disp(t)
	else
		error(t)
	end
end
y = real(x);

author: thomaskuestner


text:
%% write_mff_history.m
%  Matlab File
%  author Colin Davey
%  date 4/15/2014
%  Copyright 2014 EGI. All rights reserved.
%
%  This function adds an item to the list of history items. It is intended
%  to be run after running code that either a) creates a new mff file based
%  on an existing one, or b) modifies an existing one. The history item can
%  be seen in Net Station through the File Info->History feature.
%
%  filePath ? The path to the .mff file.
%
%  entryStruct has the following elements:
%
%  name - string that corresponds to the name of the tool specification in
%  NetStation. For example, if you create a segmentation specification for
%  a VTD experiment, and call it "VTD Seg", then the string "VTD Seg" would
%  go here.
%
%  method - The name of the tool (not the tool specification.) In the above
%  example, the string "Segmentation" would go here.
%
%  version - The version of the tool, eg "1.0".
%
%  beginTime - The time the tool started running. The following matlab
%  command creates a time in the proper format:
%    sprintf('%d-%02d-%02dT%02d:%02d:%02.5f',clock);
%
%  endTime - The time the tool finished running. See beginTime, above, for
%  the matlab command that creates a time in the proper format.
%
% sourceFileList - The list of files that were processed by the tool to
% generate the resulting file. For example, a tool that filters a single
% file to create a filtered version of the data would have just one item in
% this list. A tool that creates a grand average by averaging multiple
% single-subject files would have multiple items in this list.
%
% settingList - A cell array of zero or more strings that express the
% settings/parameters for the tool. These are optional, and are intended to
% be read by humans, not machines.
%
% resultList - A cell array of zero or more strings that express the
% results of the tool. These are optional, and are intended to be read by
% humans, not machines.
%%
function write_mff_history(filePath, entryStruct)
try
    mff_valid(filePath);
catch theException
    throw(theException);
end

%% 1) Read in history resource.
histObj = mff_getObject(com.egi.services.mff.api.MFFResourceType.kMFF_RT_History, 'history.xml', filePath);

%% 2) Get the list of current entires (it might be null, in which case
% create a new one and set it).
if isempty(histObj)
    delegate = javaObject('com.egi.services.mff.api.LocalMFFFactoryDelegate');
    factory = javaObject('com.egi.services.mff.api.MFFFactory', delegate);
    resourceVal = com.egi.services.mff.api.MFFResourceType.kMFF_RT_History;
    resourceType = javaObject('com.egi.services.mff.api.MFFResourceType', resourceVal);
    factory.createResourceAtURI([filePath filesep 'history.xml'], resourceType);
    histObj = factory.openResourceAtURI([filePath filesep 'history.xml'], resourceType);

    entryList = javaObject('java.util.ArrayList');
    histObj.setEntries(entryList);
else
    entryList = histObj.getEntries;
    if isempty(entryList) % Is this feasible? ie history exists, but entryList is empty?
        entryList = javaObject('java.util.ArrayList');
        histObj.setEntries(entryList);
    end
end

%% 3) Add an Entry object to the list.
entry = javaObject('com.egi.services.mff.api.Entry');
tool = javaObject('com.egi.services.mff.api.Tool');
tool.setName(entryStruct.name);
tool.setKind('Transformation');
tool.setMethod(entryStruct.method);
tool.setVersion(entryStruct.version);
tool.setBeginTime(entryStruct.beginTime);
tool.setEndTime(entryStruct.endTime);

sourceFileList = javaObject('java.util.ArrayList');
for p=1:size(entryStruct.sourceFileList,2)
    filePathObj = javaObject('com.egi.services.mff.api.FilePath');
    filePathObj.setFilePath(entryStruct.sourceFileList{p});
    sourceFileList.add(filePathObj);
    filePathObj = [];
end
java.lang.Runtime.getRuntime.freeMemory;
tool.setSourceFiles(sourceFileList);

settingList = javaObject('java.util.ArrayList');
for p=1:size(entryStruct.settingList,2)
    settingList.add(entryStruct.settingList{p});
end
tool.setSettings(settingList);

resultList = javaObject('java.util.ArrayList');
for p=1:size(entryStruct.resultList,2)
    resultList.add(entryStruct.resultList{p});
end
tool.setResults(resultList);

entry.setEntry(tool);
entry.setType('tool');

entryList.add(entry);

%% 4) Save the history resource in the standard way all resources are saved.
histObj.setEntries(entryList);
histObj.saveResource;


text:
function Model = buildRxnEquations(Model)
% Takes a COBRA model and creates reaction equations based
% on the metabolites list, `S` matrix and bound information. Returns the
% model with the additinal field `rxnEquations`. Is called by `cleanTmodel`, `compareCbModels`, `readCbTmodel`, `verifyModel`, `addSEEDInfo`.
%
% USAGE:
%
%    Model = buildRxnEquations(Model)
%
% INPUTS:
%    Model:     COBRA format model or `Tmodel`
%
% OUTPUTS:
%    Model:     Model with corrected `rxnEquations` field.
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

if isstruct(Model.lb) % If the function is fed Tmodel, use reaction bounds from first model.
    modelNames = fieldnames(Model.Models) ;
    firstModel = modelNames{1} ;
    lb = Model.lb.(firstModel) ;
    ub = Model.ub.(firstModel) ;
else
    lb = Model.lb ;
    ub = Model.ub ;
end

%% Make the equations.
for iRxn = 1:length(Model.rxns)
    rxnCode = [] ;

    educts   = find(Model.S(:,iRxn) < 0) ;
    products = find(Model.S(:,iRxn) > 0) ;
    if lb(iRxn) < 0 && ub(iRxn) > 0
        arrow = '<==>' ;
    elseif lb(iRxn) >= 0 && ub(iRxn) > 0
        arrow = '-->' ;
    elseif lb(iRxn) < 0 && ub(iRxn) <= 0
        arrow = '<--' ;
    elseif lb(iRxn) == 0 && ub(iRxn) == 0
        arrow = '|||' ;
    else
        disp(['Do not understand directionality of ' Model.rxns{iRxn}])
    end

    if ~isempty(educts)
        for ie = 1:length(educts)
            if Model.S(educts(ie), iRxn) ~= -1
                rxnCode = [rxnCode '(' num2str(-Model.S(educts(ie),iRxn)) ...
                          ') ' Model.mets{educts(ie)} ] ;
            else
                rxnCode = [rxnCode Model.mets{educts(ie)}] ;
            end
            rxnCode = [rxnCode ' + '] ;
        end
        % Note how just 2 positions are substracted from the end, leaving a
        % space.
        rxnCode = [rxnCode(1:end - 2) arrow] ;
    end
    if ~isempty(products)
        rxnCode = [rxnCode ' '] ;
        for ip = 1:length(products)
            if Model.S(products(ip), iRxn) ~= 1
                rxnCode = [rxnCode '(' num2str(Model.S(products(ip), iRxn))...
                          ') ' Model.mets{products(ip)}] ;
            else
                rxnCode = [rxnCode Model.mets{products(ip)}] ;
            end
            rxnCode = [rxnCode ' + '] ;
        end
        rxnCode = rxnCode(1:end - 3) ;
    end

    Model.rxnEquations{iRxn, 1} = rxnCode ;
    clear educts ; clear products; clear ie ; clear ip ; clear rxcode ;
end

author: opencobra

repo: cobratoolbox


text:
function events = in_events_nicolet(sFile, EventFile)
% IN_EVENTS_NICOLET: Open a text file with events exported from the Nicolet viewer.

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

% Open and read file
fid = fopen(EventFile,'r');
% Markers list
Markers = {};
isMarkerSection = 0;
% Read file line by line
while 1
    % Read one line
    newLine = fgetl(fid);
    if ~ischar(newLine)
        break;
    end
    % Lines to skip
    if isempty(newLine)
        continue;
    elseif ~isempty(strfind(newLine, 'Name')) && ~isempty(strfind(newLine, 'Time')) && ~isempty(strfind(newLine, 'Duration'))
        isMarkerSection = 1;
        continue
    elseif ~isempty(strfind(newLine, 'Exam Start')) || (nnz(newLine == 9) >= 2)    % Lines with tabs = events
        isMarkerSection = 1;
    elseif ~isMarkerSection
        continue;
    end
    % Split with tabs (ASCII #9): Name, Time, Duration
    splitLine = str_split(newLine, char(9), 0);
    if (length(splitLine) ~= 3) || any(cellfun(@isempty, splitLine))
        continue;
    end
    % Get the time and duration
    try
        vecStart = datevec(splitLine{2}, 'HH:MM:SS');
        tStart = vecStart(4)*3600 + vecStart(5)*60 + vecStart(6);
        vecDuration = datevec(splitLine{3}, 'MM:SS');
        tDuration = vecDuration(5)*60 + vecDuration(6);
    catch
        continue;
    end
    % Add markers entry: {name, type, start, length}
    Markers(end+1,:) = {splitLine{1}, tStart, tDuration};
end
% Close file
fclose(fid);

% List of events
if isempty(Markers)
    uniqueEvt = [];
else
    uniqueEvt = unique(Markers(:,1)');
end
% Initialize returned structure
events = repmat(db_template('event'), [1, length(uniqueEvt)]);
% Create events list
for iEvt = 1:length(uniqueEvt)
    % Find all the occurrences of event #iEvt
    iMrk = find(strcmpi(Markers(:,1)', uniqueEvt{iEvt}));
    % Add event structure
    events(iEvt).label   = uniqueEvt{iEvt};
    events(iEvt).epochs  = ones(1, length(iMrk));
    times = [Markers{iMrk,2}];
    if any([Markers{iMrk,3}] > 0)
        times(2,:) = [Markers{iMrk,2}] + [Markers{iMrk,3}];
    end
    events(iEvt).times      = round(times .* sFile.prop.sfreq) ./ sFile.prop.sfreq;
    events(iEvt).reactTimes = [];
    events(iEvt).select     = 1;
    events(iEvt).channels   = [];
    events(iEvt).notes      = [];
end

author: brainstorm-tools