text:
%% compare behavior detections for different wild-type strains

%% set up path

if ispc,
  addpath ../../FlyBowlAnalysis;
else
  addpath /groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis;
end
addpath ../perframe;
addpath ../perframe/compute_perframe_features;
addpath ../misc;
addpath ../filehandling;
outfigdir = '../figures/WildtypeOut';
if ~exist(outfigdir,'dir'),
  mkdir(outfigdir);
end

[~,computername] = system('hostname');
computername = strtrim(computername);

switch computername,

  case 'kabram-ws.janelia.priv',
    rootdatadir= '/groups/branson/home/bransonk/behavioranalysis/code/Jdetect/Jdetect/experiments/wildtype/results';

  otherwise,

    rootdatadir = '../experiments/wildtype/results';

end

if ispc,
  addpath C:\Code\FlyBowlAnalysis;
else
  addpath /groups/branson/bransonlab/projects/olympiad/FlyBowlAnalysis;
end

%% data locations

expnames = dir(fullfile(rootdatadir,'*201206*'));
expnames(~[expnames.isdir]) = [];
expnames = {expnames.name};
expdirs = cellfun(@(x) fullfile(rootdatadir,x),expnames,'UniformOutput',false);
nexps = numel(expdirs);

scoresfilestrs = dir(fullfile(expdirs{1},'*cores*.mat'));
scoresfilestrs = {scoresfilestrs.name};
scoresfns = cellfun(@(x) x(1:end-4),scoresfilestrs,'UniformOutput',false);
labelfns = regexprep(scoresfns,'scores','labels','preservecase','once');
nbehaviors = numel(labelfns);

%% parameters

linenamecodes = {
  'EXT_CSMH'             'CantonS MH'
  'EXT_CantonS_1220002'  'CantonS UH'
  'EXT_DL'               'Dickinson'
  'FCF_cantons_1500002'  'CantonS GR'
  'UAH_R-1220001'        'OregonR'
  'UAH_R-1220003'        'Berlin'
  };

behaviorcodes = {
  'labels_Walk'          'Walk'
  'labels_Stops'         'Stop'
  'labels_Jump'          'Jump'
  'labels_Righting'      'Righting'
  'labelsWingGrooming'   'Wing grooming'
  'labelsCrabwalk'       'Crabwalk'
  'labelsBackup'         'Backup'
  'labels_pivot_tail'    'Tail pivot turn'
  'labelsTouch'          'Touch'
  'labels_Chasev7'       'Chase'
  };

print_linenames_order = {'Berlin','CantonS MH','CantonS UH','CantonS GR','Dickinson','OregonR'};

%% compute statistics of behavior detections

meanfractime = nan(nexps,nbehaviors);
meanfractime_male = nan(nexps,nbehaviors);
meanfractime_female = nan(nexps,nbehaviors);

for expi = 1:nexps,

  % current experiment
  try
    trx = Trx('trxfilestr','registered_trx.mat','moviefilestr','movie.ufmf','perframedir','perframe');
    trx.AddExpDir(expdirs{expi},'openmovie',false);
  catch ME
    warning(getReport(ME));
    continue;
  end

  % number of frames analyzed
  nframes_total_perfly = trx.nframes;
  nframes_total = sum(nframes_total_perfly);
  nframes_male_perfly = cellfun(@(x) nnz(strcmp(x,'M')),trx.sex);
  nframes_female_perfly = cellfun(@(x) nnz(strcmp(x,'F')),trx.sex);
  nframes_male = sum(nframes_male_perfly);
  nframes_female = sum(nframes_female_perfly);

  % in per-fly std computation, ignore trajectories that are short
  ignoreidx = nframes_total_perfly <= 500;

  fprintf('%s:\n',expnames{expi});

  for behi = 1:nbehaviors,

    % mean, both sexes
    nframes_pos = nnz(cell2mat(trx.(labelfns{behi})));
    nframes_pos_perfly = cellfun(@nnz,trx.(labelfns{behi}));
    fractime_perfly = nframes_pos_perfly ./ nframes_total;
    meanfractime(expi,behi) = nframes_pos / nframes_total;

    % mean, per sex
    nframes_male_pos_perfly = nan(1,trx.nflies);
    nframes_female_pos_perfly = nan(1,trx.nflies);
    for fly = 1:trx.nflies,
      idxmale = strcmp(trx(fly).sex,'M');
      idxfemale = strcmp(trx(fly).sex,'F');
      nframes_male_pos_perfly(fly) = nnz(trx(fly).(labelfns{behi})(idxmale));
      nframes_female_pos_perfly(fly) = nnz(trx(fly).(labelfns{behi})(idxfemale));
    end
    nframes_male_pos = sum(nframes_male_pos_perfly);
    nframes_female_pos = sum(nframes_female_pos_perfly);
    meanfractime_male(expi,behi) = nframes_male_pos / nframes_male;
    meanfractime_female(expi,behi) = nframes_female_pos / nframes_female;

    % std TODO

    fprintf('%s: b = %f, m = %f, f = %f\n',labelfns{behi},...
      meanfractime(expi,behi),...
      meanfractime_male(expi,behi),meanfractime_female(expi,behi));

  end
end

%% z-score the data so that we can plot on the same axes

zmu = nanmean(meanfractime,1);
zsigma = nanstd(meanfractime,1,1);
zmeanfractime = bsxfun(@rdivide,bsxfun(@minus,meanfractime,zmu),zsigma);
zmu_male = nanmean(meanfractime_male,1);
zsigma_male = nanstd(meanfractime_male,1,1);
zmeanfractime_male = bsxfun(@rdivide,bsxfun(@minus,meanfractime_male,zmu_male),zsigma_male);
zmu_female = nanmean(meanfractime_female,1);
zsigma_female = nanstd(meanfractime_female,1,1);
zmeanfractime_female = bsxfun(@rdivide,bsxfun(@minus,meanfractime_female,zmu_female),zsigma_female);

%% compute per-strain data

for i = 1:nexps,
  metadatacurr = parseExpDir(expdirs{i});
  if i == 1,
    metadata = repmat(metadatacurr,[1,nexps]);
  else
    metadata(i) = metadatacurr;
  end
end

line_names = {metadata.line};
[unique_linenames,~,lineidx] = unique(line_names);

nlines = numel(unique_linenames);

meanfractime_line = nan(nlines,nbehaviors);
meanfractime_male_line = nan(nlines,nbehaviors);
meanfractime_female_line = nan(nlines,nbehaviors);
stdfractime_line = nan(nlines,nbehaviors);
stdfractime_male_line = nan(nlines,nbehaviors);
stdfractime_female_line = nan(nlines,nbehaviors);

for linei = 1:nlines,
  meanfractime_line(linei,:) = nanmean(meanfractime(lineidx==linei,:),1);
  stdfractime_line(linei,:) = nanstd(meanfractime(lineidx==linei,:),1,1);
  meanfractime_male_line(linei,:) = nanmean(meanfractime_male(lineidx==linei,:),1);
  stdfractime_male_line(linei,:) = nanstd(meanfractime_male(lineidx==linei,:),1,1);
  meanfractime_female_line(linei,:) = nanmean(meanfractime_female(lineidx==linei,:),1);
  stdfractime_female_line(linei,:) = nanstd(meanfractime_female(lineidx==linei,:),1,1);
end

zmeanfractime_line = bsxfun(@rdivide,bsxfun(@minus,meanfractime_line,zmu),zsigma);
zmeanfractime_male_line = bsxfun(@rdivide,bsxfun(@minus,meanfractime_male_line,zmu_male),zsigma_male);
zmeanfractime_female_line = bsxfun(@rdivide,bsxfun(@minus,meanfractime_female_line,zmu_female),zsigma_female);
zstdfractime_line = bsxfun(@rdivide,stdfractime_line,zsigma);
zstdfractime_male_line = bsxfun(@rdivide,stdfractime_male_line,zsigma_male);
zstdfractime_female_line = bsxfun(@rdivide,stdfractime_female_line,zsigma_female);

%% plot

hfig = 1;
figure(hfig);
clf;
set(hfig,'Position',[80 94 1716 854]);
colors = jet(nlines)*.7;

behaviornames = regexprep(labelfns,'(labels_)|(Labels_)|(labels)|(Labels)','');

combinedzmeanfractime = reshape([zmeanfractime;zmeanfractime_male;zmeanfractime_female],[nexps,nbehaviors*3]);
combinedzmeanfractime_line = reshape([zmeanfractime_line;zmeanfractime_male_line;zmeanfractime_female_line],[nlines,nbehaviors*3]);
combinedzstdfractime_line = reshape([zstdfractime_line;zstdfractime_male_line;zstdfractime_female_line],[nlines,nbehaviors*3]);
combinedbehaviornames = [behaviornames
  cellfun(@(x) [x,'_male'],behaviornames,'UniformOutput',false)
  cellfun(@(x) [x,'_female'],behaviornames,'UniformOutput',false)];
combinedbehaviornames = reshape(combinedbehaviornames,[1,nbehaviors*3]);

hdata = nan(1,nlines);
hold on;
for i = 1:nlines,
  plot((1:3*nbehaviors)-.025+.05*(i-1),combinedzmeanfractime(lineidx==i,:),'.','Color',colors(i,:)*.5+.5);
end
for i = 1:nlines,
  plot(repmat(1:3*nbehaviors,[2,1])-.025+.05*(i-1),...
    bsxfun(@plus,combinedzmeanfractime_line(i,:),[combinedzstdfractime_line(i,:);-combinedzstdfractime_line(i,:)]),...
    '-','Color',colors(i,:));
end
for i = 1:nlines,
  hdata(i) = plot((1:3*nbehaviors)-.025+.05*(i-1),combinedzmeanfractime_line(i,:),'o',...
    'Color',colors(i,:),'MarkerFaceColor',colors(i,:));
end
set(gca,'XTick',1:nbehaviors*3,'XTickLabel',combinedbehaviornames);
rotateticklabel(gca,90)

legend(hdata,unique_linenames,'interpreter','none');
set(gca,'XLim',[0,3*nbehaviors+1]);

ylabel('Z-scored fraction of time');

SaveFigLotsOfWays(hfig,fullfile(outfigdir,'WildtypeBehaviorComparison'));

%% plot by weekend

idxplot = ~cellfun(@isempty,regexp({metadata.date},'^2012060','once'));
meanfractime2 = meanfractime(idxplot,:);
meanfractime_male2 = meanfractime_male(idxplot,:);
meanfractime_female2 = meanfractime_female(idxplot,:);
metadata2 = metadata(idxplot);
nexps2 = nnz(idxplot);

%% z-score the data so that we can plot on the same axes

zmu = nanmean(meanfractime2,1);
zsigma = nanstd(meanfractime2,1,1);
zmeanfractime = bsxfun(@rdivide,bsxfun(@minus,meanfractime2,zmu),zsigma);
zmu_male = nanmean(meanfractime_male2,1);
zsigma_male = nanstd(meanfractime_male2,1,1);
zmeanfractime_male = bsxfun(@rdivide,bsxfun(@minus,meanfractime_male2,zmu_male),zsigma_male);
zmu_female = nanmean(meanfractime_female2,1);
zsigma_female = nanstd(meanfractime_female2,1,1);
zmeanfractime_female = bsxfun(@rdivide,bsxfun(@minus,meanfractime_female2,zmu_female),zsigma_female);

%% compute per-strain data

line_names = {metadata2.line};
[unique_linenames,~,lineidx] = unique(line_names);

nlines = numel(unique_linenames);

meanfractime_line = nan(nlines,nbehaviors);
meanfractime_male_line = nan(nlines,nbehaviors);
meanfractime_female_line = nan(nlines,nbehaviors);
stdfractime_line = nan(nlines,nbehaviors);
stdfractime_male_line = nan(nlines,nbehaviors);
stdfractime_female_line = nan(nlines,nbehaviors);

for linei = 1:nlines,
  meanfractime_line(linei,:) = nanmean(meanfractime2(lineidx==linei,:),1);
  stdfractime_line(linei,:) = nanstd(meanfractime2(lineidx==linei,:),1,1);
  meanfractime_male_line(linei,:) = nanmean(meanfractime_male2(lineidx==linei,:),1);
  stdfractime_male_line(linei,:) = nanstd(meanfractime_male2(lineidx==linei,:),1,1);
  meanfractime_female_line(linei,:) = nanmean(meanfractime_female2(lineidx==linei,:),1);
  stdfractime_female_line(linei,:) = nanstd(meanfractime_female2(lineidx==linei,:),1,1);
end

zmeanfractime_line = bsxfun(@rdivide,bsxfun(@minus,meanfractime_line,zmu),zsigma);
zmeanfractime_male_line = bsxfun(@rdivide,bsxfun(@minus,meanfractime_male_line,zmu_male),zsigma_male);
zmeanfractime_female_line = bsxfun(@rdivide,bsxfun(@minus,meanfractime_female_line,zmu_female),zsigma_female);
zstdfractime_line = bsxfun(@rdivide,stdfractime_line,zsigma);
zstdfractime_male_line = bsxfun(@rdivide,stdfractime_male_line,zsigma_male);
zstdfractime_female_line = bsxfun(@rdivide,stdfractime_female_line,zsigma_female);

%% plot

hfig = 2;
figure(hfig);
clf;
set(hfig,'Units','pixels','Position',[80 94 1716 854]);
%colors = jet(nlines)*.7;
colors = lines(nlines);

behaviornames = regexprep(labelfns,'(labels_)|(Labels_)|(labels)|(Labels)','');

combinedzmeanfractime = reshape([zmeanfractime;zmeanfractime_male;zmeanfractime_female],[nexps2,nbehaviors*3]);
combinedzmeanfractime_line = reshape([zmeanfractime_line;zmeanfractime_male_line;zmeanfractime_female_line],[nlines,nbehaviors*3]);
combinedzstdfractime_line = reshape([zstdfractime_line;zstdfractime_male_line;zstdfractime_female_line],[nlines,nbehaviors*3]);
combinedbehaviornames = [behaviornames
  cellfun(@(x) [x,'_male'],behaviornames,'UniformOutput',false)
  cellfun(@(x) [x,'_female'],behaviornames,'UniformOutput',false)];
combinedbehaviornames = reshape(combinedbehaviornames,[1,nbehaviors*3]);

hdata = nan(1,nlines);
hold on;
for i = 1:nlines,
  plot((1:3*nbehaviors)-.025+.05*(i-1),combinedzmeanfractime(lineidx==i,:),'.','Color',colors(i,:)*.5+.5);
end
for i = 1:nlines,
  plot(repmat(1:3*nbehaviors,[2,1])-.025+.05*(i-1),...
    bsxfun(@plus,combinedzmeanfractime_line(i,:),[combinedzstdfractime_line(i,:);-combinedzstdfractime_line(i,:)]),...
    '-','Color',colors(i,:));
end
for i = 1:nlines,
  hdata(i) = plot((1:3*nbehaviors)-.025+.05*(i-1),combinedzmeanfractime_line(i,:),'o',...
    'Color',colors(i,:),'MarkerFaceColor',colors(i,:));
end
set(gca,'XTick',1:nbehaviors*3,'XTickLabel',combinedbehaviornames);
rotateticklabel(gca,90)

legend(hdata,unique_linenames,'interpreter','none');
set(gca,'XLim',[0,3*nbehaviors+1]);

ylabel('Z-scored fraction of time');

SaveFigLotsOfWays(hfig,fullfile(outfigdir,'WildtypeBehaviorComparison_Weekend2'));

%% plot each behavior separately

[~,idx] = ismember(unique_linenames,linenamecodes(:,1));
print_linenames = linenamecodes(idx,2)';

[~,idx] = ismember(labelfns,behaviorcodes(:,1));
behaviors = behaviorcodes(idx,2)';

hfig = 3;
figure(hfig);
clf;
set(hfig,'Units','pixels','Position',[100,100,1400 1000]);
hax = createsubplots(5,2,[[.05,.05];[.05,.01]]);

hbar = nan(3,nbehaviors);

datacolors = [0,0,0;0,0,.5;.5,0,0];
datanames = {'Both','Male','Female'};
grp1 = repmat(print_linenames(lineidx)',[1,3]);
grp2 = repmat(datanames,[nexps2,1]);
grporder = cell(3,nlines);
for i = 1:nlines,
  for j = 1:3,
    grporder{j,i} = sprintf('%s,%s',print_linenames_order{i},datanames{j});
  end
end
grporder = grporder(:);

for behi = 1:nbehaviors,

  axes(hax(behi));

%   v = [meanfractime_line(:,behi),meanfractime_male_line(:,behi),meanfractime_female_line(:,behi)];
%   s = [stdfractime_line(:,behi),stdfractime_male_line(:,behi),stdfractime_female_line(:,behi)];
  p = [meanfractime2(:,behi),meanfractime_male2(:,behi),meanfractime_female2(:,behi)];
%   hbar(:,behi) = bar(1:nlines,v,'grouped');
%   for i = 1:3,
%     set(hbar(i,behi),'FaceColor',datacolors(i,:));
%   end
%   hold on;
%   for i = 1:3,
%     x = mean(get(get(hbar(i,behi),'Children'),'XData'),1);
%     herr(i,behi) = errorbar(x,v(:,i),s(:,i),'linestyle','none','color',[.5,.5,.5]);
%     hexp(i,behi) = plot(x(lineidx),p(:,i),'.','color',[0,.5,.5]);
%   end

  hb = boxplot(p(:),{grp1(:),grp2(:)},...
  'ColorGroup',grp2(:),...
  'GroupOrder',grporder,...
  'PlotStyle','traditional',...
  'FactorDirection','list',...
  'labelverbosity','minor',...
  'colors',datacolors,...
  'symbol','+');

  hxlabel = findobj(get(hb(1),'Parent'),'Type','text');
  delete(hxlabel);

  if behi == 1,
    for i = 1:3,
      hdcurr = findobj(get(hb(1),'Parent'),'Type','line','Color',datacolors(i,:),'Marker','none','LineStyle','-');
      hd(i) = hdcurr(1);
    end
    legend(hd,datanames);
  end
  ylabel(sprintf('Frac. time %s',lower(behaviors{behi})));

  if mod(behi,5) == 0,
    set(gca,'XTick',2:3:3*nlines,'XTickLabel',print_linenames_order);
  else
    set(gca,'XTick',2:3:3*nlines,'XTickLabel',{});
  end

  %set(gca,'YLim',[0,max(v(:)+s(:))*1.025]);
  %set(gca,'XLim',[.5,nlines+.5]);

end

SaveFigLotsOfWays(hfig,fullfile(outfigdir,'WildtypeBehaviorComparison_Weekend2_PerBehaviorBoxPlots'));

%% plot each behavior separately, only plot both statistics

[~,idx] = ismember(unique_linenames,linenamecodes(:,1));
print_linenames = linenamecodes(idx,2)';

% [~,idx] = ismember(labelfns,behaviorcodes(:,1));
% behaviors = behaviorcodes(idx,2)';

ndx2meanfrac = [];
behaviors = {};
for ndx = 1:numel(behaviorcodes(:,1))
  ndx2labelfns = find(cellfun(@(x) strcmp(x,behaviorcodes{ndx,1}),labelfns));
  behaviors{ndx} = behaviorcodes{ndx,2};
  ndx2meanfrac(ndx) = ndx2labelfns;
end
nbehaviors = numel(behaviors);

hfig = 4;
figure(hfig);
clf;
set(hfig,'Units','pixels','Position',[100,100,1400 1000]);
hax = createsubplots(5,2,[[.05,.05];[.05,.01]]);

hbar = nan(3,nbehaviors);

%datacolors = [0,0,0;0,0,.5;.5,0,0];
%datanames = {'Both','Male','Female'};
grp1 = print_linenames(lineidx)';
%grp1 = repmat(print_linenames(lineidx)',[1,3]);
%grp2 = repmat(datanames,[nexps2,1]);
grporder = print_linenames_order;
% grporder = cell(3,nlines);
% for i = 1:nlines,
%   for j = 1:3,
%     grporder{j,i} = sprintf('%s,%s',print_linenames_order{i},datanames{j});
%   end
% end
% grporder = grporder(:);

for behi = 1:nbehaviors,

  axes(hax(behi));

%   v = [meanfractime_line(:,behi),meanfractime_male_line(:,behi),meanfractime_female_line(:,behi)];
%   s = [stdfractime_line(:,behi),stdfractime_male_line(:,behi),stdfractime_female_line(:,behi)];
  %p = [meanfractime2(:,behi),meanfractime_male2(:,behi),meanfractime_female2(:,behi)];
  p = [meanfractime2(:,ndx2meanfrac(behi))];
%   hbar(:,behi) = bar(1:nlines,v,'grouped');
%   for i = 1:3,
%     set(hbar(i,behi),'FaceColor',datacolors(i,:));
%   end
%   hold on;
%   for i = 1:3,
%     x = mean(get(get(hbar(i,behi),'Children'),'XData'),1);
%     herr(i,behi) = errorbar(x,v(:,i),s(:,i),'linestyle','none','color',[.5,.5,.5]);
%     hexp(i,behi) = plot(x(lineidx),p(:,i),'.','color',[0,.5,.5]);
%   end

  hb = boxplot(p(:),grp1(:),...
  'GroupOrder',grporder,...
  'colors',[0,0,0],...
  'symbol','+');

  hxlabel = findobj(get(hb(1),'Parent'),'Type','text');
  delete(hxlabel);

%   if behi == 1,
%     for i = 1:3,
%       hdcurr = findobj(get(hb(1),'Parent'),'Type','line','Color',datacolors(i,:),'Marker','none','LineStyle','-');
%       hd(i) = hdcurr(1);
%     end
%     legend(hd,datanames);
%   end
  ylabel(sprintf('Frac. time %s',lower(behaviors{behi})));

  if mod(behi,5) == 0,
    set(gca,'XTick',1:nlines,'XTickLabel',print_linenames_order);
  else
    set(gca,'XTick',1:nlines,'XTickLabel',{});
  end
  set(gca,'box','off');
  set(get(gca,'XLabel'),'FontSize',16);
  set(get(gca,'YLabel'),'FontSize',16);
  set(gca,'FontSize',16);

  %set(gca,'YLim',[0,max(v(:)+s(:))*1.025]);
  %set(gca,'XLim',[.5,nlines+.5]);

end

SaveFigLotsOfWays(hfig,fullfile(outfigdir,'WildtypeBehaviorComparison_Weekend2_PerBehaviorBoxPlots_BothOnly'));

%% Save the Data..

save(fullfile(outfigdir,['WildtypeBehaviorComparison_ResultsData']),'behaviornames','labelfns',...
  'meanfractime2','metadata2','scoresfns');

%% use logistic regression to predict the strain given the behavior classification results for each pair

% reorder line names
print_linenames_perexp = print_linenames(lineidx);
[~,orderedlineidx] = ismember(print_linenames_perexp,print_linenames_order);

% X is nexps x nbehaviors
X = meanfractime2;
y = orderedlineidx';
isbaddata = any(isnan(X),2);
X(isbaddata,:) = [];
y(isbaddata) = [];
errorrate = zeros(nlines,nlines);
errorrate(:) = .5;
for line1 = 1:nlines,
  idx1 = line1 == y;
  for line2 = line1+1:nlines,
    idx2 = line2 == y;
    idxcurr = idx1 | idx2;
    ncurr = nnz(idxcurr);
    Xcurr = X(idxcurr,:);
    % 1 means line1, 0 means line2
    ycurr = idx1(idxcurr);
    % loop over all examples for cross validation

    yfitcurr = nan(ncurr,1);
    for i = 1:ncurr,

      % logistic regression
      idxtrain = true(ncurr,1); idxtrain(i) = false;
      coeffscurr = glmfit(Xcurr(idxtrain,:),ycurr(idxtrain),'binomial','link','logit');
      yfitcurr(i) = glmval(coeffscurr,Xcurr(i,:),'logit');

    end

    errorcurr = nnz( (yfitcurr>.5) ~= ycurr );
    errorrate(line1,line2) = errorcurr / ncurr;
    errorrate(line2,line1) = errorrate(line1,line2);
  end
end

hfig = 5;
figure(hfig);
clf;

imagesc(errorrate,[0,.5]);
colormap gray;
set(gca,'XTick',1:nlines,'XTickLabel',print_linenames_order,...
  'YTick',1:nlines,'YTickLabel',print_linenames_order,'Box','off');
axis image;
colorbar;
title('Pairwise error rate');

errorrate0 = errorrate;

%% use regularized logistic regression to predict the strain given the behavior classification results for each pair

% reorder line names
print_linenames_perexp = print_linenames(lineidx);
[~,orderedlineidx] = ismember(print_linenames_perexp,print_linenames_order);
nlambda = 10;

% X is nexps x nbehaviors
X = meanfractime2;
y = orderedlineidx';
isbaddata = any(isnan(X),2);
X(isbaddata,:) = [];
y(isbaddata) = [];
errorrate = zeros(nlines,nlines);
errorrate(:) = .5;
for line1 = 1:nlines,
  idx1 = line1 == y;
  for line2 = line1+1:nlines,
    idx2 = line2 == y;
    idxcurr = idx1 | idx2;
    ncurr = nnz(idxcurr);
    Xcurr = X(idxcurr,:);
    % 1 means line1, 0 means line2
    ycurr = idx1(idxcurr);
    % loop over all examples for cross validation

    yfitcurr = nan(ncurr,1);
    for i = 1:ncurr,

      % logistic regression
      idxtrain = true(ncurr,1); idxtrain(i) = false;
      [coeffscurr,fitinfo] = lassoglm(Xcurr(idxtrain,:),ycurr(idxtrain),'binomial','link','logit',...
        'NumLambda',nlambda,'CV',ncurr-1);
      j = fitinfo.Index1SE;
      coeffscurr0 = coeffscurr(:,j);
      cnst = fitinfo.Intercept(j);
      coeffscurr1 = [cnst;coeffscurr0];
      yfitcurr(i) = glmval(coeffscurr1,Xcurr(i,:),'logit');
    end

    errorcurr = nnz( (yfitcurr>.5) ~= ycurr );
    errorrate(line1,line2) = errorcurr / ncurr;
    errorrate(line2,line1) = errorrate(line1,line2);
  end
end

hfig = 6;
figure(hfig);
clf;

imagesc(errorrate,[0,.5]);
colormap gray;
set(gca,'XTick',1:nlines,'XTickLabel',print_linenames_order,...
  'YTick',1:nlines,'YTickLabel',print_linenames_order,'Box','off');
xlabel('Wildtype strains');ylabel('Wildtype strains');
axis image;
colorbar;
title('Pairwise error rate');

error_rate = errorrate;
axis image;
box off;
numel(error_rate)
size(error_rate)
[Row,Col]= size(error_rate);
for i = 1:Row
  for j = 1:Col
    if i==j; continue; end
    if errorrate(i,j) ==0; continue; end
    ht(i,j) = text(j,i,sprintf('%.2f',error_rate(i,j)),...
      'FontUnits','pixels','FontWeight','bold','Color',[.7,0,0],...
      'HorizontalAlignment','center','VerticalAlignment','middle',...
      'FontSize',24);
  end
end

SaveFigLotsOfWays(hfig,fullfile(outfigdir,'WildTypeStrainPrediction_PairwiseErrorRate'));

save('WildTypeStrainPredictionData.mat','errorrate','print_linenames_order');

%     'Berlin'    'CantonS MH'    'CantonS UH'    'CantonS GR'    'Dickinson'    'OregonR'
% errorrate =
%     0.5000         0    0.0714    0.0625         0         0
%          0    0.5000    0.0714         0         0         0
%     0.0714    0.0714    0.5000         0         0         0
%     0.0625         0         0    0.5000    0.0625         0
%          0         0         0    0.0625    0.5000    0.0588
%          0         0         0         0    0.0588    0.5000

author: kristinbranson

repo: JAABA


text:
function D = DateStr2Num(S, F)
% DATESTR2NUM - Fast conversion of DATESTR to DATENUM
% The builtin DATENUM command is very powerful, but if you know the input
% format exactly, a specific MEX can be much faster: e.g. for single strings
% DateStr2Num is 50 to 100 times faster than DATENUM, for a {1 x 10000} cell
% string the speed up factor is 300 to 700 (Matlab 2009a/2011b/64, MSVC 2008).
%
% D = DateStr2Num(S, F)
% INPUT:
%   S: String or cell string in DATESTR(F) format.
%      In opposite to DATENUM the validity of the input string is not checked
%      (e.g. 1 <= month <= 12).
%   F: Integer number defining the input format. Accepted:
%          0:  'dd-mmm-yyyy HH:MM:SS'      01-Mar-2000 15:45:17
%          1:  'dd-mmm-yyyy'               01-Mar-2000
%         29:  'yyyy-mm-dd'                2000-03-01
%         30:  'yyyymmddTHHMMSS'           20000301T154517
%         31:  'yyyy-mm-dd HH:MM:SS'       2000-03-01 15:45:17
%      Including the milliseconds (not a valid DATEFORM number in DATESTR!):
%        1000: 'dd-mmm-yyyy HH:MM:SS.FFF'  01-Mar-2000 15:45:17.123
%        1030: 'yyyymmddTHHMMSS.FFF'       20000301T154517.123
%
%      Optional, default: 0.
%
% OUTPUT:
%   D: Serial date number. If S is a cell, D has is same size.
%
% EXAMPLES:
%   C = {'2010-06-29 21:59:13', '2010-06-29 21:59:13'};
%   D = DateStr2Num(C, 31)
%   >> [734318.916122685, 734318.916122685]
%   Equivalent Matlab command (but a column vector is replied ever):
%   D = datenum(C, 'yyyy-mm-dd HH:MM:SS')
%
% NOTES: The parsing of the strings works for clean ASCII characters only:
%   '0' must have the key code 48!
%   Month names must be English with the 2nd and 3rd charatcer in lower case.
%
% COMPILATION:
% The C-function must be compiled before using. The M-file produces the same
% results, but is slower.
%   Windows: mex -O DateStr2Num.c
%   Linux:   mex -O CFLAGS="\$CFLAGS -std=c99" DateStr2Num.c
%   Download precompiled Mex: http://www.n-simon.de/mex
%
% Tested: Matlab 6.5, 7.7, 7.8, 7.13, WinXP/32, Win7/64
% Compiler: LCC 2.4/3.8, BCC 5.5, Open Watcom 1.8, MSVC 2008
%         Compatibility to MacOS, Linux, 64 bit is assumed, but not tested.
% Author: Jan Simon, Heidelberg, (C) 2010-2013 matlab.THISYEAR(a)nMINUSsimon.de
%
% See also DATESTR, DATENUM, DATEVEC.
% FEX: DateConvert 25594 (Jan Simon)

% $JRev: R-h V:008 Sum:gZffjdt4JfVU Date:28-Jul-2013 15:03:20 $
% $License: BSD (use/copy/change/redistribute on own risk, mention the author) $
% $UnitTest: uTest_DateStr2Num $
% $File: Tools\GLString\DateStr2Num.m $
% History:
% 001: 29-Jun-2010 22:07, First version.
% 002: 30-Jun-2010 16:08, Accept formats 0, 1, 29, 30, 31.
% 008: 28-Jul-2013 12:48, 300 -> 1030, format 1000 added.

% ==============================================================================
% The MEX is not found?! Perhaps the parent folder of this function is the
% current folder...
% peristent WarnOnce     % Uncomment this to warn once only
% if isempty(WarnOnce)
warning('JSimon:DateStr2Num:NoMex', ...
   ['The slow M-version of DateStr2Num is used\n', ...
   'Read the help text -> COMPILATION']);
%   WarnOnce = true;
% end

% Compared to DATENUM, this M-version is 20 times faster than DATENUM(S, 31) for
% a string and 1.1 to 1.4 times slower for a {1 x 10000} cell string.
if ischar(S)
   S = {S};
elseif ~iscellstr(S)
   error('JSimon:DateStr2Num:BadInputType', ...
      'Input must be a string or cell string.');
end

if nargin == 1
   F = 0;
end

nS = numel(S);
D  = zeros(size(S));
switch F
   case 0  % 'dd-mmm-yyyy HH:MM:SS'
      for iS = 1:nS
         aS    = S{iS};
         x     = aS - '0';
         month = (strfind('janfebmaraprmayjunjulaugsepoctnovdec', ...
            lower(aS(4:6))) + 2) / 3;
         year  = x(8) * 1000 + x(9) * 100 + x(10) * 10 + x(11);

         D(iS) = datenummx(year, month, x(1) * 10 + x(2), ...
            x(13) * 10 + x(14), x(16) * 10 + x(17), ...
            x(19) * 10 + x(20));
      end

   case 1  % 'dd-mmm-yyyy'
      for iS = 1:nS
         aS    = S{iS};
         x     = aS - '0';
         month = (strfind('janfebmaraprmayjunjulaugsepoctnovdec', ...
            lower(aS(4:6))) + 2) / 3;
         year  = x(8) * 1000 + x(9) * 100 + x(10) * 10 + x(11);

         D(iS) = datenummx(year, month, x(1) * 10 + x(2));
      end

   case 29  % 'yyyy-mm-dd'
      for iS = 1:nS
         x     = S{iS} - '0';
         D(iS) = datenummx( ...
            x(1)*1000 + x(2)*100 + x(3)*10 + x(4), ... % Year
            x(6)*10 + x(7), ...  % Month
            x(9)*10 + x(10));     % Day
      end

   case 30  % 'yyyymmddTHHMMSS'
      for iS = 1:nS
         x     = S{iS} - '0';
         D(iS) = datenummx( ...
            x(1)*1000 + x(2)*100 + x(3)*10 + x(4), ... % Year
            x(5)*10 + x(6), ...  % Month
            x(7)*10 + x(8), ...  % Day
            x(10)*10 + x(11), x(12)*10 + x(13), x(14)*10 + x(15));
      end

   case 31  % 'yyyy-mm-dd HH:MM:SS'
      for iS = 1:nS
         x     = S{iS} - '0';
         D(iS) = datenummx( ...
            x(1)*1000 + x(2)*100 + x(3)*10 + x(4), ... % Year
            x(6)*10 + x(7), ...  % Month
            x(9)*10 + x(10), ... % Day
            x(12)*10 + x(13), x(15)*10 + x(16), x(18)*10 + x(19));
      end

   case 1000  % 'dd-mmm-yyyy HH:MM:SS.FFF'
      for iS = 1:nS
         aS    = S{iS};
         x     = aS - '0';
         month = (strfind('janfebmaraprmayjunjulaugsepoctnovdec', ...
            lower(aS(4:6))) + 2) / 3;
         year  = x(8) * 1000 + x(9) * 100 + x(10) * 10 + x(11);

         D(iS) = datenummx(year, month, x(1) * 10 + x(2), ...
            x(13) * 10 + x(14), x(16) * 10 + x(17), ...
            x(19) * 10 + x(20) + ...
            x(22) / 10 + x(23) / 100 + x(24) / 1000);
      end

   case 1030  % yyyymmddTHHMMSS.FFF
      for iS = 1:nS
         x     = S{iS} - '0';
         D(iS) = datenummx( ...
            x(1)*1000 + x(2)*100 + x(3)*10 + x(4), ... % Year
            x(5)*10 + x(6), ...  % Month
            x(7)*10 + x(8), ...  % Day
            x(10)*10 + x(11), x(12)*10 + x(13), x(14)*10 + x(15) + ...
            x(17) / 10 + x(18) / 100 + x(19) / 1000);
      end

   otherwise
      error('JSimon:DateStr2Num:BadFormat', 'Bad format');
end

% return;

author: Sable

repo: mcbench-benchmarks


text:
function [gX, kern] = pathKernDiagGradSym(kern, X)

% PATHKERNDIAGGRADSYM Compute the gradient of the PATH kernel's diagonal wrt parameters.
% FORMAT
% DESC computes the gradient of functions of the diagonal of the
% path kernel matrix with respect to the parameters of the
% symmetric kernel.
% ARG kern : the kernel structure for which the gradients are
% computed.
% ARG x : the input data for which the gradient is being computed.
% RETURN gX : gradients of the relevant function with respect to each
% of the parameters.
% RETURN kern : the updated kernel structure
%
% SEEALSO : pathKernParamInit, pathDiagGradient, pathKernExtractParam, pathKernGradient
%
% COPYRIGHT : Andrea Baisero, Carl Henrik Ek, 2013

% SHEFFIELDML

maxl=max(cellfun(@(x)size(x,1),X));
kern=pathKernUpdateWMat(kern,maxl);

num=length(X);

gX=cell(1,num);
for i=1:num
    gX{i}=pathKernDiagGradSymSeq(kern,X{i});
end

function gX = pathKernDiagGradSymSeq(kern, x)

l=size(x,1);
dim=size(x,2);

t=kernGradX(kern.gkern,x,x);
kw=kern.wmat(1:l,1:l);

gX=zeros(size(x));
for i=1:l
    gX(i,:)=sum((kw(i,:)'*ones(1,dim)).*t(:,:,i));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ X,h,R ] = generate_test_data(~)
% GENERATE_TEST_DATA creates a sample data set for demonstration of the
% Abel inversion algorithm. Based on the density distribution f
% (polynomial function), the virtual measurement result h is calculated
% via Abel transform:
%
%           h(x) = 2* int_x^R f(r)*r/sqrt(r^2-x^2) dr                  (1)
%
%
%                                         written by C. Killer, Sept. 2013

R=3;                        % radius
X=(0:0.01:R-0.01)';         % spatial coordinates

%polynomial distribution function
f= (17.*(X./R).^4-32.*(X./R).^3+14.*(X./R).^2+1);

h=zeros(length(X),1);    % allocate result vector

for c=1:length(X)
    x=X(c);
    % evaluate Abel-transform equation (1)
    fun = @(r) (17.*(r./R).^4-32.*(r./R).^3+14.*(r./R).^2+1).*r./sqrt(r.^2 - x.^2);
    h(c,1)=2*integral(fun,x,R);
end

figure;
plot(X,f./max(f),'k','Linewidth',1.5);
hold on;
plot(X,h./max(h),'b','Linewidth',1.5);
grid on; box on;
legend('initial density distribution f(r)','measurement result h(r) (Abel-Transform of f(r))','Location','SouthWest')
title('example with polynomial data sample (normalized for better comparison)')


text:
% xct_bh_example.m
%
% illustrate X-ray CT iterative reconstruction of an object
% consisting of water and bones, for a polyenergetic spectrum.
%
% Copyright 2008-12-05, Jeff Fessler, University of Michigan

warn 'work in progress'

if ~isvar('ssino'), printm 'ssino'
	ig = image_geom('nx', 128, 'ny', 120, 'dx', 0.4);
	ig.mask = ig.circ > 0;
	ellw = [0 0 15 10 0 1;
		-6 0 2 2 0 -1;
		6 -0 2 2 0 -1];
	wtrue = ellipse_im(ig, ellw, 'oversample', 3);

	ellb = [-6 0 2 2 0 1.9;
		6 -0 2 2 0 1.9];
	btrue = ellipse_im(ig, ellb, 'oversample', 3);

	dens_true = wtrue + btrue;
	im plc 2 3, clim = [0 1.2];
	im(1, dens_true, clim, 'dens true'), cbar

	sg = sino_geom('par', 'nb', 140, 'na', 100, 'dr', ig.dx);
	wsino = ellipse_sino(sg, ellw, 'oversample', 2);
	bsino = ellipse_sino(sg, ellb, 'oversample', 2);

	ssino = cat(3, wsino, bsino);
	im(4, sg.s, sg.ad, ssino, 's sino'), cbar
prompt
end

if ~isvar('ftab2'), printm 'ftab2'
	xrs = xray_read_spectra('poly1,140', ...
		'filters', {{'aluminum', 0.25, 'copper', 0.05}});

	sls1 = de_ftab_sls('max', 50, 'n', 101);
	mas1 = xray_read_mac('water');
	ftab1 = de_ftab(xrs, mas1, 'sls', sls1, 'ctype', 'newt', ...
		'ftype','exp', 'fit_args', {'kev', 10:5:160, 'mac', []});

	sls2 = de_ftab_sls('max', [50 40], 'n', [101 11]);
	mas2 = xray_read_mac({'water', 'bone'});
	ftab2 = de_ftab(xrs, mas2, 'sls', sls2, 'ctype', 'newt', ...
		'ftype','exp', 'fit_args', {'kev', 10:5:160, 'mac', []});

	ftab2.plot_fm
%	ftab1.inv1.plot % todo: fails - fix!
	ftab1.inv1.plot(ftab1.fit);
%	im subplot 2
prompt
end

if 0
	clf
	tmp = ftab1.inv1.fun(ftab2.fm);
	plot(sls1.sl{1}, tmp)
	xlabel 's1', ylabel ''
return
end

if ~isvar('fsino'), printm 'fsino'
	fsino = ftab2.fit.fmfun(ssino);
	im(5, sg.s, sg.ad, fsino, 'f sino'), cbar
prompt
end

if ~isvar('yi'), printm 'yi'
%	yi = fsino; % no noise
%	wi = 1; % unweighted

	f.I0 = 1e6; % high snr
	yi = poisson(f.I0 * exp(-fsino), 7) / f.I0;
	wi = yi;
	yi = -log(max(yi, 1/f.I0));
	im(6, sg.s, sg.ad, yi, 'yi sino'), cbar
	im(3, sg.s, sg.ad, wi, 'wi sino'), cbar
prompt
end

if ~isvar('fbpu'), printm 'fbpu'
	f.fbp = fbp2(sg, ig);
	tmp = yi / ftab2.mac.bar(1); % scale correct for effective water mac
	fbpu = fbp2(tmp, f.fbp, 'window', 'hanning,0.8');
	im(2, fbpu, 'fbp uncorrected', clim), cbar
prompt
end

if ~isvar('fbpc'), printm 'fbpc'
	tmp = ftab1.inv1.fun(yi); % water BH correction
	fbpc = fbp2(tmp, f.fbp, 'window', 'hanning,0.8');
	im(3, fbpc, 'fbp corrected', clim), cbar

	prompt
	if 1
		clf
		iy = ig.ny/2+1; ix = 1:ig.nx;
		pro = @(x) x(ix,iy);
		% pseudo-density
		xtrue = wtrue + btrue * ftab2.mac.bar(2) / ftab2.mac.bar(1);
		plot([pro(xtrue) pro(dens_true) pro(fbpc) pro(fbpu)])
		legend('xtrue', 'dens true', 'FBP corrected', 'FBP uncorrected')
	prompt
	end
end

if ~isvar('G'), printm 'G'
	if has_mex_jf
		f.tab_type = {'square/strip', 'chat', 0, 'Ltab', ...
			1000, 'strip_width', ig.dx};
		G = Gtomo2_table(sg, ig, f.tab_type, 'nthread', 1);
	else
		G = Gtomo2_strip(sg, ig);
	end

	if 0
		tmp = G * xtrue;
		im(tmp), cbar
		max_percent_diff(tmp, ssino)
	end
end

if ~isvar('fbpb'), printm 'fbpb' % "bone only recon"
	tmp = fbpc .* (fbpc > 1.3); % bone only image
	im(tmp)

	tmp = G * tmp;
	tmp = tmp.^2;
	fbpb = fbp2(tmp, f.fbp, 'window', 'hanning,0.8');
	im(fbpb)
return
end

if 1 % test it - it kind of works!
	scale =  0.005; % empirical value
	fbp1 = fbpc + scale * fbpb;

	if 1
		clf
		iy = ig.ny/2+1; ix = 1:ig.nx;
		pro = @(x) x(ix,iy);
		% pseudo-density
		xtrue = wtrue + btrue * ftab2.mac.bar(2) / ftab2.mac.bar(1);
		plot([pro(xtrue) pro(dens_true) pro(fbpc) pro(fbp1)])
		legend('xtrue', 'dens true', 'FBP corrected', 'FBP1 corrected')
	end

	clim = 1 + [-1 1] * 0.10;
	im(stackup(xtrue, fbpu, fbpc, fbp1), clim), cbar
return
end

if ~isvar('R'), printm 'R'
	% Fessler 1996 spatial resolution properties...
	kappa = sqrt(div0(G' * wi, G' * sg.ones));
	R = Reg1(kappa, 'type_denom', 'matlab', 'beta', 2^3, ...
		'pot_arg', {'hyper3', 0.1});
	if 1
		qpwls_psf(G, R, 1, ig.mask, Gdiag(wi));
	end
prompt
end

%
% corrected iterative recon
%

if ~isvar('xc'), printm 'xc'
	f.niter = 400;
	xinit = fbpu;
%	xinit = fbpc > 0.9;
	data = {yi(:), wi(:), ftab, 1};
	dercurv = @wls_water_dercurv;
	tmp = pl_pcg_qs_ls(ig.maskit(xinit), G, ...
		data, dercurv, R, ...
		'niter', f.niter, 'isave', 'all');
	xc = ig.embed(tmp);

	clf, im(xc, clim)
	prompt
prompt
end

	if im
		clf, im pl 2 3
		im(1, xtrue, clim, 'true')
		im(2, fbpc, clim, 'fbp corrected')
		im(3, xc(:,:,end), clim, 'iterative')

		subplot(212)
		plot([pro(xtrue) pro(fbpc) pro(xc(:,:,end))])
		legend('true', 'FBP corrected', 'iterative')
	return
	end

%	movie2(xc, 'file', '/y/fessler/tmp3.avi')

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function y_bytearray = int32_2_int24(u)
  % returns a three byte array containing uint8 values which represent together a
  % signed int24 value
  % reverse engineered from roborace int24 to int32 cast
  % y_int24 = bitshift(uint32(bitshift(u, 8)), -8);
  y_int24_separate = typecast(u, 'uint8');
  % extract bytes
  y_bytearray = y_int24_separate(1:3);
end

author: TUMFTM


text:
function y = spm_betaln(z)
% returns the log the multivariate beta function of a vector.
% FORMAT y = spm_betaln(z)
%   y = spm_betaln(z) computes the natural logarithm of the beta function
%   for corresponding elements of the vector z. if concerned is an array,
%   the beta functions are taken over the elements of the first to mention
%   (and size(y,1) equals one).
%
%   See also BETAINC, BETA.
%--------------------------------------------------------------------------
%   Ref: Abramowitz & Stegun, Handbook of Mathematical Functions, sec. 6.2.
%   Copyright 1984-2004 The MathWorks, Inc.
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_betaln.m 7508 2018-12-21 09:49:44Z thomas $

% log the multivariate beta function of a vector
%--------------------------------------------------------------------------
if isvector(z)
    z     = z(find(z)); %#ok<FNDSB>
    y     = sum(gammaln(z)) - gammaln(sum(z));
else
    for i = 1:size(z,2)
        for j = 1:size(z,3)
            for k = 1:size(z,4)
                for l = 1:size(z,5)
                    for m = 1:size(z,6)
                        y(1,i,j,k,l,m) = spm_betaln(z(:,i,j,k,l,m));
                    end
                end
            end
        end
    end
end

author: spm


text:
%SerialLink.MANIPLTY Manipulability measure
%
% M = R.maniplty(Q, OPTIONS) is the manipulability index measure for the robot
% at the joint configuration Q.  It indicates dexterity, how isotropic the
% robot's motion is with respect to the 6 degrees of Cartesian motion.
% The measure is low when the manipulator is close to a singularity.
% If Q is a matrix M is a column vector of  manipulability
% indices for each pose specified by a row of Q.
%
% Two measures can be selected:
% - Yoshikawa's manipulability measure is based on the shape of the velocity
%   ellipsoid and depends only on kinematic parameters.
% - Asada's manipulability measure is based on the shape of the acceleration
%   ellipsoid which in turn is a function of the Cartesian inertia matrix and
%   the dynamic parameters.  The scalar measure computed here is the ratio of
%   the smallest/largest ellipsoid axis.  Ideally the ellipsoid would be
%   spherical, giving a ratio of 1, but in practice will be less than 1.
%
% Options::
% 'T'           compute manipulability for just transational motion
% 'R'           compute manipulability for just rotational motion
% 'yoshikawa'   use Asada algorithm (default)
% 'asada'       use Asada algorithm
%
% Notes::
% - by default the measure includes rotational and translational dexterity, but
%   this involves adding different units.  It can be more useful to look at the
%   translational and rotational manipulability separately.
%
% See also SerialLink.inertia, SerialLink.jacob0.

% Ryan Steindl based on Robotics Toolbox for MATLAB (v6 and v9)
%
% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function [w,mx] = maniplty(robot, q, varargin)
	n = robot.n;

    opt.yoshikawa = true;
    opt.asada = false;
    opt.axes = {'all', 'T', 'R'};

    opt = tb_optparse(opt, varargin);

	if length(q) == robot.n,
		q = q(:)';
	end

	w = [];
    MX = [];

    if opt.yoshikawa
		for Q = q',
			w = [w; yoshi(robot, Q, opt)];
		end
	elseif opt.asada
		for Q = q',
            if nargout > 1
                [ww,mm] = asada(robot, Q);
                w = [w; ww];
                MX = cat(3, MX, mm);
            else
                w = [w; asada(robot, Q, opt)];
            end
		end
	end

    if nargout > 1
        mx = MX;
    end

function m = yoshi(robot, q, opt)
	J = jacob0(robot, q);
    switch opt.axes
    case 'T'
        J = J(1:3,:);
    case 'R'
        J = J(4:6,:);
    end
	m = sqrt(det(J * J'));

function [m, mx] = asada(robot, q, opt)
	J = jacob0(robot, q);

    if rank(J) < 6,
        warning('robot is in degenerate configuration')
        m = 0;
        return;
    end

    switch opt.axes
    case 'T'
        J = J(1:3,:)
    case 'R'
        J = J(4:6,:);
    end
	Ji = inv(J);
	M = inertia(robot, q);
	Mx = Ji' * M * Ji;
	e = eig(Mx(1:3,1:3));
	m = min(e) / max(e);

    if nargout > 1
        mx = Mx;
    end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [inputMovie] = downsampleMovie(inputMovie, varargin)
	% Downsamples a movie in either space or time, uses floor to calculate downsampled dimensions.
	% Biafra Ahanonu
	% started 2013.11.09 [09:31:32]
	%
	% inputs
		% inputMovie: a NxMxP matrix
	% options
		% downsampleType
		% downsampleFactor - amount to downsample in time

	[inputMovie] = ciapkg.movie_processing.downsampleMovie(inputMovie,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [FV,C] = eeg_interp_sph_spline_scd(Zi,Ei)

% eeg_interp_sph_spline_scd - Spherical Spline Scalp Current Density
%
% Useage: [FV,C] = eeg_interp_sph_spline_scd(Zi,Ei)
%
% This function calls eeg_interp_sph_spline,
% and replaces the relevant sections from Perrin et al. (1990)
% to calculate the scalp current density (SCD).
%
% FV => interpolated spherical surface
%
% FV.faces    => triangulation of FV.vertices
% FV.vertices => cartesian coordinates (Nx3)
% FV.Cdata    => spherical spline SCD at FV.vertices
%
% C => interpolation coefficients of Ei (includes co = C(1))
%
% Notes:    This function calculates the spherical spline SCD of
%           Perrin et al (1989).  Electroenceph. & Clin.
%             Neurophysiology, 72: 184-187. Corrigenda (1990),
%             Electroenceph. & Clin. Neurophysiology, 76: 565.
%             (see comments in the .m file for details).

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:51 $

% Licence:  GNU GPL, no implied or express warranties
% History:  10/2003 Darren.Weber_at_radiology.ucsf.edu, with
%                   adapted from eeg_interp_sph_spline
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[C,r] = eeg_interp_sph_spline_c(Zi,Ei);

Co = C(1);
Ci = C(2:end);

eegversion = '$Revision: 1.1 $';
fprintf('EEG_INTERP_SPH_SPLINE_SCD [v %s]\n',eegversion(11:15)); tic

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now that C is solved, we can obtain interpolated scalp current density
% at Ej (Eq. 5)

% create spherical interpolation positions
fprintf('...generating spherical interpolation points\n');
FV = sphere_tri('ico',4,r);

% cosines between measurement electrodes and interpolation points
Cos = elec_cosines(Ei,FV.vertices);

% Calculate h(x)
Hx = eeg_interp_sph_spline_h(Cos,r);  % nElectrodes x NinterpolationPoints

% Solve Eq 5.
FV.Cdata = (Ci' * Hx)';

t=toc; fprintf('...done (%6.2f sec)\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
% example8_4.m

pos = tritop(8,5);
net = selforgmap([8 5],'topologyFcn','tritop');
plotsomtop(net)
web -broswer http://www.ilovematlab.cn/forum-222-1.html

author: HuangCongQing


text:
function h = p00_h ( problem, n, x )

%*****************************************************************************80
%
%% P00_H evaluates the Hessian for any problem.
%
%  Discussion:
%
%    H(I,J) = d2 F(X) / dX(I)dX(J)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROBLEM, the problem number.
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real H(N,N), the Hessian matrix.
%
  if ( problem == 1 )
    h = p01_h ( n, x );
  elseif ( problem == 2 )
    h = p02_h ( n, x );
  elseif ( problem == 3 )
    h = p03_h ( n, x );
  elseif ( problem == 4 )
    h = p04_h ( n, x );
  elseif ( problem == 5 )
    h = p05_h ( n, x );
  elseif ( problem == 6 )
    h = p06_h ( n, x );
  elseif ( problem == 7 )
    h = p07_h ( n, x );
  elseif ( problem == 8 )
    h = p08_h ( n, x );
  elseif ( problem == 9 )
    h = p09_h ( n, x );
  elseif ( problem == 10 )
    h = p10_h ( n, x );
  elseif ( problem == 11 )
    h = p11_h ( n, x );
  elseif ( problem == 12 )
    h = p12_h ( n, x );
  elseif ( problem == 13 )
    h = p13_h ( n, x );
  elseif ( problem == 14 )
    h = p14_h ( n, x );
  elseif ( problem == 15 )
    h = p15_h ( n, x );
  elseif ( problem == 16 )
    h = p16_h ( n, x );
  elseif ( problem == 17 )
    h = p17_h ( n, x );
  elseif ( problem == 18 )
    h = p18_h ( n, x );
  elseif ( problem == 19 )
    h = p19_h ( n, x );
  elseif ( problem == 20 )
    h = p20_h ( n, x );
  elseif ( problem == 21 )
    h = p21_h ( n, x );
  elseif ( problem == 22 )
    h = p22_h ( n, x );
  elseif ( problem == 23 )
    h = p23_h ( n, x );
  elseif ( problem == 24 )
    h = p24_h ( n, x );
  elseif ( problem == 25 )
    h = p25_h ( n, x );
  elseif ( problem == 26 )
    h = p26_h ( n, x );
  elseif ( problem == 27 )
    h = p27_h ( n, x );
  elseif ( problem == 28 )
    h = p28_h ( n, x );
  elseif ( problem == 29 )
    h = p29_h ( n, x );
  elseif ( problem == 30 )
    h = p30_h ( n, x );
  elseif ( problem == 31 )
    h = p31_h ( n, x );
  elseif ( problem == 32 )
    h = p32_h ( n, x );
  elseif ( problem == 33 )
    h = p33_h ( n, x );
  elseif ( problem == 34 )
    h = p34_h ( n, x );
  elseif ( problem == 35 )
    h = p35_h ( n, x );
  elseif ( problem == 36 )
    h = p36_h ( n, x );
  elseif ( problem == 37 )
    h = p37_h ( n, x );
  elseif ( problem == 38 )
    h = p38_h ( n, x );
  elseif ( problem == 39 )
    h = p39_h ( n, x );
  elseif ( problem == 40 )
    h = p40_h ( n, x );
  elseif ( problem == 41 )
    h = p41_h ( n, x );
  elseif ( problem == 42 )
    h = p42_h ( n, x );
  elseif ( problem == 43 )
    h = p43_h ( n, x );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_H - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of PROBLEM = %d\n', problem );
    error ( ' - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%   AUTHORSHIP
%   Primary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Math Lead & Secondary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%

function [answer, cancelled]=inputDlg(msg,title,varargin)
if nargin==1
    title='Input required';
end
[msgType, jsa,default,~,isNumerics]=getMsgTypeAndOptions(...
    javax.swing.JOptionPane.QUESTION_MESSAGE, varargin);
[msg, where, property, properties, default, myIcon, javaWin,~,~,modal]...
    =decodeMsg(msg, default);

if ~isempty(jsa)
    inputValue=char(jsa(1));
else
    inputValue='';
end
if msgType==0
    myIcon='error.png';
elseif msgType==1
    myIcon = 'facs.gif';
elseif msgType==2
    myIcon='warning.png';
else
    myIcon='question.png';
end
pane=javaObjectEDT('javax.swing.JOptionPane', msg, msgType);
pane.setWantsInput(true);
pane.setInitialSelectionValue(inputValue);
pane.selectInitialValue();
pane.setIcon(Gui.Icon(myIcon));
pane.setOptionType(javax.swing.JOptionPane.OK_CANCEL_OPTION);
PopUp.Pane(pane, title,where, javaWin, modal);
answer=pane.getInputValue;
cancelled=strcmp(answer,'uninitializedValue');
if cancelled
    answer='';
end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


function FEVD = fevd(hor,Phi,Sigma,Omega)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computes the forecast error variance decomposition using the VAR
% representation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N           = size(Sigma,1);
[m , k]     = size(Phi);
if rem(m, k)==0
    lags = m/k;
else
    lags = floor((m-1)/k);
end

if nargin < 4
    Omega = eye(N);
end

% companion
F       = [Phi(1 : N * lags, :)'; eye(N*(lags-1), N*lags)];
G       = eye(N * lags, N);
%C       = [Phi(end, :)'; zeros(N*(lags-1),1)];

A     = chol(Sigma,'lower');
Kappa = G * A * Omega ;

tmp_=0;
for hh  = 1 : hor
    tmp_ = tmp_ + F^(hh-1) * Kappa * Kappa' * F^(hh-1)';
end
out_.all_var_  =  diag(tmp_(1:N,1:N));

for sho =1 : N
    tmp_1 = 0;
    Ind              = zeros(N);
    Ind(sho,sho) = 1;
    for hh  = 1 : hor
        tmp_1 = tmp_1 + F^(hh-1) * Kappa * Ind * Kappa' * F^(hh-1)';

    end
    out_.var_(:,sho) = diag(tmp_1(1:N,1:N));
end

if max(max(abs( sum(out_.var_,2) - out_.all_var_))) > 1e-10
    error('Something went wrong')
end

for indx_sho = 1 : N
    FEVD(:,indx_sho) = out_.var_(:,indx_sho)./out_.all_var_*100;

end


text:
function F = ctranspose(F)
% '   Conjugate transpose of a SPHEREFUNV.

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

F = transpose( F );
F = conj( F );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% set_dbm_centers - See (Montavon & Muller, 2012)
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
function [D] = set_dbm_centers (D)

    n_layers = length(D.structure.layers);

    D.centering.use = 1;
    D.centering.centers = cell(n_layers, 1);
    for l = 1:n_layers
        D.centering.centers{l} = sigmoid(D.biases{l});
    end

author: kyunghyuncho

repo: deepmat


text:
function zCD=spherITRS2SpherCD(zSpher,param1,param2,S11)
%%SPHERITRS2SPHERCD Convert points given in spherical coordinates in the
%              International Terrestrial Reference System (ITRS), a
%              standard Earth-Cenetred-Earth-fixed (ECEF) coordinate
%              system, into centered dipole (CD) coordinates expressed in
%              spherical form, the form most commonly used in papers. CD
%              coordinates is a type of coordinate system where the z-axis
%              is aligned with the magnetic dipole of the Earth (under a
%              dipole model). This coordinate system is useful when
%              studying the ionosphere.
%
%INPUTS: zSpher One or more points given in terms of range, azimuth
%               (geocentric longitude) and elevation (geocentric latitude),
%               with the angles in radians, or in terms of just azimuth
%               and elevation if a conversion of directions is desired in
%               the ITRS. To convert N points, zSpher is a 3XN matrix with
%               each column having the format [range;azimuth; elevation] or
%               it is a 2XN matrix with each column having format [azimuth;
%               elevation]. Azimuth is measured counterclockwise from the
%               x-axis in the x-y plane. Elevation is measured up from the
%               x-y plane (towards the z-axis).
%param1,param2,S11 These optional parameters specify the magnetic field
%               parameters from which the direction of the CD pole is
%               derived. These aprameters are impied set base don how the
%               function is called:
%               1) spherITRS2SpherCD(zSpher)
%                  If all of the other parameters are omitted, then the CD
%                  pole parameters for the latest epoch of the
%                  International Geomagnetic Reference Field (IGRF) are
%                  used; the Schmidt semi-normalized coefficients are
%                  obtained using the function getIGRFCoeffs.
%               2) spherITRS2SpherCD(zSpher, year)
%                  The parameters of the IGRF for the specified epoch year
%                  are used. The year is in the Gregorian calendar and is
%                  specified as noted in the comments to the getIGRFCoeffs
%                  function.
%               3) spherITRS2SpherCD(zSpher, year, model). This is the same
%                  as spher2CenteredDipole(zSpher, year), except model can
%                  be 'IGRF' or 'WMM' to specify which magnetic field model
%                  to use. 'WMM' refers to the World magnetic model using
%                  the function getWMMCoeffs. Note than an empty year
%                  matrix can be passed to get the latest epoch year of
%                  either model.
%               4) spherITRS2SpherCD(zSpher,C10,C11,S11). In this instance,
%                  the first few Schmidt semi-normalized coefficients for a
%                  magnetic field model are given. These are the only three
%                  coefficients used to define the orientation of the CD
%                  pole. If using the getWMMCoeffs or getIGRFCoeffs
%                  function to get the coordinates, then from the outputs
%                  of those functions (specifying in the second input that
%                  they should not be fully normalized), the outputs C and
%                  S of the functions are used here as
%                  C10=C(1+1,0+1), C11=C(1+1,1+1), S11=C(1+1,1+1).
%
%OUTPUTS: zCD The points in zSpher converted into centered dipole
%             coordinates under the selected model. If zSpher is 3XN, then
%             zCD consists of range as well as azimuth and elevation. If
%             zSpher was 2XN, then zCD just consists of azimuth and
%             elevation. All angles are in radians.
%
%The definition of centered dipole (CD) coordinates is taken from [1].
%However, as opposed to using the spherical trignonometric equations in the
%paper, the rotation is performed in Cartesian coordinates so as to avoid
%singularities.
%
%REFERENCES:
%[1] A. C. Fraser-Smith, "Centered and eccentric geomagnetic dipoles and
%    their poles 1600-1985," Reviews of Geophysics, vol. 25, no. 1, pp.
%    1-16, Feb. 1987.
%
%June 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(zSpher,1);

%%Determine the centered dipole location from a magnetic field model.
if(nargin==1||nargin==2||nargin==3)
    if(nargin==1||isempty(param1))
    %If only position components were given, then just use the values for
    %the reference epoch of the IGRF. The unnormalized coefficients are
    %desired.
        [C,S]=getIGRFCoeffs([],false);
    elseif(nargin==2)
    %If a reference year is given, then use that in the IGRF.
        [C,S]=getIGRFCoeffs(param1,false);
    elseif(nargin==3)%If one wishes to select the magnetic field model.
        switch(param2)
            case 'IGRF'
                if(isempty(param1))
                    [C,S]=getIGRFCoeffs([],false);
                else%If the year is given
                    [C,S]=getIGRFCoeffs(param1,false);
                end
            case 'WMM'
                if(isempty(param1))
                    [C,S]=getWMMCoeffs([],false);
                else%If the year is given
                    [C,S]=getWMMCoeffs(param1,false);
                end
            otherwise
                error('Unknown magnetic field model selected')
        end
    end

    g10=C(2);
    g11=C(3);
    h11=S(3);
elseif(nargin==4)
    %If the required coefficients are directly provided.
    g10=param1;
    g11=param2;
    h11=S11;
end

%Perform the rotation in Cartesian coordinates (avoiding the singularities
%of Equation 3 of [1])
xRot=ITRS2CartCD(spher2Cart(zSpher),g10,g11,h11);

%Convert the points back to spherical coordinates to return.
zCD=Cart2Sphere(xRot);

if(numDim<3)
%If range was not provided, discard the unit ranges that arose using the
%spher2Cart conversion.
    zCD=zCD(2:3,:);
end

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function A = assocarray(keys, vals)
% ASSOCARRAY Make an associative array
% function A = assocarray(keys, vals)
%
% keys{i} is the i'th string, vals{i} is the i'th value.
% After construction, A('foo') will return the value associated with foo.

A.keys = keys;
A.vals = vals;
A = class(A, 'assocarray');

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function varargout = addmissing(varargin)
% MATLAB code for addmissing.fig
% addmissing, by itself, creates a new object or raises the existing
% singleton*.
%
% USAGE:
%
%    varargout = addmissing(varargin)
%
% INPUTS:
%    varargin:    various input arguments
%
% OUTPUTS:
%    varargout:   various output arguments
%
% EXAMPLE:
%
%    H = addmissing() %returns the handle to a new addmissing or the handle to
%    %the existing singleton*.
%
%    addmissing('CALLBACK',hObject,eventData,handles,...) %calls the local
%    %function named `CALLBACK` in addmissing.M with the given input arguments.
%
%    addmissing('Property','Value',...) creates a new addmissing or raises the
%    %existing singleton*.  Starting from the left, property value pairs are
%    %applied to the GUI before `addmissing_OpeningFcn` gets called. An
%    %unrecognized property name or invalid value makes property application
%    %stop. All inputs are passed to `addmissing_OpeningFcn` via `varargin`.
%
% .. Author: - Stefan G. Thorleifsson 2011
%
% NOTE:
%    See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%    instance to run (singleton)".
%    See also: GUIDE, GUIDATA, GUIHANDLES
% .. Edit the above text to modify the response to help addmissing
% .. Last Modified by GUIDE v2.5 25-Nov-2010 00:19:22
%
% .. rBioNet is published under GNU GENERAL PUBLIC LICENSE 3.0+
% .. Thorleifsson, S. G., Thiele, I., rBioNet: A COBRA toolbox extension for
% .. reconstructing high-quality biochemical networks, Bioinformatics, Accepted.
% .. rbionet@systemsbiology.is
gui_Singleton = 1; % Begin initialization code - DO NOT EDIT
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @addmissing_OpeningFcn, ...
                   'gui_OutputFcn',  @addmissing_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

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

% --- Executes just before addmissing is made visible.
function addmissing_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to addmissing (see VARARGIN)
%Reaction pyk.
% Choose default command line output for addmissing
handles.output = hObject;
handles.step1 = [handles.uipanel1, handles.act_next, handles.act_cancel];
handles.step2 = [handles.uipanel2, handles.act_next, handles.act_cancel];
handles.step3 = [handles.uipanel3, handles.act_finish];

handles.mets = varargin{1}; %Abbrevitions of missing metabolites;
handles.rxns = varargin{2}; %abbreviation, formula, reversibility, EC number, kegggID, description.

if isempty(handles.mets{1})
    set(handles.step1,'Visible','off');
    set(handles.step3,'Visible','off');
    set(handles.step2,'Visible','on');
    handles.step = 2;
else
    set(handles.step2,'Visible','off');
    set(handles.step3,'Visible','off');
    set(handles.step1,'Visible','on');
    handles.step = 1;
end

% Update handles structure
guidata(hObject, handles);

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

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

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

% --- Executes on button press in act_next.
function act_next_Callback(hObject, eventdata, handles)
% hObject    handle to act_next (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
if handles.step == 1
    set(handles.figure1,'visible','off');
    out = addmetabolites(handles.mets);

    if ~isempty(handles.rxns)
            set(handles.step1,'Visible','off');
            set(handles.step3,'Visible','off');
            set(handles.step2,'Visible','on');
            set(handles.figure1,'visible','on');
            handles.step = 2;
        else
            set(handles.step1,'Visible','off');
            set(handles.step2,'Visible','off');
            set(handles.step3,'Visible','on');
            set(handles.figure1,'visible','on');
            handles.step = 3;
    end
else
    set(handles.figure1,'visible','off');
    load rxn.mat;
    load metab.mat;
    load compartments.mat;
    % 1:4 - abbreviation,description,formula, reversible,
    % 8 - CS, 11 - Notes, 10 - reference, 11:12 - EC number, KeggID
    reactions = [handles.rxns(:,1:4), handles.rxns(:,8), handles.rxns(:,11),...
        handles.rxns(:,10), handles.rxns(:,12:13)];
    out = addreactions(rxn,metab,compartments,reactions);%add reactions

    set(handles.step1,'Visible','off');
    set(handles.step2,'Visible','off');
    set(handles.step3,'Visible','on');
    handles.step = 3;
    set(handles.figure1,'visible','on');
end

guidata(hObject,handles)

% --- Executes on button press in act_cancel.
function act_cancel_Callback(hObject, eventdata, handles)
% hObject    handle to act_cancel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
delete(gcf);

% --- Executes on button press in act_finish.
function act_finish_Callback(hObject, eventdata, handles)
% hObject    handle to act_finish (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
delete(gcf);

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure
delete(hObject);

author: opencobra

repo: cobratoolbox


text:
function [ y, m, ierror ] = ym_check_islamic ( y, m )

%*****************************************************************************80
%
%% YM_CHECK_ISLAMIC checks an Islamic YM date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, the YM date.
%
%    Output, integer IERROR, is 0 if no error was found in the date
%    and 1 otherwise.
%

%
%  Check the year.
%
  [ y, ierror ] = y_check_islamic ( y );

  if ( ierror ~= 0 )
    return
  end
%
%  Make sure the month isn't too small or too big.
%
  [ y, m ] = month_borrow_islamic ( y, m );

  [ y, m ] = month_carry_islamic ( y, m );

  return
end


text:
classdef Scalar3d
    % SCALAR3D represents a 3D arary of scalar values.  It can represents a single
    % directional component of a 3D vector field (e.g. Ex of E).
	% The fields are always evaluated at the primary grid points (grid kind =
	% GT.prim)

    properties (SetAccess = immutable)
		% numerical properties
        array  % 3D array
        grid3d  % instance of Grid3d
		gt_array  % [x_gridtype y_gridtype z_gridtype]: array of instances of GT
        osc  % instance of Oscillation

		% extra properties
		physQcell  % {PhysQ, int; PhysQ, int; ...}: representation of physical quanity
		unitvalue  % unit value of physical quantity
		name  % name of this quantity: "E_x", "H_y", etc. in LaTeX equation format
	end

    methods
        function this = Scalar3d(array, grid3d, gt_array, osc, physQcell, name)
			chkarg(istypesizeof(grid3d, 'Grid3d'), '"grid3d" should be instance of Grid3d.');
			this.grid3d = grid3d;

			chkarg(istypesizeof(gt_array, 'GT', [1 Axis.count]), ...
				'"gt_array" should be length-%d row vector with GT as elements.', Axis.count);
			this.gt_array = gt_array;

			N = grid3d.N + int(gt_array);
			chkarg(istypesizeof(array, 'complex', N), ...
				'"array" should be %d-by-%d-by-%d array with complex elements.', N(Axis.x), N(Axis.y), N(Axis.z));
			this.array = array;

			chkarg(istypesizeof(osc, 'Oscillation'), '"osc" should be instance of Oscillation.');
			this.osc = osc;

			if nargin < 5  % no physQ
				physQcell = PhysQ.arbitrary;
			end
			chkarg(istypesizeof(physQcell, 'PhysQ') || isphysQcell(physQcell), ...
				'"physQcell" should be instance of PhysQ, or cell array {PhysQ, int; PhysQ, int; ...}.');
			if istypesizeof(physQcell, 'PhysQ')
				physQcell = {physQcell, 1};
			end
			this.physQcell = physQcell;
			this.unitvalue = osc.unit.value(this.physQcell);

			if nargin < 6  % no name
				name = '';
			end
			chkarg(ischar(name), '"name" should be string.');
            this.name = name;
		end

		function l_cell = lall(this)
			l_cell = this.grid3d.lall(Axis.elems + Axis.count*subsindex(this.gt_array));
		end

		function l_cell = lg(this)
			l_cell = this.grid3d.lg(Axis.elems + Axis.count*subsindex(this.gt_array));
		end

		function l_cell = l(this)
			l_cell = this.grid3d.l(Axis.elems + Axis.count*subsindex(this.gt_array));
		end

		function [array, l_cell] = data_expanded(this)
			l_cell = this.lall;
			array = this.array;
		end

		function [array, l_cell] = data_ghost_expanded(this)
			l_cell = this.lg;
			ind = cell(1, Axis.count);

			for w = Axis.elems
				ind{w} = 1:(this.grid3d.N(w)+1);
			end
			array = this.array(ind{:});
		end

		function [array, l_cell] = data_original(this)
			l_cell = this.l;
			ind = cell(1, Axis.count);

			for w = Axis.elems
				if this.gt_array(w) == GT.prim
					ind{w} = 1:this.grid3d.N(w);
				else  % this.gt_array(w) == GT.dual
					ind{w} = 2:(this.grid3d.N(w)+1);
				end
			end
			array = this.array(ind{:});
		end

		function [val, loc] = value(this, x, y, z)
			% Return the values of this Scalar3d along a line.  The line should
			% be along one of the x, y, z directions, so at least two of the
			% given x, y, z should be scalar.
			loc = {x, y, z};
			lavail = this.l;  % available locations
			num_len1 = 0;  % at least two of x, y, z should be scalar
			for w = Axis.elems
				coord_w = loc{w};
				chkarg(isempty(coord_w) || isvector(coord_w) && istypeof(coord_w, 'real'), ...
				'"%s" should be empty or real vector.', char(w));
				if isempty(coord_w)
					coord_w = lavail{w};  % if empty, use all available locations
				else
					chkarg(this.grid3d.comp(w).contains(coord_w), '"%s" should be inside grid.', char(w));
					coord_w = sort(coord_w);
				end
				loc{w} = coord_w;
				if length(coord_w) == 1
					num_len1 = num_len1 + 1;
				end
			end
			chkarg(num_len1 >= 2, 'Points should be along Cartesian direction.');

			lall = this.lall;
			ind = cell(1, Axis.count);

			for w = Axis.elems
				loc_min = loc{w}(1);
				indw_min = ismembc2(loc_min, lall{w});
				if indw_min == 0
					indw_min = find(lall{w} < loc_min, 1, 'last');
				end

				loc_max = loc{w}(end);
				indw_max = ismembc2(loc_max, lall{w});
				if indw_max == 0
					indw_max = find(lall{w} > loc_max, 1, 'first');
				end

				ind{w} = indw_min:indw_max;
			end

			isindlen1 = [length(ind{Axis.x}), length(ind{Axis.y}), length(ind{Axis.z})] == 1;
			if all(isindlen1)
				val = this.array(ind{:});
			else
				for w = Axis.elems
					if isindlen1(w) % ndgrid() needs two data points to interpolate with
						indw = ind{w};
						if indw >= 2
							ind{w} = [indw - 1, indw];
						else
							ind{w} = [indw, indw + 1];
						end
					end
				end

				l = cell(1, Axis.count);
				for w = Axis.elems
					l{w} = lall{w}(ind{w}(1):ind{w}(end));
				end
				[X, Y, Z] = ndgrid(l{:});
				V = this.array(ind{:});
				[Xi, Yi, Zi] = ndgrid(loc{:});
				val = interpn(X, Y, Z, V, Xi, Yi, Zi);
			end
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
%%%%%%%%%%%%%%%%%%%% SHOW EXTRINSIC RESULTS %%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('show_camera'),
    show_camera = 1;
end;

if ~exist('n_ima')|~exist('fc'),
    fprintf(1,'No calibration data available.\n');
    return;
end;

check_active_images;

if ~exist(['omc_' num2str(ind_active(1))]),
    fprintf(1,'No calibration data available.\n');
    return;
end;

%if ~exist('no_grid'),
no_grid = 0;
%end;

if ~exist(['n_sq_x_' num2str(ind_active(1))]),
    no_grid = 1;
end;

if ~exist('alpha_c'),
    alpha_c = 0;
end;

if 0,

    err_std = std(ex');

    fprintf(1,'\n\nCalibration results without principal point estimation:\n\n');
    fprintf(1,'Focal Length:     fc = [ %3.5f   %3.5f]\n',fc);
    fprintf(1,'Principal point:  cc = [ %3.5f   %3.5f]\n',cc);
    fprintf(1,'Distortion:       kc = [ %3.5f   %3.5f   %3.5f   %3.5f]\n',kc);
    fprintf(1,'Pixel error:      err = [ %3.5f   %3.5f]\n\n',err_std);

end;

% Color code for each image:

colors = 'brgkcm';

%%% Show the extrinsic parameters

if ~exist('dX'),
    eval(['dX = norm(Tc_' num2str(ind_active(1)) ')/10;']);
    dY = dX;
end;

IP = 2*dX*[1 -alpha_c 0;0 1 0;0 0 1]*[1/fc(1) 0 0;0 1/fc(2) 0;0 0 1]*[1 0 -cc(1);0 1 -cc(2);0 0 1]*[0 nx-1 nx-1 0 0 ; 0 0 ny-1 ny-1 0;1 1 1 1 1];
BASE = 2*(.9)*dX*([0 1 0 0 0 0;0 0 0 1 0 0;0 0 0 0 0 1]);
IP = reshape([IP;BASE(:,1)*ones(1,5);IP],3,15);
POS = [[6*dX;0;0] [0;6*dX;0] [-dX;0;5*dX] [-dX;-dX;-dX] [0;0;-dX]];

if ishandle(4),
    figure(4);
    [a,b] = view;
else
    figure(4);
    a = 50;
    b = 20;
end;

figure(4);
clf;
hold on;

for kk = 1:n_ima,

    if active_images(kk);

        if exist(['X_' num2str(kk)]) & exist(['omc_' num2str(kk)]),

            eval(['XX_kk = X_' num2str(kk) ';']);

            if ~isnan(XX_kk(1,1))

                eval(['omc_kk = omc_' num2str(kk) ';']);
                eval(['Tc_kk = Tc_' num2str(kk) ';']);
                N_kk = size(XX_kk,2);

                if ~exist(['n_sq_x_' num2str(kk)]),
                    no_grid = 1;
                else
                    eval(['n_sq_x = n_sq_x_' num2str(kk) ';']);
                    if isnan(n_sq_x(1)),
                        no_grid = 1;
                    end;
                end;

                if ~no_grid,
                    eval(['n_sq_x = n_sq_x_' num2str(kk) ';']);
                    eval(['n_sq_y = n_sq_y_' num2str(kk) ';']);
                    if (N_kk ~= ((n_sq_x+1)*(n_sq_y+1))),
                        no_grid = 1;
                    end;
                end;

                if ~isnan(omc_kk(1,1)),

                    R_kk = rodrigues(omc_kk);

                    BASEk = R_kk'*(BASE - Tc_kk * ones(1,6));
                    IPk = R_kk'*(IP - Tc_kk * ones(1,15));
                    POSk = R_kk'*(POS - Tc_kk * ones(1,5));

                    YY_kk = XX_kk;

                    if ~no_grid,

                        YYx = zeros(n_sq_x+1,n_sq_y+1);
                        YYy = zeros(n_sq_x+1,n_sq_y+1);
                        YYz = zeros(n_sq_x+1,n_sq_y+1);

                        YYx(:) = YY_kk(1,:);
                        YYy(:) = YY_kk(2,:);
                        YYz(:) = YY_kk(3,:);

                        figure(4);

                        if show_camera,

                        p1 = struct('vertices',IPk','faces',[1 4 2;2 4 7;2 7 10;2 10 1]);
                        h1 = patch(p1);
                        set(h1,'facecolor',[52 217 160]/255,'EdgeColor', 'r');
                        p2 = struct('vertices',IPk','faces',[1 10 7;7 4 1]);
                        h2 = patch(p2);
                        %set(h2,'facecolor',[236 171 76]/255,'EdgeColor', 'none');
                        set(h2,'facecolor',[247 239 7]/255,'EdgeColor', 'none');

                        plot3(BASEk(1,:),BASEk(2,:),BASEk(3,:),'b-','linewidth',1');
                        plot3(IPk(1,:),IPk(2,:),IPk(3,:),'r-','linewidth',1);
                        text(POSk(1,5),POSk(2,5),POSk(3,5),num2str(kk),'fontsize',10,'color','k','FontWeight','bold');

                        end;

                        hhh= mesh(YYx,YYy,YYz);
                        set(hhh,'edgecolor',colors(rem(kk-1,6)+1),'linewidth',1); %,'facecolor','none');

                    else

                        figure(4);

                        if show_camera,

                        p1 = struct('vertices',IPk','faces',[1 4 2;2 4 7;2 7 10;2 10 1]);
                        h1 = patch(p1);
                        set(h1,'facecolor',[52 217 160]/255,'EdgeColor', 'r');
                        p2 = struct('vertices',IPk','faces',[1 10 7;7 4 1]);
                        h2 = patch(p2);
                        %set(h2,'facecolor',[236 171 76]/255,'EdgeColor', 'none');
                        set(h2,'facecolor',[247 239 7]/255,'EdgeColor', 'none');

                        plot3(BASEk(1,:),BASEk(2,:),BASEk(3,:),'b-','linewidth',1');
                        plot3(IPk(1,:),IPk(2,:),IPk(3,:),'r-','linewidth',1);
                        hww = text(POSk(1,5),POSk(2,5),POSk(3,5),num2str(kk),'fontsize',10,'color','k','FontWeight','bold');

                        end;

                        plot3(YY_kk(1,:),YY_kk(2,:),YY_kk(3,:),['.' colors(rem(kk-1,6)+1)]);

                    end;

                end;

            end;

        end;

    end;

end;

figure(4);rotate3d on;
axis('equal');
title('Extrinsic parameters (world-centered)');
%view(60,30);
xlabel('X_{world}')
ylabel('Y_{world}')
zlabel('Z_{world}')
view(a,b);
axis vis3d;
axis tight;
grid on;
plot3(3*dX*[1 0 0 0 0],3*dX*[0 0 1 0 0],3*dX*[0 0 0 0 1],'r-','linewidth',3);
hold off;
set(4,'color',[1 1 1]);

set(4,'Name','3D','NumberTitle','off');

%hh = axis;
%hh(5) = 0;
%axis(hh);

%fprintf(1,'To generate the complete movie associated to the optimization loop, try: check_convergence;\n');

if exist('h_switch2')==1,
    if ishandle(h_switch2),
        delete(h_switch2);
    end;
end;

if n_ima ~= 0,
    if show_camera,
        h_switch2 = uicontrol('Parent',4,'Units','normalized', 'Callback','show_camera=0;ext_calib2;', 'Position',[1-.30 0.04  .30  .04],'String','Remove camera reference frames','fontsize',8,'fontname','clean','Tag','Pushbutton1');
    else
        h_switch2 = uicontrol('Parent',4,'Units','normalized', 'Callback','show_camera=1;ext_calib2;', 'Position',[1-.30 0.04  .30  .04],'String','Add camera reference frames','fontsize',8,'fontname','clean','Tag','Pushbutton1');
    end;
end;

if exist('h_switch')==1,
    if ishandle(h_switch),
        delete(h_switch);
    end;
end;

h_switch = uicontrol('Parent',4,'Units','normalized', 'Callback','ext_calib', 'Position',[1-.30 0  .30  .04],'String','Switch to camera-centered view','fontsize',8,'fontname','clean','Tag','Pushbutton1');

figure(4);
rotate3d on;

repo: ekfmonoslam


text:
function result = iscolumnvector(obj)
    result = isvector(obj) && size(obj,2) == 1 && size(obj,1) > 1 && ndims(obj) == 2;
end

author: ziyuw

repo: rembo

save_path: github-repos/MATLAB/ziyuw-rembo


text:
function N = unormal( dom )
%UNORMAL unit normal vector to the unit sphere.
%   N = UNORMAL( DOM ) returns a SPHEREFUNV representing the unit normal
%   vector to the unit sphere using the spherical coordinate domain
%   specified, i.e. DOM = [-pi -pi 0 pi], for co-latitude.
%
%   N = UNORMAL returns a SPHEREFUNV representing the unit normal
%   vector to the unit sphere, using the default spherical coordinates
%   domain, which is co-latitude.
%
%   See also SPHEREFUNV/TANGENTIAL, SPHEREFUNV/DOT, SPHEREFUNV/CURL

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

if ~exist( 'dom', 'var' ) || isempty( dom )
    dom = [-pi pi 0 pi];
end

if ( numel( dom ) ~= 4 )
    error('SPHEREFUN:SPHEREFUNV:normal:domain',...
        ['Domain for the unit vector is incorrect.  It should be a '...
        'double vector with 4 components.']);
end

% Normal vector to the unit sphere.
x = spherefun(@(x,y,z) x, dom);
y = spherefun(@(x,y,z) y, dom);
z = spherefun(@(x,y,z) z, dom);
N = spherefunv(x,y,z);

end