text:
function [outexpdir] = ConvertMWT2Trx(expdir,rootoutputdir,varargin)

%% set parameters

% default parameters
blobsext = 'blobs';
datext = 'dat';
spineext = 'spine';
trxfilestr = 'trx.mat';
perframedirstr = 'perframe';
MINDT = .001;
MAXTIMESTAMPERR = .01;
DEBUG = false;
makesoftlink = true;

perframedict = {'x','x_mm'
  'y','y_mm'
  'area','area_mm'};
trxfns = {'x','y','a','b','theta','x_mm','y_mm','a_mm','b_mm','theta_mm'};

[~,experiment_name] = myfileparts(expdir);

[blobsext,datext,spineext,trxfilestr,perframedirstr,...
  DEBUG,makesoftlink,perframedict,trxfns] = ...
  myparse(varargin,...
  'blobsext',blobsext,...
  'datext',datext,...
  'spineext',spineext,...
  'trxfilestr',trxfilestr,...
  'perframedirstr',perframedirstr,...
  'debug',DEBUG,...
  'makesoftlink',makesoftlink,...
  'perframedict',perframedict,...
  'trxfns',trxfns);

%% read in blobs files

% find blobs files
files = dir(fullfile(expdir,[experiment_name,'*.',blobsext]));
blobsnames = cellfun(@(x) fullfile(expdir,x),{files.name},'UniformOutput',false);

if isempty(blobsnames),
  error('Could not find any blobs files in %s',expdir);
end
fprintf('Reading %d blobs files...\n',numel(blobsnames));
trx = ReadMWTBlobFile(blobsnames);
nflies = numel(trx);

% make sure all the timestamps agree with each other
timestamps = unique([trx.timestamp]);
dts = diff(timestamps);
% sanity check to look for rounding errors
if any(dts < MINDT),
  error('Rounding error check failed: there are timestamps < %f apart',MINDT);
end

% add dt, fps
fps = nanmean(dts);
for i = 1:nflies,
  trx(i).dt = diff(trx(i).timestamp);
  trx(i).fps = fps;
end

%% read in spine files

haveremoved = false;
removeid = false(1,numel(trx));
trxids = [trx.id];

% find dat files
files = dir(fullfile(expdir,[experiment_name,'*.',spineext]));
if isempty(files),
  warning('Spine file does not exist in directory %s',expdir);
else
  spinename = fullfile(expdir,files(1).name);

  fprintf('Reading spine file %s...\n',spinename);

  % read in data
  datacurr = ReadChoreSpineFile(spinename);

  % match up with identities
  dataids = [datacurr.id];
  newremoveid = ~ismember(trxids,dataids);
  if haveremoved && any(newremoveid ~= removeid),
    warning('removeid and newremoveid do not match');
  end
  removeid = removeid | newremoveid;
  haveremoved = true;
  [oldid,idx] = ismember(dataids,trxids);
  if any(~oldid),
    warning('ids found in spine file %s that are not in trx',spinename);
    idx(~oldid) = [];
  end
  for jj = 1:numel(idx),
    j = idx(jj);
    [~,f0] = min(abs(trx(j).timestamp-datacurr(jj).timestamp(1)));
    [~,f1] = min(abs(trx(j).timestamp-datacurr(jj).timestamp(end)));
    if f1-f0+1 ~= size(datacurr(jj).xspine,2),
      error('timestamps don''t match up for spine');
    end
    maxerr = max(abs(trx(j).timestamp(f0:f1)-datacurr(jj).timestamp));
    if maxerr > MAXTIMESTAMPERR,
      error('timestamps don''t match up for spine');
    end
    trx(j).xspine_mm = nan(size(datacurr(jj).xspine,1),trx(j).nframes);
    trx(j).yspine_mm = nan(size(datacurr(jj).xspine,1),trx(j).nframes);
    trx(j).xspine_mm(:,f0:f1) = datacurr(jj).xspine;
    trx(j).yspine_mm(:,f0:f1) = datacurr(jj).yspine;
  end

end

%% read in dat files

% find dat files
files = dir(fullfile(expdir,[experiment_name,'*.',datext]));
datnames = {};
perframefns = cell(0,2);

for i = 1:numel(files),
  match = regexp(files(i).name,['\.(.+)\.',datext,'$'],'tokens','once');
  if isempty(match),
    continue;
  end
  datnames{end+1} = fullfile(expdir,files(i).name); %#ok<AGROW>
  fn = match{1};
  j = find(strcmp(fn,perframedict(:,1)),1);
  if ~isempty(j),
    fn = perframedict{j,2};
  end
  if ismember(fn,trxfns),
    perframefn = fn;
  else
    perframefn = ['perframe_',fn];
  end
  perframefns(end+1,:) = {fn,perframefn}; %#ok<AGROW>

  fprintf('Reading in dat file %s...\n',datnames{end});

  % read in data
  datacurr = ReadChoreDatFile(datnames{end});

  % match up with identities
  dataids = [datacurr.id];
  newremoveid = ~ismember(trxids,dataids);
  if haveremoved && any(newremoveid ~= removeid),
    warning('removeid and newremoveid do not match');
  end
  removeid = removeid | newremoveid;
  haveremoved = true;
  [oldid,idx] = ismember(dataids,trxids);
  if any(~oldid),
    warning('ids found in dat file %s that are not in trx',datnames{end});
    idx(~oldid) = [];
  end
  for jj = 1:numel(idx),
    j = idx(jj);
    [~,f0] = min(abs(trx(j).timestamp-datacurr(jj).timestamp(1)));
    [~,f1] = min(abs(trx(j).timestamp-datacurr(jj).timestamp(end)));
    if f1-f0+1 ~= numel(datacurr(jj).value),
      error('timestamps don''t match up for %s',fn);
    end
    maxerr = max(abs(trx(j).timestamp(f0:f1)-datacurr(jj).timestamp));
    if maxerr > MAXTIMESTAMPERR,
      error('timestamps don''t match up for %s',fn);
    end
    trx(j).(perframefn) = nan(1,trx(j).nframes);
    trx(j).(perframefn)(f0:f1) = datacurr(jj).value;
  end

end

if any(removeid),
  trx(removeid) = [];
end

% replace first frames with nans
for i = 1:size(perframefns,1),
  fn = perframefns{i,2};
  off = inf;
  for j = 1:numel(trx),
    if isempty(trx(j).(fn)),
      continue;
    end
    if trx(j).(fn)(1) ~= 0,
      off = 0;
      break;
    end
    if all(trx(j).(fn) == 0),
      continue;
    end
    k = find(trx(j).(fn) ~= 0,1);
    off = min(off,k-1);
  end
  if off > 0,
    fprintf('First %d frames of data for %s are always 0, setting to nan\n',off,fn);
  end
  for j = 1:numel(trx),
    trx(j).(fn)(1:min(numel(trx(j).(fn)),off)) = nan;
  end
end

% replace last frames with nans
for i = 1:size(perframefns,1),
  fn = perframefns{i,2};
  off = inf;
  for j = 1:numel(trx),
    if isempty(trx(j).(fn)),
      continue;
    end
    if trx(j).(fn)(end) ~= 0,
      off = 0;
      break;
    end
    if all(trx(j).(fn) == 0),
      continue;
    end
    k = find(trx(j).(fn) ~= 0,1,'last');
    off = min(off,numel(trx(j).(fn))-k);
  end
  if off > 0,
    fprintf('Last %d frames of data for %s are always 0, setting to nan\n',off,fn);
  end
  for j = 1:numel(trx),
    trx(j).(fn)(max(1,numel(trx(j).(fn))-off+1):end) = nan;
  end
end

%% compute pxpermm

pxpermm = nanmean([trx.x trx.y] ./ [trx.x_mm trx.y_mm]);
for i = 1:numel(trx),
  trx(i).pxpermm = pxpermm;
  if ~isfield(trx,'x_mm'),
    trx(i).x_mm = trx(i).x / pxpermm;
  end
  if ~isfield(trx,'y_mm'),
    trx(i).y_mm = trx(i).y / pxpermm;
  end
  trx(i).a_mm = trx(i).a / pxpermm;
  trx(i).b_mm = trx(i).b / pxpermm;
  trx(i).theta_mm = trx(i).theta;
  trx(i).xspine = trx(i).xspine_mm * pxpermm;
  trx(i).yspine = trx(i).yspine_mm * pxpermm;
end

%% plot the trajectories

if DEBUG,

ncontours = numel(trx);
nframes = numel(timestamps);
colors = jet(ncontours)*.7;
for t = 1:nframes,
  hold off;
  for i = 1:ncontours,
    if t > trx(i).endframe || t < trx(i).firstframe,
      continue;
    end
    j = t + trx(i).off;
    plot(trx(i).xcontour{j},trx(i).ycontour{j},'.-','color',colors(i,:));
    hold on;
    drawellipse(trx(i).x(j),trx(i).y(j),trx(i).theta(j),trx(i).a(j)*2,trx(i).b(j)*2,'color','r');
  end
  axis equal;
  drawnow;
end

end

%% create output directory

fprintf('Saving...\n');

% make sure the root output directory exists
if ~exist(rootoutputdir,'dir'),
  mkdir(rootoutputdir);
end

% create the experiment directory
outexpdir = fullfile(rootoutputdir,experiment_name);
if ~exist(outexpdir,'dir'),
  mkdir(outexpdir);
end

%% create per-frame directory

perframedir = fullfile(outexpdir,perframedirstr);
if ~exist(perframedir,'dir'),
  mkdir(perframedir);
end

%% save the per-frame data

for i = 1:size(perframefns,1),
  fn = perframefns{i,1};
  perframefn = perframefns{i,2};
  if ~ismember(fn,trxfns),
    data = {trx.(perframefn)}; %#ok<NASGU>
    % units not set yet
    units = parseunits('unit'); %#ok<NASGU>
    outfile = fullfile(perframedir,[fn,'.mat']);
    save(outfile,'data','units');
    trx = rmfield(trx,perframefn);
  end
end

for i = 1:numel(trxfns),
  fn = trxfns{i};
  if isfield(trx,fn),
    data = {trx.(fn)}; %#ok<NASGU>
    % units not set yet
    units = parseunits('unit'); %#ok<NASGU>
    outfile = fullfile(perframedir,[fn,'.mat']);
    save(outfile,'data','units');
  end
end
%% create the trx file

outmatname = fullfile(outexpdir,trxfilestr);
save(outmatname,'trx');


text:
function  [img,coords] = makeMontage(imgCube, sliceList, fileName, numAcross, backVal)
%img = makeMontage(imgCube, [sliceList], [fileName], [numAcross], [backVal])
%
% Compiles a montage image from the images in imgCube.
% (imgCube is x by y by sliceNum)
% sliceList, if specified, determines which slices to extract.
% (it defaults to all slices if it's empty or omitted)
% The image can be displayed with something like:
%
%  figure; imagesc(img); colormap(gray); axis equal; axis off;
%
% If a fileName is specified, the image will be saved to that file.
% Note that the image file format type depends on the filename extension
% (see imwrite).
%
% see also makeMontageIfiles
if(~exist('sliceList', 'var')) sliceList = []; end
if(~exist('numAcross', 'var')) numAcross = []; end
if(~exist('fileName', 'var')) fileName = []; end
if(~exist('backVal', 'var')) backVal = []; end
if(size(imgCube,4)==3)
   [img(:,:,1),coords] = makeMontage(imgCube(:,:,:,1), sliceList, fileName, numAcross,backVal);
   img(:,:,2) = makeMontage(imgCube(:,:,:,2), sliceList, fileName, numAcross,backVal);
   img(:,:,3) = makeMontage(imgCube(:,:,:,3), sliceList, fileName, numAcross,backVal);
   return;
end

% Sanity check
if(any(size(imgCube)>10000))
    error('At least one dimension of input image is >10,000- refusing to continue...');
end

header = 0;
if(isempty(sliceList))
   sliceList = [1:size(imgCube, 3)];
end
sz = size(imgCube);
r = sz(1); c = sz(2);

nImages = length(sliceList);
if(isempty(numAcross))
    numAcross = ceil(sqrt(nImages)*sqrt((r/c)));
end
numDown = ceil(nImages/numAcross);
if(isempty(backVal)) backVal = 0; end
img = ones(r*numDown,c*numAcross)*backVal;
eval(['img = ',class(imgCube),'(img);']);
count = 0;
for ii = 1:length(sliceList)
    curSlice = sliceList(ii);
    count = count+1;
    x = rem(count-1, numAcross)*c;
    y = floor((count-1)/ numAcross)*r;
    img(y+1:y+r,x+1:x+c) = imgCube(:,:,curSlice);
    coords(ii,:) = [x+1,y+1];
end

if(~isempty(fileName))
    imgMin = min(img(:));
    imgMax = max(img(:));
    img2 = uint8(floor((double(img)+imgMin)./(imgMax-imgMin)*255));

    %figure; colormap(repmat([0:1/255:1]',1,3)); image(img2); axis image; axis off;
    [p,f,e] = fileparts(fileName);
    if(isempty(e))
        fileName = [fileName, '.jpg'];
    end
    imwrite(img2, fileName);
    disp(['Wrote montage to ' fullfile(pwd, fileName)]);
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%
% LibQPEP: A Library for Globally Optimal Solving Quadratic Pose Estimation Problems (QPEPs),
%          It also gives highly accurate uncertainty description of the solutions.
%
%
% Article:
%      Wu, J., Zheng, Y., Gao, Z., Jiang, Y., Hu, X., Zhu, Y., Jiao, J., Liu, M. (2020)
%           Quadratic Pose Estimation Problems: Unified Solutions,
%           Solvability/Observability Analysis and Uncertainty Description
%           in A Globally Optimal Framework.
%
%
% Authors:      Jin Wu and Ming Liu
% Affiliation:  Hong Kong University of Science and Technology (HKUST)
% Emails:       jin_wu_uestc@hotmail.com; eelium@ust.hk
% Websites:     https://zarathustr.github.io
%               https://ram-lab.com

function x = sym2vec(X)
n = size(X, 1);
x = [];
for i = 1 : n
    x = [x; X(i, i : n).'];
end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [c4,m4,c2] = cum4(X,prewhiten)
%CUM4 Fourth-order cumulant tensor.
%   [c4,m4,c2] = cum4(X) computes the second-order cumulant (covariance
%   matrix) c2, fourth-order moment m4 and fourth-order cumulant
%   (quadricovariance tensor) c4 of a matrix X in which each row is an
%   observation and each column is a variable. Herein,
%
%      c2(i,j)     = E[xi.*conj(xj)]
%      m4(i,j,k,l) = E[xi.*conj(xj).*conj(xk).*xl]
%      c4(i,j,k,l) = E[xi.*conj(xj).*conj(xk).*xl] ...
%                    - E[xi.*conj(xj)]*E[conj(xk).*xl] ...
%                    - E[xi.*conj(xk)]*E[conj(xj).*xl] ...
%                    - E[xi.*xl]*E[conj(xj).*conj(xk)]
%
%   where the expectation E is approximated by the arithmetic mean and xi
%   is the i-th mean centered variable, X(:,i)-mean(X(:,i)) (and
%   analogously for xj, xk and xl).
%
%   [c4,m4,c2] = cum4(X,'prewhiten') applies a linear transformation to the
%   columns of X so that the covariance matrix of the new matrix is the
%   identity matrix before computing its fourth-order cumulant.
%
%   See also cov, scov.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] P. McCullagh, "Tensor Methods in Statistics," Chapman and Hall,
%       London, 1987.
%   [2] C. Nikias, A. Petropulu, "Higher-Order Spectra Analysis: A
%       Nonlinear Signal Processing Framework," Prentice Hall, 1993.

% Check the prewhiten option.
if nargin < 3, prewhiten = false; end
if ischar(prewhiten), prewhiten = strcmpi(prewhiten,'prewhiten'); end

% Center the variables.
X = bsxfun(@minus,X,mean(X));

% Apply a prewhitening to X if requested.
n = size(X,1);
if prewhiten
    [U,S,~] = svd(X,'econ');
    X = U*(S*pinv(S))*sqrt(n);
end

% Compute c2 = E[xi*conj(xj)] and r2 = E[xi*xj].
c2 = conj(X'*X)/n;
r2 = (X.'*X)/n;

% Compute m4 = E[xi*conj(xj)*conj(xk)*xl].
% Introduce singleton dimensions 3 and 4.
m4 = bsxfun(@times,permute(X,[2 3 4 5 1]),permute(conj(X),[3 2 4 5 1]));
m4 = mean(bsxfun(@times,m4,permute(m4,[3 4 2 1 5])),5);

% Compute c4(i,j,k,l) = m4 - E[xi.*conj(xj)]*E[conj(xk).*xl] - ...
% E[xi.*conj(xk)]*E[conj(xj).*xl] - E[xi.*xl]*E[conj(xj).*conj(xk)].
c4 = m4-bsxfun(@times,c2,permute(c2,[3 4 2 1])) ...
       -bsxfun(@times,permute(c2,[1 3 2 4]),permute(c2,[3 2 4 1])) ...
       -bsxfun(@times,permute(r2,[1 3 4 2]),permute(conj(r2),[3 1 2 4]));

author: andrewssobral

repo: lrslibrary


text:
function [blocks,idx] = my_im2col(I,blkSize,slidingDis);
if (slidingDis==1)
    blocks = im2col(I,blkSize,'sliding');
    idx = [1:size(blocks,2)];
    return
end

idxMat = zeros(size(I)-blkSize+1);
idxMat([[1:slidingDis:end-1],end],[[1:slidingDis:end-1],end]) = 1; % take blocks in distances of 'slidingDix', but always take the first and last one (in each row and column).
idx = find(idxMat);
[rows,cols] = ind2sub(size(idxMat),idx);
blocks = zeros(prod(blkSize),length(idx));
for i = 1:length(idx)
    currBlock = I(rows(i):rows(i)+blkSize(1)-1,cols(i):cols(i)+blkSize(2)-1);
    blocks(:,i) = currBlock(:);
end

author: lbasek

repo: image-denoising-benchmark


text:
%%*****************************************************************************
%% sqlp: solve an semidefinite-quadratic-linear program
%%       by infeasible path-following method.
%%
%%  [obj,X,y,Z,info,runhist] = sqlp(blk,At,C,b,OPTIONS,X0,y0,Z0);
%%
%%  Input: blk: a cell array describing the block diagonal structure of SQL data.
%%          At: a cell array with At{p} = [svec(Ap1) ... svec(Apm)]
%%         b,C: data for the SQL instance.
%%  (X0,y0,Z0): an initial iterate (if it is not given, the default is used).
%%     OPTIONS: a structure that specifies parameters required in sqlp.m,
%%              (if it is not given, the default in sqlparameters.m is used).
%%
%%  Output: obj  = [<C,X> <b,y>].
%%          (X,y,Z): an approximately optimal solution or a primal or dual
%%                   infeasibility certificate.
%%          info.termcode = termination-code
%%          info.iter     = number of iterations
%%          info.obj      = [primal-obj, dual-obj]
%%          info.cputime  = total-time
%%          info.gap      = gap
%%          info.pinfeas  = primal_infeas
%%          info.dinfeas  = dual_infeas
%%          runhist.pobj    = history of primal objective value.
%%          runhist.dobj    = history of dual   objective value.
%%          runhist.gap     = history of <X,Z>.
%%          runhist.pinfeas = history of primal infeasibility.
%%          runhist.dinfeas = history of dual   infeasibility.
%%          runhist.cputime = history of cputime spent.
%%----------------------------------------------------------------------------
%%  The OPTIONS structure specifies the required parameters:
%%      vers  gam  predcorr  expon  gaptol  inftol  steptol
%%      maxit  printlevel  scale_data ...
%%      (all have default values set in sqlparameters.m).
%%
%%*************************************************************************
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*************************************************************************

  function [obj,X,y,Z,info,runhist] = sqlp(blk,At,C,b,OPTIONS,X0,y0,Z0);
%%
%%-----------------------------------------
%% get parameters from the OPTIONS structure.
%%-----------------------------------------
%%
   global matlabversion ispc_hp_ibm
   global spdensity  iter  solve_ok  switch2LU  depconstr
   global cachesize  smallblkdim  printlevel
   global schurfun   schurfun_par

   warning off;
   matlabversion = sscanf(version,'%f');
   matlabversion = matlabversion(1);
   ispc_hp_ibm = strncmp(computer,'PC',2) | strncmp(computer,'HP',2) | ...
                 strncmp(computer,'IBM',3);

   vers        = 1;
   predcorr    = 1;
   gam         = 0;
   expon       = 1;
   gaptol      = 1e-8;
   inftol      = 1e-8;
   steptol     = 1e-6;
   maxit       = 100;
   printlevel  = 3;
   stoplevel   = 1;
   scale_data  = 0;
   spdensity   = 0.5;
   rmdepconstr = 0;
   cachesize   = 256;
   smallblkdim = 15;
   schurfun     = cell(size(blk,1),1);
   schurfun_par = cell(size(blk,1),1);
   if exist('OPTIONS')
      if isfield(OPTIONS,'vers');        vers     = OPTIONS.vers; end
      if isfield(OPTIONS,'predcorr');    predcorr = OPTIONS.predcorr; end
      if isfield(OPTIONS,'gam');         gam      = OPTIONS.gam; end
      if isfield(OPTIONS,'expon');       expon    = OPTIONS.expon; end
      if isfield(OPTIONS,'gaptol');      gaptol   = OPTIONS.gaptol; end
      if isfield(OPTIONS,'inftol');      inftol   = OPTIONS.inftol; end
      if isfield(OPTIONS,'steptol');     steptol  = OPTIONS.steptol; end
      if isfield(OPTIONS,'maxit');       maxit    = OPTIONS.maxit; end
      if isfield(OPTIONS,'printlevel');  printlevel  = OPTIONS.printlevel; end
      if isfield(OPTIONS,'stoplevel');   stoplevel   = OPTIONS.stoplevel; end
      if isfield(OPTIONS,'scale_data');  scale_data  = OPTIONS.scale_data; end
      if isfield(OPTIONS,'spdensity');   spdensity   = OPTIONS.spdensity; end
      if isfield(OPTIONS,'rmdepconstr'); rmdepconstr = OPTIONS.rmdepconstr; end
      if isfield(OPTIONS,'cachesize');   cachesize   = OPTIONS.cachesize; end
      if isfield(OPTIONS,'smallblkdim'); smallblkdim = OPTIONS.smallblkdim; end
      if isfield(OPTIONS,'schurfun');
         schurfun = OPTIONS.schurfun;
         if ~isempty(schurfun); scale_data = 0; end
      end
      if isfield(OPTIONS,'schurfun_par'); schurfun_par = OPTIONS.schurfun_par; end
      if isempty(schurfun); schurfun = cell(size(blk,1),1); end
      if isempty(schurfun_par); schurfun_par = cell(size(blk,1),1); end
   end
%%
   if all(vers-[1 2]); error('*** vers must be 1 or 2 ***'); end;
%%
%%-----------------------------------------
%% convert matrices to cell arrays.
%%-----------------------------------------
%%
   if ~iscell(At); At = {At}; end;
   if ~iscell(C);  C = {C}; end;
   m = length(b);
   if all(size(At) == [size(blk,1), m]);
      convertyes = zeros(size(blk,1),1);
      for p = 1:size(blk,1)
         if strcmp(blk{p,1},'s') & all(size(At{p,1}) == sum(blk{p,2}))
            convertyes(p) = 1;
         end
      end
      if any(convertyes)
         if (printlevel); fprintf('\n sqlp: converting At into required format'); end
         At = svec(blk,At,ones(size(blk,1),1));
      end
   end
   if (nargin <= 5) | (isempty(X0) | isempty(y0) | isempty(Z0));
      [X0,y0,Z0] = infeaspt(blk,At,C,b);
   end
   X = X0; y = y0; Z = Z0;
   if ~iscell(X);  X = {X}; end;
   if ~iscell(Z);  Z = {Z}; end;
%%
%%-----------------------------------------
%% validate SQLP data.
%%-----------------------------------------
%%
   tstart = cputime;
   [blk,At,C,b,dim,numblk,X,Z] = validate(blk,At,C,b,X,y,Z);
   if (printlevel>=2)
      fprintf('\n num. of constraints = %2.0d',length(b));
      if dim(1);
         fprintf('\n dim. of sdp    var  = %2.0d,',dim(1));
         fprintf('   num. of sdp  blk  = %2.0d',numblk(1));
      end
      if dim(2);
         fprintf('\n dim. of socp   var  = %2.0d,',dim(2));
         fprintf('   num. of socp blk  = %2.0d',numblk(2));
      end
      if dim(3); fprintf('\n dim. of linear var  = %2.0d',dim(3)); end
      if dim(4); fprintf('\n dim. of free   var  = %2.0d',dim(4)); end
   end
%%
%%-----------------------------------------
%% convert unrestricted blk to linear blk.
%%-----------------------------------------
%%
   ublkidx = zeros(size(blk,1),1);
   for p = 1:size(blk,1)
      if strcmp(blk{p,1},'u')
         ublkidx(p) = 1;
         n = 2*blk{p,2};
         blk{p,1} = 'l';
         blk{p,2} = n;
         At{p} = [At{p}; -At{p}];
         C{p} = [C{p}; -C{p}];
         b2 = 1 + abs(b');
         normC = 1+norm(C{p});
         normA = 1+sqrt(sum(At{p}.*At{p}));
         X{p} = max(1,max(b2./normA)) *ones(n,1);
         Z{p} = max(1,max([normA,normC])/sqrt(n)) *ones(n,1);
      end
   end
%%
%%-----------------------------------------
%% check whether {A1,...,Am} is
%% linearly independent.
%%-----------------------------------------
%%
   m0 = length(b);
   [At,b,y,indeprows,depconstr,feasible] = checkdepconstr(blk,At,b,y,rmdepconstr);
   if (~feasible)
      fprintf('\n sqlp: SQLP is not feasible'); return;
   end
%%
%%-----------------------------------------
%% scale SQLP data. Note: must be done only
%% after checkdepconstr
%%-----------------------------------------
%%
   normC2 = zeros(length(C),1);
   for p = 1:length(C); normC2(p) = max(max(abs(C{p}))); end
   normC2 = 1+max(normC2);
   normb2 = 1+max(abs(b));
   normX0 = 1+ops(X0,'norm'); normZ0 = 1+ops(Z0,'norm');
   if (scale_data)
      [At,C,b,X,y,Z,normA,normC,normb] = scaling(blk,At,C,b,X,y,Z);
   else
      normA = 1; normC = 1; normb = 1;
   end
%%
%%-----------------------------------------
%% find the combined list of non-zero
%% elements of Aj, j = 1:k, for each k.
%%-----------------------------------------
%%
   m = length(b);
   [At,C,X,Z,par.permA,par.permZ] = sortA(blk,At,C,b,X,Z);
   [par.isspA,par.nzlistA,par.nzlistAsum,par.isspAy,par.nzlistAy] = nzlist(blk,At,m);
%%
%%-----------------------------------------
%% initialization
%%-----------------------------------------
%%
   [Xchol,indef(1)] = blkcholfun(blk,X);
   [Zchol,indef(2)] = blkcholfun(blk,Z);
   if any(indef)
      if (printlevel); fprintf('\n Stop: X, Z are not both positive definite'); end
      termcode = -3;
      return;
   end
   nn = 0;
   for p = 1:size(blk,1);
      pblk = blk(p,:);
      if strcmp(pblk{1},'s') | strcmp(pblk{1},'q') | strcmp(pblk{1},'l')
         nn = nn + sum(pblk{2});
      end
   end
   AX = AXfun(blk,At,par.permA,X);
   rp = b-AX;
   ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,y));
   ZpATynorm = ops(ZpATy,'norm');
   Rd = ops(C,'-',ZpATy);
   obj = (normb*normC)*[blktrace(blk,C,X), b'*y];
   trXZ = blktrace(blk,X,Z);
   gap = (normb*normC)*trXZ;
   mu  = trXZ/nn;
   rel_gap = gap/(1+sum(abs(obj)));
   prim_infeas = norm(rp)/normb2;
   dual_infeas = ops(Rd,'norm')/normC2;
   infeas_meas = max(prim_infeas,dual_infeas);
   if (scale_data)
      infeas_org(1) = norm(normA.*rp)*normb/normb2;
      infeas_org(2) = ops(Rd,'norm')*normC/normC2;
   else
      infeas_org = [0,0];
   end
   termcode = -6;
   pstep = 0; dstep = 0; pred_convg_rate = 1; corr_convg_rate = 1;
   prim_infeas_bad = 0;  homRd = inf; homrp = inf;
   runhist.pobj = obj(1);
   runhist.dobj = obj(2);
   runhist.gap  = gap;
   runhist.relgap  = rel_gap;
   runhist.pinfeas = prim_infeas;
   runhist.dinfeas = dual_infeas;
   runhist.infeas  = infeas_meas;
   runhist.step    = 0;
   runhist.cputime = cputime-tstart;
   ttime.preproc   = runhist.cputime;
   ttime.pred = 0; ttime.pred_pstep = 0; ttime.pred_dstep = 0;
   ttime.corr = 0; ttime.corr_pstep = 0; ttime.corr_dstep = 0;
   ttime.pchol = 0; ttime.dchol = 0; ttime.misc = 0;
%%
%%-----------------------------------------
%% display parameters and initial info
%%-----------------------------------------
%%
   if (printlevel >= 2)
      fprintf('\n********************************************');
      fprintf('***********************\n');
      fprintf('   SDPT3: Infeasible path-following algorithms');
      fprintf('\n********************************************');
      fprintf('***********************\n');
      [hh,mm,ss] = mytimed(ttime.preproc);
      if (printlevel>=3)
         fprintf(' version  predcorr  gam  expon  scale_data\n');
         if (vers == 1); fprintf('   HKM '); elseif (vers == 2); fprintf('    NT '); end
         fprintf('     %1.0f      %4.3f',predcorr,gam);
         fprintf('   %1.0f        %1.0f    %1.0f\n',expon,scale_data);
         fprintf('\nit  pstep dstep p_infeas d_infeas  gap')
         fprintf('     mean(obj)    cputime\n');
         fprintf('------------------------------------------------');
         fprintf('-------------------\n');
         fprintf('%2.0f  %4.3f %4.3f %2.1e %2.1e',0,0,0,prim_infeas,dual_infeas);
         fprintf('  %2.1e %- 7.6e  %d:%d:%d',gap,mean(obj),hh,mm,ss);
      end
   end
%%
%%---------------------------------------------------------------
%% start main loop
%%---------------------------------------------------------------
%%
   param.inftol      = inftol;
   param.normA       = normA;
   param.normC       = normC;
   param.normb       = normb;
   param.normX0      = normX0;
   param.normZ0      = normZ0;
   param.m0          = m0;
   param.indeprows   = indeprows;
   param.gaptol      = gaptol;
   param.inftol      = inftol;
   param.scale_data  = scale_data;
   param.printlevel  = printlevel;
%%
   for iter = 1:maxit;

       update_iter = 0; breakyes = 0; pred_slow = 0; corr_slow = 0; step_short = 0;
       tstart  = cputime;
       timeold = cputime;
%%
%%---------------------------------------------------------------
%% predictor step.
%%---------------------------------------------------------------
%%
       if (predcorr)
          sigma = 0;
       else
          sigma = 1-0.9*min(pstep,dstep);
          if (iter == 1); sigma = 0.5; end;
       end
       sigmu = sigma*mu;

       invXchol = cell(size(blk,1),1);
       invZchol = ops(Zchol,'inv');
       if (vers == 1);
          [par,dX,dy,dZ,coeff,L,hRd] = ...
           HKMpred(blk,At,par,rp,Rd,sigmu,X,Z,invZchol);
       elseif (vers == 2);
          [par,dX,dy,dZ,coeff,L,hRd] = ...
           NTpred(blk,At,par,rp,Rd,sigmu,X,Z,Zchol,invZchol);
       end
       if (solve_ok <= 0)
          fprintf('\n Stop: difficulty in computing predictor directions');
          runhist.cputime(iter+1) = cputime-tstart;
          termcode = -4;
          break;
       end
       timenew = cputime;
       ttime.pred = ttime.pred + timenew-timeold; timeold = timenew;
%%
%%-----------------------------------------
%% step-lengths for predictor step
%%-----------------------------------------
%%
      if (gam == 0)
         gamused = 0.9 + 0.09*min(pstep,dstep);
      else
         gamused = gam;
      end
      [Xstep,invXchol] = steplength(blk,X,dX,Xchol,invXchol);
      pstep = min(1,gamused*Xstep);
      if (Xstep > .99e12) & (blktrace(blk,C,dX) < -1e-3) & (prim_infeas < 1e-3)
         pstep = Xstep;
         if (printlevel); fprintf('\n Predictor: dual seems infeasible.'); end
      end
      timenew = cputime;
      ttime.pred_pstep = ttime.pred_pstep + timenew-timeold; timeold = timenew;
      Zstep = steplength(blk,Z,dZ,Zchol,invZchol);
      dstep = min(1,gamused*Zstep);
      if (Zstep > .99e12) & (b'*dy > 1e-3) & (dual_infeas < 1e-3)
         dstep = Zstep;
         if (printlevel); fprintf('\n Predictor: primal seems infeasible.'); end
      end
      trXZpred = trXZ + pstep*blktrace(blk,dX,Z) + dstep*blktrace(blk,X,dZ) ...
                 + pstep*dstep*blktrace(blk,dX,dZ);
      gappred = (normb*normC)*trXZpred;
      mupred  = trXZpred/nn;
      mupredhist(iter) = mupred;
      timenew = cputime;
      ttime.pred_dstep = ttime.pred_dstep + timenew-timeold; timeold = timenew;
%%
%%-----------------------------------------
%%  stopping criteria for predictor step.
%%-----------------------------------------
%%
      if (min(pstep,dstep) < steptol) & (stoplevel)
         if (printlevel)
            fprintf('\n Stop: steps in predictor too short:');
            fprintf(' pstep = %3.2e,  dstep = %3.2e\n',pstep,dstep);
         end
         runhist.cputime(iter+1) = cputime-tstart;
         termcode = -2;
         breakyes = 1;
      end
      if (iter >= 2)
         idx = [max(2,iter-2) : iter];
         pred_slow = all(mupredhist(idx)./mupredhist(idx-1) > 0.4);
         idx = [max(2,iter-5) : iter];
         pred_convg_rate = mean(mupredhist(idx)./mupredhist(idx-1));
         pred_slow = pred_slow + (mupred/mu > 5*pred_convg_rate);
      end
      if (~predcorr)
         if (max(mu,infeas_meas) < 1e-6) & (pred_slow) & (stoplevel)
            if (printlevel)
               fprintf('\n  Stop: lack of progress in predictor:');
               fprintf(' mupred/mu = %3.2f, pred_convg_rate = %3.2f.',...
                         mupred/mu,pred_convg_rate);
            end
            runhist.cputime(iter+1) = cputime-tstart;
            termcode = -1;
            breakyes = 1;
         else
            update_iter = 1;
         end
      end
%%
%%---------------------------------------------------------------
%% corrector step.
%%---------------------------------------------------------------
%%
      if (predcorr) & (~breakyes)
         step_pred = min(pstep,dstep);
         if (mu > 1e-6)
            if (step_pred < 1/sqrt(3));
               expon_used = 1;
            else
               expon_used = max(expon,3*step_pred^2);
            end
         else
            expon_used = max(1,min(expon,3*step_pred^2));
         end
         sigma = min( 1, (mupred/mu)^expon_used );
         sigmu = sigma*mu;
%%
         if (vers == 1)
            [dX,dy,dZ] = HKMcorr(blk,At,par,rp,Rd,sigmu,hRd,...
             dX,dZ,coeff,L,X,Z);
         elseif (vers == 2)
            [dX,dy,dZ] = NTcorr(blk,At,par,rp,Rd,sigmu,hRd,...
             dX,dZ,coeff,L,X,Z);
         end
         if (solve_ok <= 0)
            fprintf('\n Stop: difficulty in computing corrector directions');
            runhist.cputime(iter+1) = cputime-tstart;
            termcode = -4;
            break;
         end
         timenew = cputime;
         ttime.corr = ttime.corr + timenew-timeold; timeold = timenew;
%%
%%-----------------------------------
%% step-lengths for corrector step
%%-----------------------------------
%%
         if (gam == 0)
            gamused = 0.9 + 0.09*min(pstep,dstep);
         else
            gamused = gam;
         end
         Xstep = steplength(blk,X,dX,Xchol,invXchol);
         pstep = min(1,gamused*Xstep);
         if (Xstep > .99e12) & (blktrace(blk,C,dX) < -1e-3) & (prim_infeas < 1e-3)
            pstep = Xstep;
            if (printlevel); fprintf('\n Corrector: dual seems infeasible.'); end
         end
         timenew = cputime;
         ttime.corr_pstep = ttime.corr_pstep + timenew-timeold; timeold = timenew;
         Zstep = steplength(blk,Z,dZ,Zchol,invZchol);
         dstep = min(1,gamused*Zstep);
         if (Zstep > .99e12) & (b'*dy > 1e-3) & (dual_infeas < 1e-3)
            dstep = Zstep;
            if (printlevel); fprintf('\n Corrector: primal seems infeasible.'); end
         end
         trXZcorr = trXZ + pstep*blktrace(blk,dX,Z) + dstep*blktrace(blk,X,dZ)...
                    + pstep*dstep*blktrace(blk,dX,dZ);
         gapcorr = (normb*normC)*trXZcorr;
         mucorr  = trXZcorr/nn;
         timenew = cputime;
         ttime.corr_dstep = ttime.corr_dstep + timenew-timeold; timeold = timenew;
%%
%%-----------------------------------------
%%  stopping criteria for corrector step
%%-----------------------------------------
%%
         if (iter >= 2)
            idx = [max(2,iter-2) : iter];
            corr_slow = all(runhist.gap(idx)./runhist.gap(idx-1) > 0.8);
            idx = [max(2,iter-5) : iter];
            corr_convg_rate = mean(runhist.gap(idx)./runhist.gap(idx-1));
            corr_slow = corr_slow + (mucorr/mu > max(min(1,5*corr_convg_rate),0.8));
         end
	 if (max(mu,infeas_meas) < 1e-6) & (iter > 10) & (corr_slow) & (stoplevel)
   	    if (printlevel)
               fprintf('\n  Stop: lack of progress in corrector:');
               fprintf(' mucorr/mu = %3.2f, corr_convg_rate = %3.2f',...
                         mucorr/mu,corr_convg_rate);
            end
            runhist.cputime(iter+1) = cputime-tstart;
            termcode = -1;
            breakyes = 1;
         else
            update_iter = 1;
         end
      end
%%
%%---------------------------------------------------------------
%% udpate iterate
%%---------------------------------------------------------------
%%
      indef = [1,1];
      if (update_iter)
         for t = 1:5
            [Xchol,indef(1)] = blkcholfun(blk,ops(X,'+',dX,pstep));
            timenew = cputime;
            ttime.pchol = ttime.pchol + timenew-timeold; timeold = timenew;
            if (indef(1)); pstep = 0.8*pstep; else; break; end
         end
	 if (t > 1); pstep = gamused*pstep; end
	 for t = 1:5
            [Zchol,indef(2)] = blkcholfun(blk,ops(Z,'+',dZ,dstep));
            timenew = cputime;
            ttime.dchol = ttime.dchol + timenew-timeold; timeold = timenew;
            if (indef(2)); dstep = 0.8*dstep; else; break; end
         end
	 if (t > 1); dstep = gamused*dstep; end
         AXtmp = AX + pstep*AXfun(blk,At,par.permA,dX);
         prim_infeasnew = norm(b-AXtmp)/normb2;
         if any(indef)
            if (printlevel); fprintf('\n Stop: X, Z not both positive definite'); end
            termcode = -3;
            breakyes = 1;
         elseif (prim_infeasnew > max([rel_gap,20*prim_infeas,1e-8])) ...
	    | (prim_infeasnew > max([1e-4,20*prim_infeas]) & (switch2LU))
            if (stoplevel) & (max(pstep,dstep)<=1)
               if (printlevel)
                  fprintf('\n Stop: primal infeas has deteriorated too much, %2.1e',prim_infeasnew);
               end
               termcode = -7;
               breakyes = 1;
            end
         else
            X = ops(X,'+',dX,pstep);
            y = y+dstep*dy;
            Z = ops(Z,'+',dZ,dstep);
         end
      end
%%---------------------------------------------------------------
%% adjust linear blk arising from unrestricted blk
%%---------------------------------------------------------------
%%
      for p = 1:size(blk,1)
         if (ublkidx(p) == 1)
            len = blk{p,2}/2;
            alpha = 0.8;
            xtmp = min(X{p}([1:len]),X{p}(len+[1:len]));
            X{p}([1:len]) = X{p}([1:len]) - alpha*xtmp;
            X{p}(len+[1:len]) = X{p}(len+[1:len]) - alpha*xtmp;
            if (mu < 1e-8)
               Z{p} = 0.5*mu./max(1,X{p});
	    else
               ztmp = min(1,max(Z{p}([1:len]),Z{p}(len+[1:len])));
               beta1 = xtmp'*(Z{p}([1:len])+Z{p}(len+[1:len]));
               beta2 = (X{p}([1:len])+X{p}(len+[1:len])-2*xtmp)'*ztmp;
               beta = max(0.1,min(beta1/beta2,0.5));
               Z{p}([1:len]) = Z{p}([1:len]) + beta*ztmp;
               Z{p}(len+[1:len]) = Z{p}(len+[1:len]) + beta*ztmp;
            end
         end
      end
%%
%%---------------------------------------------------------------
%% compute rp, Rd, infeasibities, etc.
%%---------------------------------------------------------------
%%
      trXZ = blktrace(blk,X,Z);
      gap = (normb*normC)*trXZ;
      mu  = trXZ/nn;
      AX  = AXfun(blk,At,par.permA,X);
      rp  = b-AX;
      ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,y));
      ZpATynorm = ops(ZpATy,'norm');
      Rd  = ops(C,'-',ZpATy);
      obj = (normb*normC)*[blktrace(blk,C,X), b'*y];
      rel_gap = gap/(1+sum(abs(obj)));
      prim_infeas  = norm(rp)/normb2;
      dual_infeas = ops(Rd,'norm')/normC2;
      if (scale_data)
         infeas_org(1) = norm(normA.*rp)*normb/normb2;
         infeas_org(2) = ops(Rd,'norm')*normC/normC2;
      end
      infeas_meas = max(prim_infeas,dual_infeas);
      if (obj(2) > 0); homRd = ZpATynorm/(obj(2)/(normb*normC)); else; homRd = inf; end
      if (obj(1) < 0); homrp = norm(AX)/(-obj(1)/(normb*normC)); else; homrp = inf; end
      runhist.pobj(iter+1)  = obj(1);
      runhist.dobj(iter+1)  = obj(2);
      runhist.gap(iter+1)   = gap;
      runhist.relgap(iter+1)  = rel_gap;
      runhist.pinfeas(iter+1) = prim_infeas;
      runhist.dinfeas(iter+1) = dual_infeas;
      runhist.infeas(iter+1)  = infeas_meas;
      runhist.step(iter+1)    = min(pstep,dstep);
      runhist.cputime(iter+1) = cputime-tstart;
      timenew = cputime;
      ttime.misc = ttime.misc + timenew-timeold; timeold = timenew;
      [hh,mm,ss] = mytimed(sum(runhist.cputime));
      if (printlevel>=3)
         fprintf('\n%2.0f  %4.3f %4.3f',iter,pstep,dstep);
         fprintf(' %2.1e %2.1e  %2.1e',prim_infeas,dual_infeas,gap);
         fprintf(' %- 7.6e  %d:%d:%d',mean(obj),hh,mm,ss);
      end
%%
%%--------------------------------------------------
%% check convergence.
%%--------------------------------------------------
%%
      param.iter        = iter;
      param.obj         = obj;
      param.rel_gap     = rel_gap;
      param.gap         = gap;
      param.mu          = mu;
      param.prim_infeas = prim_infeas;
      param.dual_infeas = dual_infeas;
      param.homRd       = homRd;
      param.homrp       = homrp;
      param.AX          = AX;
      param.ZpATynorm   = ZpATynorm;
      param.normX       = ops(X,'norm');
      param.normZ       = ops(Z,'norm');
      param.termcode    = termcode;
      param.stoplevel   = stoplevel;
      param.prim_infeas_bad = prim_infeas_bad;
%%
      if (~breakyes)
         [termcode,breakyes,prim_infeas_bad,restart] = sqlpcheckconvg(param,runhist);
      end
      if (breakyes); break; end
      if (restart)
         [X,y,Z] = infeaspt(blk,At,C,b,2,1e5);
         trXZ = blktrace(blk,X,Z);
         gap = (normb*normC)*trXZ;
         mu  = trXZ/nn;
         rp  = b-AXfun(blk,At,par.permA,X);
         ZpATy = ops(Z,'+',Atyfun(blk,At,par.permA,par.isspAy,y));
         Rd  = ops(C,'-',ZpATy);
         prim_infeas  = norm(rp)/normb2;
         dual_infeas = ops(Rd,'norm')/normC2;
         infeas_meas = max(prim_infeas,dual_infeas);
         [Xchol,indef(1)] = blkcholfun(blk,X);
         [Zchol,indef(2)] = blkcholfun(blk,Z);
         stoplevel = 3;
      end
   end
%%
%%---------------------------------------------------------------
%% end of main loop
%%---------------------------------------------------------------
%%
%%---------------------------------------------------------------
%% unscale and produce infeasibility certificates if appropriate
%%---------------------------------------------------------------
%%
   if (iter >= 1)
      param.termcode = termcode;
      [X,y,Z,termcode,resid,reldist] = sqlpmisc(blk,At,C,b,X,y,Z,par.permZ,param);
   end
%%
%%---------------------------------------------------------------
%% recover unrestricted blk from linear blk
%%---------------------------------------------------------------
%%
   for p = 1:size(blk,1)
      if (ublkidx(p) == 1)
         n = blk{p,2}/2;
         X{p} = X{p}(1:n)-X{p}(n+[1:n]);
         Z{p} = Z{p}(1:n);
      end
   end
%%
%%---------------------------------------------------------------
%% print summary
%%---------------------------------------------------------------
%%
   if (scale_data)
      dimacs = [infeas_org(1); 0; infeas_org(2); 0];
   else
      dimacs = [prim_infeas; 0; dual_infeas; 0];
   end
   dimacs = [dimacs; [-diff(obj); gap]/(1+sum(abs(obj)))];
   info.dimacs   = dimacs;
   info.termcode = termcode;
   info.iter     = iter;
   info.obj      = obj;
   info.gap      = gap;
   info.relgap   = rel_gap;
   info.pinfeas  = prim_infeas;
   info.dinfeas  = dual_infeas;
   info.cputime  = sum(runhist.cputime);
   info.resid    = resid;
   info.reldist  = reldist;
%%
   nnorm.b = norm(b); nnorm.C = ops(C,'norm'); nnorm.A = ops(At,'norm');
   nnorm.X = ops(X,'norm'); nnorm.y = norm(y); nnorm.Z = ops(Z,'norm');
   sqlpsummary(info,ttime,infeas_org,nnorm,printlevel);
%%*****************************************************************************

author: yu-jiang

repo: radpbook

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


text:
function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.

% Create New Figure
figure; hold on;

% ====================== YOUR CODE HERE ======================
% Instructions: Plot the positive and negative examples on a
%               2D plot, using the option 'k+' for the positive
%               examples and 'ko' for the negative examples.
%

pos = find(y == 1);
neg = find(y == 0);
plot(X(pos, 1), X(pos, 2), 'k+' ,'LineWidth', 2, 'MarkerSize', 7);
plot(X(neg, 1), X(neg, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);

% =========================================================================

hold off;

end

author: scruel

repo: Notes-ML-AndrewNg


text:
% distance to wall
function [data,units] = compute_distnose2wall_rect(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  xnose = trx(fly).x_mm + 2*trx(fly).a_mm.*cos(trx(fly).theta_mm);
  ynose = trx(fly).y_mm + 2*trx(fly).a_mm.*sin(trx(fly).theta_mm);

  % The ordering of the arguments is to account for weird y directions of images.
  dtop = getDist(trx.landmark_params{n}.tl_x(fly),trx.landmark_params{n}.tl_y(fly), trx.landmark_params{n}.tr_x(fly),trx.landmark_params{n}.tr_y(fly),...
    xnose,ynose);
  dleft = getDist(trx.landmark_params{n}.bl_x(fly),trx.landmark_params{n}.bl_y(fly), trx.landmark_params{n}.tl_x(fly),trx.landmark_params{n}.tl_y(fly),...
    xnose,ynose);
  dright = getDist(trx.landmark_params{n}.tr_x(fly),trx.landmark_params{n}.tr_y(fly), trx.landmark_params{n}.br_x(fly),trx.landmark_params{n}.br_y(fly),...
    xnose,ynose);
  dbottom = getDist(trx.landmark_params{n}.br_x(fly),trx.landmark_params{n}.br_y(fly), trx.landmark_params{n}.bl_x(fly),trx.landmark_params{n}.bl_y(fly),...
    xnose,ynose);

  data{i} = min([dtop;dleft; dright; dbottom],[],1);

end
units = parseunits('mm');

function d = getDist(p1_x,p1_y,p2_x,p2_y,p_x,p_y)

dp1p2 = (p1_x-p2_x).^2 + (p1_y-p2_y).^2;

dotpr_x = (p1_x - p_x).*(p1_x-p2_x);
dotpr_y = (p1_y - p_y).*(p1_y-p2_y);

t = (dotpr_x+dotpr_y)/dp1p2;
proj_x = p1_x + t.*(p2_x-p1_x);
proj_y = p1_y + t.*(p2_y-p1_y);
d = sqrt(  (p_x-proj_x).^2 + (p_y-proj_y).^2);

% Find whether the mice is in the interior of the square.
side = sign((p2_x-p1_x)*(p_y-p1_y)-(p2_y-p1_y)*(p_x-p1_x));
d = side.*d;

author: kristinbranson

repo: JAABA


text:
function lik = lik_inputdependentnoise(varargin)
%lik_inputdependentnoise    Create input-dependent noise likelihood structure
%
%  Description
%    LIK = LIK_INPUTDEPENDENTNOISE('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a Gaussian likelihood with input dependent noise structure
%    in which the named parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    LIK = LIK_INPUTDEPENDENTNOISE(LIK,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a likelihood function structure with the named
%    parameters altered with the specified values.
%
%    Parameters for Gaussian likelihood function [default]
%      sigma2       - variance [0.1]
%      sigma2_prior - prior for sigma2 [prior_logunif]
%      n            - number of observations per input (See using average
%                     observations below)
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%    The likelihood is defined as follows:
%                    __ n
%      p(y|f1, f2) = || i=1 N(y_i | f1_i, sigma2*exp(f2_i))
%
%      where f1 is the first latent variable defining the mean of the
%      gaussian distribution, f2 is the second latent variable defining
%      the noise structure and sigma2 is coefficient for noise.
%
%  See also
%    GP_SET, LIK_*, PRIOR_*
%
% Copyright (c) 2007-2010 Jarno Vanhatalo & Jouni Hartikainen
% Copyright (c) 2010 Aki Vehtari
% Copyright (c) 2011 Ville Tolvanen

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'LIK_INPUTDEPENDENTNOISE';
  ip.addOptional('lik', [], @isstruct);
  ip.addParamValue('sigma2',0.1, @(x) isscalar(x) && x>0);
  ip.addParamValue('sigma2_prior',prior_logunif(), @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.nondiagW=true;
    lik.type = 'Inputdependentnoise';
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'Inputdependentnoise')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('sigma2',ip.UsingDefaults)
    lik.sigma2 = ip.Results.sigma2;
  end
  % Initialize prior structure
  if init
    lik.p=[];
  end
  if init || ~ismember('sigma2_prior',ip.UsingDefaults)
    lik.p.sigma2=ip.Results.sigma2_prior;
  end

  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_inputdependentnoise_pak;
    lik.fh.unpak = @lik_inputdependentnoise_unpak;
    lik.fh.lp = @lik_inputdependentnoise_lp;
    lik.fh.lpg = @lik_inputdependentnoise_lpg;
    lik.fh.ll = @lik_inputdependentnoise_ll;
    lik.fh.llg = @lik_inputdependentnoise_llg;
    lik.fh.llg2 = @lik_inputdependentnoise_llg2;
    lik.fh.llg3 = @lik_inputdependentnoise_llg3;
    lik.fh.predy = @lik_inputdependentnoise_predy;
    lik.fh.invlink = @lik_inputdependentnoise_invlink;
    lik.fh.predprcty = @lik_inputdependentnoise_predprcty;
    lik.fh.recappend = @lik_inputdependentnoise_recappend;
  end

  function [w,s,h] = lik_inputdependentnoise_pak(lik)
  %LIK_INPUTDEPENDENTNOISE_PAK  Combine likelihood parameters into one vector.
  %
  %  Description
  %    W = LIK_INPUTDEPENDENTNOISE_PAK(LIK) takes a likelihood structure LIK and
  %    combines the parameters into a single row vector W. This is a mandatory
  %    subfunction used for example in energy and gradient computations.
  %
  %       w = log(lik.sigma2)
  %
  %   See also
  %   LIK_INPUTDEPENDENTNOISE_UNPAK, GP_PAK

  w = []; s = {}; h=[];
  if isfield(lik.p, 'sigma2') && ~isempty(lik.p.sigma2)
    w = [w log(lik.sigma2)];
    s = [s 'log(gaussian.sigma2)'];
    h = [h 0];
    % Hyperparameters of sigma2
    [wh, sh, hh] = lik.p.sigma2.fh.pak(lik.p.sigma2);
    w = [w wh];
    s = [s sh];
    h = [h hh];
  end
  end

  function [lik, w] = lik_inputdependentnoise_unpak(lik, w)
  %LIK_INPUTDEPENDENTNOISE_UNPAK  Extract likelihood parameters from the vector.
  %
  %  Description
  %    [LIK, W] = LIK_INPUTDEPENDENTNOISE_UNPAK(W, LIK) takes a likelihood
  %    structure LIK and extracts the parameters from the vector W
  %    to the LIK structure. This is a mandatory subfunction used for
  %    example in energy and gradient computations.
  %
  %   Assignment is inverse of
  %       w = log(lik.sigma2)
  %
  %   See also
  %   LIK_INPUTDEPENDENTNOISE_PAK, GP_UNPAK

  if ~isempty(lik.p.sigma2)
    lik.sigma2 = exp(w(1));
    w = w(2:end);

    % Hyperparameters of sigma2
    [p, w] = lik.p.sigma2.fh.unpak(lik.p.sigma2, w);
    lik.p.sigma2 = p;
  end
  end

  function lp = lik_inputdependentnoise_lp(lik, varargin)
  %LIK_INPUTDEPENDENTNOISE_LP  log(prior) of the likelihood parameters
  %
  %  Description
  %    LP = LIK_INPUTDEPENDENTNOISE_LP(LIK) takes a likelihood structure
  %    LIK and returns log(p(th)), where th collects the parameters. This
  %    subfunction is needed when there are likelihood parameters.
  %
  %  See also
  %    LIK_INPUTDEPENDENTNOISE_LLG, LIK_INPUTDEPENDENTNOISE_LLG3, LIK_INPUTDEPENDENTNOISE_LLG2, GPLA_E

  % If prior for sigma2 parameter, add its contribution
  lp = 0;

  if ~isempty(lik.p.sigma2)
    likp=lik.p;
    lp = likp.sigma2.fh.lp(lik.sigma2, likp.sigma2) + log(lik.sigma2);
  end

  end

  function lpg = lik_inputdependentnoise_lpg(lik)
  %LIK_INPUTDEPENDENTNOISE_LPG  d log(prior)/dth of the likelihood
  %                parameters th
  %
  %  Description
  %    E = LIK_INPUTDEPENDENTNOISE_LPG(LIK) takes a likelihood structure
  %    LIK and returns d log(p(th))/dth, where th collects the parameters.
  %    This subfunction is needed when there are likelihood parameters.
  %
  %  See also
  %    LIK_INPUTDEPENDENTNOISE_LLG, LIK_INPUTDEPENDENTNOISE_LLG3, LIK_INPUTDEPENDENTNOISE_LLG2, GPLA_G

  lpg = [];

  if ~isempty(lik.p.sigma2)
    likp=lik.p;

    lpgs = likp.sigma2.fh.lpg(lik.sigma2, likp.sigma2);
    lpg = lpgs(1).*lik.sigma2 + 1;
    if length(lpgs) > 1
      lpg = [lpg lpgs(2:end)];
    end
  end
  end

  function ll = lik_inputdependentnoise_ll(lik, y, ff, z)
  %LIK_INPUTDEPENDENTNOISE_LL  Log likelihood
  %
  %  Description
  %    LL = LIK_INPUTDEPENDENTNOISE_LL(LIK, Y, F, Z) takes a likelihood
  %    structure LIK, incedence counts Y, expected counts Z, and
  %    latent values F. Returns the log likelihood, log p(y|f,z).
  %    This subfunction is needed when using Laplace approximation
  %    or MCMC for inference with non-Gaussian likelihoods. This
  %    subfunction is also used in information criteria (DIC, WAIC)
  %    computations.
  %
  %  See also
  %    LIK_INPUTDEPENDENTNOISE_LLG, LIK_INPUTDEPENDENTNOISE_LLG3, LIK_INPUTDEPENDENTNOISE_LLG2, GPLA_E

    f=ff(:);
    n=size(y,1);
    f1=f(1:n);
    f2=f((n+1):2*n);
    expf2=exp(f2);
    expf2(expf2<eps)=eps;
    expf2(expf2>1e6)=1e6;
    sigma2 = lik.sigma2;

    ll = sum(-0.5*log(2*pi.*sigma2.*expf2) - 1./(2.*sigma2.*expf2).*(y-f1).^2);

  end

  function llg = lik_inputdependentnoise_llg(lik, y, ff, param, z)
  %LIK_INPUTDEPENDENTNOISE_LLG  Gradient of the log likelihood
  %
  %  Description
  %    LLG = LIK_INPUTDEPENDENTNOISE_LLG(LIK, Y, F, PARAM) takes a likelihood
  %    structure LIK, incedence counts Y, expected counts Z and
  %    latent values F. Returns the gradient of the log likelihood
  %    with respect to PARAM. At the moment PARAM can be 'param' or
  %    'latent'. This subfunction is needed when using Laplace
  %    approximation or MCMC for inference with non-Gaussian likelihoods.
  %
  %  See also
  %    LIK_INPUTDEPENDENTNOISE_LL, LIK_INPUTDEPENDENTNOISE_LLG2, LIK_INPUTDEPENDENTNOISE_LLG3, GPLA_E

    f=ff(:);
    n=size(y,1);
    f1=f(1:n);
    f2=f((n+1):2*n);
    expf2 = exp(f2);
    expf2(expf2<eps)=eps;
    expf2(expf2>1e6)=1e6;
    sigma2 = lik.sigma2;

    switch param
      case 'param'

        llg=sum(-0.5./sigma2+(y-f1).^2./(2*expf2.*sigma2^2));
        % correction for the log transformation
        llg = llg.*lik.sigma2;

      case 'latent'

        llg1= (y-f1)./(sigma2*expf2);
        llg2= -0.5+(y-f1).^2./(2.*(expf2.*sigma2));

%         llg1=(-y+f1)/(sqrt(2*pi).*(sigma2.*expf2).^(3/2)).*exp(-1/(2.*sigma2*expf2).*(y-f1).^2);
%         llg2=-exp(f2-(y-f1).^2/(2.*expf2.*sigma2)).*sigma2/(2.*sqrt(2.*pi).*(expf2.*sigma2).^(3/2))+exp(-f2-(y-f1).^2/(2*expf2*sigma2)).*(y-f1).^2/(2.*sqrt(2.*pi.*expf2).*sigma2.^(3/2));

        llg=[llg1; llg2];
    end
  end

  function [llg2] = lik_inputdependentnoise_llg2(lik, y, ff, param, z)
  %function [pi_vec, pi_mat] = lik_inputdependentnoise_llg2(lik, y, ff, param, z)
  %LIK_INPUTDEPENDENTNOISE_LLG2  Second gradients of the log likelihood
  %
  %  Description
  %    LLG2 = LIK_INPUTDEPENDENTNOISE_LLG2(LIK, Y, F, PARAM) takes a likelihood
  %    structure LIK, incedence counts Y, expected counts Z, and
  %    latent values F. Returns the Hessian of the log likelihood
  %    with respect to PARAM. At the moment PARAM can be only
  %    'latent'. LLG2 is a vector with diagonal elements of the
  %    Hessian matrix (off diagonals are zero). This subfunction
  %    is needed when using Laplace approximation or EP for
  %    inference with non-Gaussian likelihoods.
  %
  %  See also
  %    LIK_INPUTDEPENDENTNOISE_LL, LIK_INPUTDEPENDENTNOISE_LLG, LIK_INPUTDEPENDENTNOISE_LLG3, GPLA_E

    f=ff(:);

    n=length(y);
    f1=f(1:n);
    f2=f((n+1):2*n);
    sigma2 = lik.sigma2;
    expf2=exp(f2);
    expf2(expf2<eps)=eps;
    expf2(expf2>1e6)=1e6;

    switch param
      case 'param'

      case 'latent'

%         llg2 = [2./(sigma2.*expf2); 3/2.*(y-f1).^2./(sigma2.*expf2)];
%         llg2mat = [diag(1./sqrt(sigma2.*expf2)); diag(-(y-f1)./sqrt(sigma2.*expf2))];

        llg2_11=-1./(sigma2.*expf2);
        llg2_12=-(y-f1)./(sigma2*expf2);
        llg2_22=-(y-f1).^2./(2.*sigma2.*expf2);

        llg2 = [llg2_11 llg2_12; llg2_12 llg2_22];

      case 'latent+param'

        llg2_1=-(y-f1)./(expf2.*sigma2^2);
        llg2_2=-(y-f1).^2./(2*sigma2.^2.*expf2);

        llg2=[llg2_1; llg2_2];
        % correction for the log transformation
        llg2 = llg2.*lik.sigma2;

    end
  end

  function llg3 = lik_inputdependentnoise_llg3(lik, y, ff, param, z)
  %LIK_INPUTDEPENDENTNOISE_LLG3  Third gradients of the log likelihood
  %
  %  Description
  %    LLG3 = LIK_INPUTDEPENDENTNOISE_LLG3(LIK, Y, F, PARAM) takes a likelihood
  %    structure LIK, incedence counts Y, expected counts Z and
  %    latent values F and returns the third gradients of the log
  %    likelihood with respect to PARAM. At the moment PARAM can be
  %    only 'latent'. LLG3 is a vector with third gradients. This
  %    subfunction is needed when using Laplace approximation for
  %    inference with non-Gaussian likelihoods.
  %
  %  See also
  %    LIK_INPUTDEPENDENTNOISE_LL, LIK_INPUTDEPENDENTNOISE_LLG, LIK_INPUTDEPENDENTNOISE_LLG2, GPLA_E, GPLA_G

    f=ff(:);

    n=size(y,1);
    f1=f(1:n);
    f2=f((n+1):2*n);
    expf2=exp(f2);
    expf2(expf2<eps)=eps;
    expf2(expf2>1e6)=1e6;
    sigma2 = lik.sigma2;

    switch param
      case 'param'

      case 'latent'
        nl=2;
        llg3=zeros(nl,nl,nl,n);

        % y=0:
        % thrid derivative derivative wrt f1 (11)
        %       llg3(1,1,1,:) = 0
        % thrid derivative derivative wrt f2 (11)
        llg3(1,1,2,:) = 1./(expf2.*sigma2);

        % thrid derivative derivative wrt f1 (12/21)
        llg3(1,2,1,:) = 1./(expf2.*sigma2);
        llg3(2,1,1,:) = llg3(1,2,1,:);
        % thrid derivative derivative wrt f2 (12/21)
        llg3(1,2,2,:) = (y-f1)./(expf2.*sigma2);
        llg3(2,1,2,:) = llg3(1,2,2,:);

        % thrid derivative derivative wrt f1 (22)
        llg3(2,2,1,:) = llg3(1,2,2,:);
        % thrid derivative derivative wrt f1 (22)
        llg3(2,2,2,:) = (y-f1).^2./(2.*expf2.*sigma2);

      case 'latent2+param'

        llg3_11=1./(expf2.*sigma2.^2);
        llg3_12=(y-f1)./(expf2.*sigma2.^2);
        llg3_22=(y-f1).^2./(2.*expf2.*sigma2.^2);

        llg3 = [diag(llg3_11) diag(llg3_12); diag(llg3_12) diag(llg3_22)];
        % correction for the log transformation
        llg3 = llg3.*lik.sigma2;

    end
  end

  function [lpy, Ey, Vary] = lik_inputdependentnoise_predy(lik, Ef, Varf, yt, z)
  %LIK_INPUTDEPENDENTNOISE_PREDY  Returns the predictive mean, variance and density of y
  %
  %  Description
  %    [LPY] = LIK_INPUTDEPENDENTNOISE_PREDY(LIK, EF, VARF, YT) takes a
  %    likelihood structure LIK, posterior mean EF, posterior
  %    variance VARF of the latent variable and observations YT and
  %    returns the logarithm of the predictive density PY of YT, that is
  %        p(yt | th) = \int p(yt | f, th) p(f|y) df.
  %    This subfunction is needed when computing posterior predictive
  %    distributions for future observations.
  %
  %    [LPY, EY, VARY] = LIK_INPUTDEPENDENTNOISE_PREDY(LIK, EF, VARF YT)
  %    Returns also the posterior predictive mean EY and variance VARY of
  %    the observations related to the latent variables. This subfunction
  %    is needed when computing posterior predictive distributions for
  %    future observations.
  %
  %  See also
  %    GPLA_PRED, GPEP_PRED, GPMC_PRED

%     ntest=size(yt,1);
    Ef = Ef(:);
    ntest = 0.5*size(Ef,1);
    Ef1=Ef(1:ntest); Ef2=Ef(ntest+1:end);
    if size(Varf,2) == size(Varf,1)
      Varf1=diag(Varf(1:ntest,1:ntest));Varf2=diag(Varf(ntest+1:end,ntest+1:end));
    else
      if size(Varf,2)==1
        Varf=reshape(Varf,ntest,2);
      end
      Varf1=Varf(:,1); Varf2=Varf(:,2);
    end
    sigma2=lik.sigma2;

    lpy = zeros(size(yt));

    Ey=Ef1;
    Vary=Varf1 + sigma2.*exp(Ef2+Varf2/2);

    if ~isempty(yt)
      for i2=1:ntest
        m1=Ef1(i2); m2=Ef2(i2);
        s1=sqrt(Varf1(i2)); s2=sqrt(Varf2(i2));
        pd=@(f1,f2) norm_pdf(yt(i2), f1, sqrt(sigma2.*exp(f2))).*norm_pdf(f1,Ef1(i2),sqrt(Varf1(i2))).*norm_pdf(f2,Ef2(i2),sqrt(Varf2(i2)));
        lpy(i2) = log(dblquad(pd, m1-6.*s1, m1+6.*s1, m2-6.*s2, m2+6.*s2));
      end
    else
      lpy=[];
    end
%     sigma2 = lik.sigma2;
%     Ef1=Ef(1:ntest);
%     Ef2=Ef((ntest+1):2*ntest);
%     Ey = Ef1;
%     Vary = Varf + sigma2.*exp(Ef2);

  end

  function prctys = lik_inputdependentnoise_predprcty(lik, Ef, Varf, zt, prcty)
  %LIK_BINOMIAL_PREDPRCTY  Returns the percentiles of predictive density of y
  %
  %  Description
  %    PRCTY = LIK_BINOMIAL_PREDPRCTY(LIK, EF, VARF YT, ZT)
  %    Returns percentiles of the predictive density PY of YT, that is
  %    This requires also the succes counts YT, numbers of trials ZT. This
  %    subfunction is needed when using function gp_predprcty.
  %
  %  See also
  %    GP_PREDPCTY

  n=size(Ef,1)./2;
  prcty = prcty./100;
  prcty = norminv(prcty, 0, 1);
  prctys = bsxfun(@plus, Ef(1:n), bsxfun(@times, sqrt(Varf(1:n) + lik.sigma2.*exp(Ef(n+1:end))), prcty));

  end

  function p = lik_inputdependentnoise_invlink(lik, f, z)
  %LIK_INPUTDEPENDENTNOISE_INVLINK  Returns values of inverse link function
  %
  %  Description
  %    P = LIK_INPUTDEPENDENTNOISE_INVLINK(LIK, F) takes a likelihood structure LIK and
  %    latent values F and returns the values of inverse link function P.
  %    This subfunction is needed when using function gp_predprcty.
  %
  %     See also
  %     LIK_INPUTDEPENDENTNOISE_LL, LIK_INPUTDEPENDENTNOISE_PREDY

    p = exp(f);
  end

  function reclik = lik_inputdependentnoise_recappend(reclik, ri, lik)
  %RECAPPEND  Append the parameters to the record
  %
  %  Description
  %    RECLIK = GPCF_INPUTDEPENDENTNOISE_RECAPPEND(RECLIK, RI, LIK) takes a
  %    likelihood record structure RECLIK, record index RI and
  %    likelihood structure LIK with the current MCMC samples of
  %    the parameters. Returns RECLIK which contains all the old
  %    samples and the current samples from LIK. This subfunction
  %    is needed when using MCMC sampling (gp_mc).
  %
  %  See also
  %    GP_MC

  % Initialize record
    if nargin == 2
      reclik.type = 'Inputdependentnoise';
      reclik.nondiagW=true;

      % Initialize parameter

      % Set the function handles
      reclik.fh.pak = @lik_inputdependentnoise_pak;
      reclik.fh.unpak = @lik_inputdependentnoise_unpak;
      reclik.fh.lp = @lik_inputdependentnoise_lp;
      reclik.fh.lpg = @lik_inputdependentnoise_lpg;
      reclik.fh.ll = @lik_inputdependentnoise_ll;
      reclik.fh.llg = @lik_inputdependentnoise_llg;
      reclik.fh.llg2 = @lik_inputdependentnoise_llg2;
      reclik.fh.llg3 = @lik_inputdependentnoise_llg3;
      reclik.fh.predy = @lik_inputdependentnoise_predy;
      reclik.fh.predprcty = @lik_inputdependentnoise_predprcty;
      reclik.fh.invlink = @lik_inputdependentnoise_invlink;
      reclik.fh.recappend = @lik_inputdependentnoise_recappend;
      reclik.p=[];
      if ~isempty(ri.p.sigma2)
        reclik.p.sigma2 = ri.p.sigma2;
      end
      return
    end

    reclik.sigma2(ri,:)=lik.sigma2;
    if ~isempty(lik.p.sigma2)
      reclik.p.sigma2 = feval(lik.p.sigma2.fh.recappend, reclik.p.sigma2, ri, lik.p.sigma2);
    end
  end
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function beta_values_test ( )

%*****************************************************************************80
%
%% BETA_VALUES_TEST demonstrates the use of BETA_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BETA_VALUES_TEST:\n' );
  fprintf ( 1, '  BETA_VALUES stores values of\n' );
  fprintf ( 1, '  the Beta function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X              Y         BETA(X,Y)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, y, fxy ] = beta_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %12f  %24.16f\n', x, y, fxy );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [Iphi,SigmaPhi,deltaMuPhi,suffStat] = VBA_Iphi_UNL(phi,y,posterior,suffStat,dim,u,options)
% Gauss-Newton update of the obs. parameters, for un-normalized likelihoods

if options.DisplayWin % Display progress
    try
        STR = 'VB Gauss-Newton on observation parameters... ';
        set(options.display.hm(1),'string',STR);
        set(options.display.hm(2),'string','0%');
        drawnow
    end
end

%  Look-up which evolution parameter to update
indIn = options.params2update.phi;

% Preallocate intermediate variables
Q = options.priors.SigmaPhi(indIn,indIn);
iQ = VBA_inv(Q);
muPhi0 = options.priors.muPhi;
Phi = muPhi0;
Phi(indIn) = phi;
dphi0 = muPhi0-Phi;

% inverse temperature
beta = posterior.a_sigma./posterior.b_sigma;

LL = 0;
dLLdP = zeros(1,options.dim.n_phi);
d2LLdP2 = zeros(options.dim.n_phi,options.dim.n_phi);
Ey = zeros(options.dim.p,options.dim.n_t);
Vy = zeros(options.dim.p,options.dim.n_t);
div = 0;

%--- Loop over time series ---%
for t=1:dim.n_t

    if ~options.isYout

        % evaluate re-normalized likelihood (as well as gradients and Hessians)
        [LLt,dLLdXt,dLLdPt,d2LLdX2t,d2LLdP2t,Ey(:,t),Vy(:,t)] = VBA_evalAL([],Phi,beta,u(:,t),y(:,t),options);

        LL = LL + LLt;
        dLLdP = dLLdP + dLLdPt;
        d2LLdP2 = d2LLdP2 + d2LLdP2t;

        % Accelerate divergent update
        if VBA_isWeird ({LL, dLLdP, d2LLdP2, Ey, Vy})
            div = 1;
            break
        end

    end

    % Display progress
    if mod(t,dim.n_t./10) < 1
        if  options.DisplayWin
            try
                set(options.display.hm(2),'string',[num2str(floor(100*t/dim.n_t)),'%']);
                drawnow
            end
        end
    end

    % Check gradients?
    if options.checkGrads
        mayPause = 0;
        if ~isempty(Phi)
            dLLdPt2 = VBA_numericDiff(@VBA_evalAL,2,[],Phi,beta,u(:,t),y(:,t),options);
            if ~ VBA_isWeird (dLLdPt2)
                [hf] = VBA_displayGrads(dLLdPt',dLLdPt2,'Gradients wrt parameters',options.g_fname,'g');
                mayPause = 1;
            else
                VBA_disp('VBA check_grads: Warning: weird numerical gradients!!!')
            end
            d2LLdP2t2 = VBA_numericDiff(@numericDiff,4,@VBA_evalAL,2,[],Phi,beta,u(:,t),y(:,t),options);
            if ~ VBA_isWeird (dLLdPt2)
                [hf] = VBA_displayGrads(d2LLdP2t,d2LLdP2t2,'Hessians wrt parameters',options.g_fname,'g');
                mayPause = 1;
            else
                VBA_disp('VBA check_grads: Warning: weird numerical gradients!!!')
            end
        end
        if mayPause
            pause
            close(setdiff(hf,0))
        end
    end

end

% Display progress
if options.DisplayWin
    try
        set(options.display.hm(2),'string','OK');
        drawnow
    end
end

% posterior covariance matrix
iSigmaPhi = iQ - d2LLdP2(indIn,indIn);
SigmaPhi = VBA_inv(iSigmaPhi);

% mode
tmp = iQ*dphi0(indIn) + dLLdP(indIn)';
deltaMuPhi = SigmaPhi*tmp;

% variational energy
Iphi = -0.5.*dphi0(indIn)'*iQ*dphi0(indIn) + LL;
if VBA_isWeird ({Iphi, SigmaPhi})
    Iphi = -Inf;
end

% update sufficient statistics
suffStat.logL = LL;
suffStat.gx = Ey; % model post-dicted observations
suffStat.dy = y - suffStat.gx;
suffStat.vy = Vy;
suffStat.dphi = dphi0;
suffStat.div = div;

author: MBB-team

repo: VBA-toolbox


text:
function tf = spm_cfg_eeg_tf
% Configuration file for M/EEG time-frequency analysis
%__________________________________________________________________________
% Copyright (C) 2010-2011 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_cfg_eeg_tf.m 6929 2016-11-14 13:07:31Z guillaume $

%--------------------------------------------------------------------------
% D
%--------------------------------------------------------------------------
D        = cfg_files;
D.tag    = 'D';
D.name   = 'File Name';
D.filter = 'mat';
D.num    = [1 1];
D.help   = {'Select the M/EEG mat file.'};

%--------------------------------------------------------------------------
% timewin
%--------------------------------------------------------------------------
timewin         = cfg_entry;
timewin.tag     = 'timewin';
timewin.name    = 'Time window';
timewin.strtype = 'r';
timewin.num     = [1 2];
timewin.val     = {[-Inf Inf]};
timewin.help    = {'Time window (ms)'};

%--------------------------------------------------------------------------
% frequencies
%--------------------------------------------------------------------------
frequencies         = cfg_entry;
frequencies.tag     = 'frequencies';
frequencies.name    = 'Frequencies of interest';
frequencies.strtype = 'r';
frequencies.num     = [0 Inf];
frequencies.val     = {[]};
frequencies.help    = {'Frequencies of interest (as a vector), if empty 1-48 with optimal frequency bins ~1 Hz or above resolution'};

%--------------------------------------------------------------------------
% phase
%--------------------------------------------------------------------------
phase        = cfg_menu;
phase.tag    = 'phase';
phase.name   = 'Save phase';
phase.help   = {'Save phase as well as power'};
phase.labels = {'yes', 'no'};
phase.values = {1, 0};
phase.val    = {0};

%--------------------------------------------------------------------------
% method
%--------------------------------------------------------------------------
method      = cfg_choice;
method.tag  = 'method';
method.name = 'Spectral estimation';
method.help = {'Spectral estimation'};

specest_funs = spm_select('List',spm('dir'),'^spm_eeg_specest_.*\.m$');
specest_funs = cellstr(specest_funs);
method.values = cell(1,numel(specest_funs));
for i = 1:numel(specest_funs)
    method.values{i} = feval(spm_file(specest_funs{i},'basename'));
end

%--------------------------------------------------------------------------
% prefix
%--------------------------------------------------------------------------
prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the output dataset.'};
prefix.strtype = 's';
prefix.num     = [0 Inf];
prefix.val     = {''};

%--------------------------------------------------------------------------
% M/EEG Time-Frequency Analysis
%--------------------------------------------------------------------------
tf = cfg_exbranch;
tf.tag = 'tf';
tf.name = 'Time-frequency analysis';
tf.val = {D, spm_cfg_eeg_channel_selector, frequencies, timewin, method, phase, prefix};
tf.help = {'Perform time-frequency analysis of epoched M/EEG data.'};
tf.prog = @eeg_tf;
tf.vout = @vout_eeg_tf;
tf.modality = {'EEG'};

%==========================================================================
% function out = eeg_tf(job)
%==========================================================================
function out = eeg_tf(job)
% construct the S struct
S   = [];
S.D = job.D{1};

S.channels = spm_cfg_eeg_channel_selector(job.channels);

S.frequencies = job.frequencies;
S.timewin = job.timewin;
S.phase = job.phase;

S.method = cell2mat(fieldnames(job.method));
S.settings = job.method.(S.method);

S.prefix = job.prefix;

[Dtf, Dtph] = spm_eeg_tf(S);

out.Dtf = Dtf;
out.Dtfname = {Dtf.fullfile};

out.Dtph = Dtph;
if ~isempty(Dtph)
    out.Dtphname = {Dtph.fullfile};
else
    out.Dtphname = {''};
end

%==========================================================================
% function dep = vout_eeg_tf(job)
%==========================================================================
function dep = vout_eeg_tf(job)
% return dependencies
dep(1)            = cfg_dep;
dep(1).sname      = 'M/EEG time-frequency power dataset';
dep(1).src_output = substruct('.','Dtf');
% this can be entered into any evaluated input
dep(1).tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2)            = cfg_dep;
dep(2).sname      = 'M/EEG time-frequency power dataset';
dep(2).src_output = substruct('.','Dtfname');
% this can be entered into any file selector
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});

if job.phase
    dep(3)            = cfg_dep;
    dep(3).sname      = 'M/EEG time-frequency phase dataset';
    dep(3).src_output = substruct('.','Dtph');
    % this can be entered into any evaluated input
    dep(3).tgt_spec   = cfg_findspec({{'strtype','e'}});

    dep(4)            = cfg_dep;
    dep(4).sname      = 'M/EEG time-frequency phase dataset';
    dep(4).src_output = substruct('.','Dtphname');
    % this can be entered into any file selector
    dep(4).tgt_spec   = cfg_findspec({{'filter','mat'}});
end

author: spm


text:
function varargout = openopt_fun(varargin)

persistent params

if nargin>1
    params = varargin{2};
    return
end

if ischar(varargin{1})
    switch varargin{1}
        case 'getStartPoint'
            varargout{1} = params.x0;
            return
        case 'getOptimPoint'
            varargout{1} =[];
            return
    end
end

xevaled = zeros(1,length(params.c));
xevaled(params.linearindicies) = varargin{1};

% Apply the precomputed evaluation scheme (if necessary)
xevaled = apply_recursive_evaluation(params,xevaled);

xevaled = xevaled(:);
varargout{1} = params.f + (params.c'+xevaled'*params.Q)*xevaled;


text:
function hail_test ( )

%*****************************************************************************80
%
%% HAIL_TEST tests HAIL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HAIL_TEST\n' );
  fprintf ( 1, '  HAIL(I) computes the length of the hail sequence\n' );
  fprintf ( 1, '  for I, also known as the 3*N+1 sequence.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I,  HAIL(I)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 20
    fprintf ( 1, '%2d  %6d\n', i, hail(i) );
  end

  return
end


text:
function [ a, b ] = p01_lim ( )

%*****************************************************************************80
%
%% P01_LIM returns the integration limits for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A(2), B(2), the lower and upper limits of integration.
%
  a(1:2) = 0.0;
  b(1:2) = 1.0;

  return
end


text:
%  Figure 9.54      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
function y = fas(u)
if u < 0.0606,
   y = 0.1 ;
elseif u < 0.0741,
   y = 0.1 + (u-0.0606)*20;
else y = 0.9;
end

author: Sable

repo: mcbench-benchmarks


text:
function [time,data,nEpochs] = readBESAsb(filename)

% readBESAsb reads information from a *.dat, i.e. a simple binary data file
%
% This function requires a file with the same basename as the data name but
% suffix '.generic' or '.gen' to be located in the same folder. This file is
% generated automatically by BESA during file export.
%
% Use as
%   [time,data,nEpochs] = readBESAsb(filename)
%
% The following output is generated:
%   time: a vector containing the time points corresponding to the data
%         points. Negative times are prestimulus
%   data: The data matrix with dimension [nChannels x nEpochs x nSamples],
%         where nChannels is the number of Channels and nSamples the number
%         of samples within one epoch
%   nEpochs (optional): The number of epochs contained in the file

% Modified February 21, 2007 Karsten Hoechstetter
% Modified April 24, 2007 Robert Oostenveld
% Modified September 24, 2009 Karsten Hoechstetter

if isempty(findstr(filename,'.dat'))
  filename = [filename,'.dat'];
end

fid=fopen([filename(1:end-4),'.generic'],'r');
if fid==-1
  fid=fopen([filename(1:end-4),'.gen'],'r');
end

fscanf(fid,'BESA Generic Data\n');
nChannels = fscanf(fid,'nChannels=%i\n');
sRate = fscanf(fid,'sRate=%f\n');
nSamples = fscanf(fid,'nSamples=%i\n');
format = fscanf(fid,'format=%s');
file = fscanf(fid,'\nfile=%s');
prestimulus = fscanf(fid,'prestimulus=%f\n');   if isempty(prestimulus),prestimulus=0;end
epochs = fscanf(fid,'epochs=%i\n');             if isempty(epochs),epochs=1;end
fclose(fid);

time=[-prestimulus:1/sRate*1000:(nSamples/epochs-1)*1000/sRate-prestimulus];

fid=fopen(filename, 'r', 'ieee-le');
xdata=fread(fid,[nChannels,nSamples],'float32');
fclose(fid);

data=zeros(nChannels,epochs,nSamples/epochs);
for i=1:epochs
  data(:,i,:)=xdata(:,1+(i-1)*nSamples/epochs:i*nSamples/epochs);
end
nEpochs=epochs;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% MCMCACF - autocorrelation plots
% Copyright (c) 1998, Harvard University. Full copyright in the file Copyright
%
% This is not a function, but a way to get acf plots for MCMC runs.
%
% S = mcmcsumm(A) contains S.acf, a matrix of the autocorrelations.
%
% So, mcmctrace(S.acf) will plot the autocorrelations.
%
% See also: MCMCSUMM, MCMCTRACE


text:
mex -lopencv_core -lopencv_imgproc -L/usr/local/lib -I/usr/local/include/ mexResize.cpp MxArray.cpp
mex gradientMex.cpp

author: flyers

repo: drone-tracking


text:
function [d,delay_offset] = driving_function_imp_wfs_pwd(x0,ppwd,xq,conf)
%DRIVING_FUNCTION_IMP_WFS_PWD driving signal for a plane wave expansion in WFS
%
%   Usage: [d,delay_offset] = driving_function_imp_wfs_pwd(x0,ppwd,[xq],conf)
%
%   Input parameters:
%       x0          - position and direction of the secondary source / m [N0x6]
%       ppwd        - plane wave coefficients [N x Npw]
%       xq          - centre of plane wave expansion, default = [0,0,0];
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       d             - driving function signals [N x N0]
%       delay_offset  - additional added delay, so you can correct it
%
%   See also: driving_function_imp_wfs_vss, driving_function_imp_localwfs_sbl

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

%% ===== Checking of input  parameters ==================================
nargmin = 3;
nargmax = 4;
narginchk(nargmin,nargmax);
isargsecondarysource(x0);
isargmatrix(ppwd);
if nargin == nargmin
    conf = xq;
    xq = [0, 0, 0];
else
    isargxs(xq);
end
isargstruct(conf);

%% ===== Computation ====================================================
% create distribution of plane waves as virtual secondary sources
Npw = size(ppwd, 2);
phipw = (0:Npw-1).'*2*pi/Npw;
xv = [cos(phipw), sin(phipw)];  % [Npw x 2]
xv(:,3) = 0;
xv(:,4:6) = xv(:,1:3);
xv(:,7) = 1./Npw;  % apply integrations weights to pwd

% shift coordinates to expansion center
x0(:,1:3) = bsxfun(@minus,x0(:,1:3),xq);
conf.xref = [0 0 0];

[d,delay_offset] = driving_function_imp_wfs_vss(x0,xv,'pw',ppwd,conf);

author: sfstoolbox

repo: sfs-matlab


text:
function [Y,Cb,Cr]=RGB2YCbCr(img_rgb)
R=img_rgb(:,:,1);
G=img_rgb(:,:,2);
B=img_rgb(:,:,3);

% rgbתihs
Y   = 0.257*R+0.564*G+0.098*B+16;
Cb = -0.148*R-0.291*G+0.439*B+128;
Cr  = 0.439*R-0.368*G-0.071*B+128;

end

author: thfylsty


text:
function [vs, output, v0] = computeConfidenceInterval2(v0, expdata, model, max_score)

majorIterationLimit = 2000;  %max number of iterations
minorIterationLimit = 1e7;   % essentially infinity
diffInterval = 1e-5;         %gradient step size.
feasibilityTolerance = max_score/20; % how close you need to be to the max score.

v0 = fitC13Data(v0,expdata,model);

if ~isfield(model, 'N')
   model.N = null(model.S);
end

x0 = model.N\v0; % back substitute

% safety check:
if (max(abs(model.S*v0))> 1e-6)
    display('v0 not quite in null space');
    pause;
end
if(max(abs(model.N*x0 - v0)) > 1e-6)
    display('null basis is weird');
    pause;
end

nalpha = size(model.N, 2);
x_L = -1000*ones(nalpha,1);
x_U = 1000*ones(nalpha,1);
Name = 't2';
[A, b_L, b_U] = defineLinearConstraints(model);

numpoints = size(x0,2);
scores = zeros(numpoints,1);

% compute scores for all points.
tProb.user.expdata = expdata;
tProb.user.model = model;
for i = 1:numpoints
    scores(i) = errorComputation2(x0(:,i),tProb);
end
valid_index = scores < max_score + feasibilityTolerance;
fprintf('found %d valid points\n', sum(valid_index));

x0_valid = x0(:,valid_index);
x0_invalid = x0;
scores_valid = scores(valid_index);

c = []; c_L = []; c_U = []; HessPattern = [];
f = 'errorComputation2';
g = 'errorComputation2_grad'; H = [];
%c = 'errorComputation2';
%c_L = 0;
%c_U = max_score;

dc = []; d2c = []; ConsPattern = [];
pSepFunc = [];
x_min = []; x_max = []; f_opt = []; x_opt = [];
Solver = 'snopt';

% pre-compute unnecesary directions.
% if checkedbefore(i) ~= 0 then direction i is redundant
% checkedbefore(i) < 0 means that the sign of all computations should be
% switched.
checkedbefore = zeros(length(model.lb),1);
for i = 2:length(model.lb)
    d = objective_coefficient(i, model);
    for j = 1:i-1
        dj = objective_coefficient(j, model);
        if max(abs(dj - d))< 1e-4
            checkedbefore(i) = j;
            break;
        elseif max(abs(dj + d))< 1e-4
            checkedbefore(i) = -j;
            break;
        end
    end
end

% initialize variables;
numpoints = size(x0, 2);

outputminv = 222*ones(length(model.lb), numpoints);
outputminexitflag = -222*ones(length(model.lb), numpoints);
outputminfinalscore = -222*ones(length(model.lb), numpoints);
outputmaxv = -222*ones(length(model.lb), numpoints);
outputmaxexitflag = -222*ones(length(model.lb), numpoints);
outputmaxfinalscore = -222*ones(length(model.lb), numpoints);
outputminstruct = cell(length(model.lb), numpoints);
outputmaxstruct = cell(length(model.lb), numpoints);

%iterate through directions
parfor i = 1:length(model.lb)
    if checkedbefore(i) ~= 0
        continue;
    end
    for j = -1:2:1 % max and min
        d = objective_coefficient(i,model)*j;
        lbprob = lpAssign(d, A, b_L, b_U, x_L, x_U, [], 'linear_bound', [],[],[],[],[],[],[]);
        lbresult = tomRun('cplex', lbprob, 0);
        fLowBnd = lbresult.f_k;
        fprintf('reaction %d of %d, direction %d, lowerbound %f\n', i,length(model.lb), j, fLowBnd);

        % short circuit if x0 already close to a bound.
        obj1 = d'*x0_valid;
        if(any(abs(obj1-fLowBnd)<.0001))
            display('short circuiting');
            [nil, index1] = min(obj1);
            if (j > 0)
                outputminv(i,:) = j*fLowBnd; % multiply by j to correct sign.
                outputminexitflag(i,:) = 111;
                outputminfinalscore(i,:) = scores_valid(index1);
            else
                outputmaxv(i,:) = j*fLowBnd; % multiply by j to correct sign.
                outputmaxexitflag(i,:) = 111;
                outputmaxfinalscore(i,:) = scores_valid(index1);
            end
        else % gotta actually do the computation.
            all_ds = d'*x0_invalid;
            [nil, index2] = sort(all_ds);

            % initialize temp variables to make parfor work.
            v = j*222*ones(1,numpoints);
            exitflag = -222*ones(1,numpoints);
            finalscore = 222*ones(1,numpoints);
            ostruct = cell(1,numpoints);

            for k = 1:length(index2)
                if exist('ttt.txt', 'file')
                    fprintf('quitting due to file found\n');
                    continue;
                end
                xinitial = x0_invalid(:,index2(k));
%                 Prob  = lpconAssign(d, x_L, x_U, Name, xinitial,...
%                                       A, b_L, b_U,...
%                                       c, dc, d2c, ConsPattern, c_L, c_U,...
%                                       fLowBnd, x_min, x_max, f_opt, x_opt);
                Prob = conAssign(f, g, H, HessPattern, x_L, x_U, Name, xinitial, ...
                                pSepFunc, fLowBnd, ...
                                A, b_L, b_U, c, dc, d2c, ConsPattern, c_L, c_U, ...
                                x_min, x_max, f_opt, x_opt);
                %Prob.NumDiff = 2; % central diff
                %Prob.optParam.CentralDiff = 1e-5;
                %pause;
                Prob.user.expdata = expdata;
                Prob.user.model = model;
                Prob.user.objective = d;
                Prob.user.max_error = max_score;
                Prob.user.diff_interval = diffInterval;
                Prob.user.multiplier = 10;

                Prob.optParam.IterPrint = 0;
                %Prob.optParam.cTol = .1*feasibilityTolerance;

                Prob.PriLevOpt = 0;
                Prob.SOL.PrintFile = strcat('temp/snoptp', num2str(i), 'x', num2str(j), 'x', num2str(k),'.txt');
                Prob.SOL.SummFile = strcat('temp/snopts', num2str(i), 'x', num2str(j), 'x', num2str(k),'.txt');
                Prob.SOL.optPar(35) = majorIterationLimit; %This is major iteration count.
                Prob.SOL.optPar(30) = minorIterationLimit; %total iteration limit;
                %Prob.SOL.optPar(11) = feasibilityTolerance; % feasibility tolerance

                Result = tomRun(Solver, Prob, 5);
                tscore = errorComputation2(Result.x_k, tProb);
                tbest = Result.f_k;

                fprintf('reaction %d (%d), x %d; x=%f (%f); score=%f (%f)\n', i,length(model.lb),j, tbest,fLowBnd, tscore, max_score)

                v(k) = j*tbest;
                exitflag(k) = Result.Inform;
                finalscore(k) = tscore;
                ostruct{k} = Result;
            end
            if (j > 0) %minimizing
                outputminv(i,:) = v; % multiply by j to correct sign.
                outputminexitflag(i,:) = exitflag;
                outputminfinalscore(i,:) = finalscore;
                outputminstruct(i,:) = ostruct;
            else
                outputmaxv(i,:) = v; % multiply by j to correct sign.
                outputmaxexitflag(i,:) = exitflag;
                outputmaxfinalscore(i,:) = finalscore;
                outputmaxstruct(i,:) = ostruct;
            end
        end
    end
end

output.minv = outputminv;
output.maxv = outputmaxv;
output.minexitflag = outputminexitflag;
output.maxexitflag = outputmaxexitflag;
output.minfinalscore = outputminfinalscore;
output.maxfinalscore = outputmaxfinalscore;
output.minstruct = outputminstruct;
output.maxstruct = outputmaxstruct;

for i = 1:length(model.lb)
    if checkedbefore(i) > 0 %short circuit if seen before.
        output.minv(i,:) = output.minv(checkedbefore(i),:);
        output.maxv(i,:) = output.maxv(checkedbefore(i),:);
        output.minexitflag(i,:) = output.minexitflag(checkedbefore(i),:);
        output.maxexitflag(i,:) = output.maxexitflag(checkedbefore(i),:);
        output.minfinalscore(i,:) = output.minfinalscore(checkedbefore(i),:);
        output.maxfinalscore(i,:) = output.maxfinalscore(checkedbefore(i),:);
        output.minstruct(i,:) = output.minstruct(checkedbefore(i),:);
        output.maxstruct(i,:) = output.maxstruct(checkedbefore(i),:);
    elseif checkedbefore(i) < 0
        output.minv(i,:) = -output.maxv(-checkedbefore(i),:);
        output.maxv(i,:) = -output.minv(-checkedbefore(i),:);
        output.minexitflag(i,:) = output.maxexitflag(-checkedbefore(i),:);
        output.maxexitflag(i,:) = output.minexitflag(-checkedbefore(i),:);
        output.minfinalscore(i,:) = output.maxfinalscore(-checkedbefore(i),:);
        output.maxfinalscore(i,:) = output.minfinalscore(-checkedbefore(i),:);
        output.minstruct(i,:) = output.maxstruct(-checkedbefore(i),:);
        output.maxstruct(i,:) = output.minstruct(-checkedbefore(i),:);
    end
end

vs = zeros(length(model.lb), 2);
for i = 1:length(model.lb)
    validindex = output.minfinalscore(i,:) < max_score + feasibilityTolerance;
    if any(validindex)
        vs(i,1) = min(output.minv(i,validindex));
    else
        vs(i,1) = 222;
    end
    validindex = output.maxfinalscore(i,:) < max_score + feasibilityTolerance;
    if any(validindex)
        vs(i,2) = max(output.maxv(i,validindex));
    else
        vs(i,2) = -222;
    end
end

return;

% function that returns the proper objective coefficient for each reaction
% takes into account the reversibility of reactinos etc.
function [d] =  objective_coefficient(i, model)
d = zeros(length(model.lb),1);
d(i) = 1;
if (model.match(i))
    d(model.match(i)) = -1;
end
d = (d'*model.N)'; % transform to null space;
return

author: opencobra

repo: cobratoolbox


text:
% DIAGPROD - Computes the diagonal elements of a matrix product.
%
% D = DIAGPROD(A,B)
%
% Evaluates diag(A*B') efficiently. A and B must have the same size.
%
% D = DIAGPROD(A)
%
% Evaluates diag(A*A') efficiently.

% Last modified 2010-06-24
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function d = diagprod(A,B)

if nargin < 2
  d = dot(A,A,2);
else
  d = dot(A,B,2);
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function wrench = fromParametrizationToForces(xi, Config)

    % delta feet size
    delta_x   =  (Config.feet_size(2,2) - Config.feet_size(2,1))/2;
    delta_x_0 = -(Config.feet_size(2,2) + Config.feet_size(2,1))/2;
    delta_y   =  (Config.feet_size(1,2) - Config.feet_size(1,1))/2;
    delta_y_0 =  (Config.feet_size(1,2) + Config.feet_size(1,1))/2;

    % linear forces
    Fz = exp(xi(3)) + Config.fZmin;
    Fx = sqrt(Config.forceFrictionCoefficient)*tanh(xi(1))*Fz/sqrt(1 + tanh(xi(2))^2);
    Fy = sqrt(Config.forceFrictionCoefficient)*tanh(xi(2))*Fz/sqrt(1 + tanh(xi(1))^2);

    % moments
    Mz = Config.torsionalFrictionCoefficient*tanh(xi(6))*Fz;
    Mx = (delta_x*tanh(xi(4)) + delta_x_0)*Fz;
    My = (delta_y*tanh(xi(5)) + delta_y_0)*Fz;

    % compute the final wrench
    wrench = [Fx; Fy; Fz; Mx; My; Mz];
end

author: robotology

repo: whole-body-controllers


text:
function [AF,BF]=gabframebounds(g,a,M,varargin)
%GABFRAMEBOUNDS  Calculate frame bounds of Gabor frame
%   Usage:  fcond=gabframebounds(g,a,M);
%           [A,B]=gabframebounds(g,a,M);
%           [A,B]=gabframebounds(g,a,M,L);
%           [A,B]=gabframebounds(g,a,M,'lt',lt);
%
%   Input parameters:
%           g     : The window function.
%           a     : Length of time shift.
%           M     : Number of channels.
%           L     : Length of transform to consider.
%           lt    : Lattice type (for non-separable lattices).
%   Output parameters:
%           fcond : Frame condition number (B/A)
%           A,B   : Frame bounds.
%
%   `gabframebounds(g,a,M)` calculates the ratio $B/A$ of the frame bounds
%   of the Gabor system with window *g*, and parameters *a*, *M*.
%
%   `[A,B]=gabframebounds(...)` returns the frame bounds *A* and *B*
%   instead of just the ratio.
%
%   The window *g* may be a vector of numerical values, a text string or a
%   cell array. See the help of |gabwin| for more details.
%
%   `gabframebounds(g,a,M,L)` will cut or zero-extend the window to length
%   *L*.
%
%   `gabframebounds(g,a,M,'lt',lt)` does the same for a non-separable
%   lattice specified by *lt*. Please see the help of |matrix2latticetype|
%   for a precise description of the parameter *lt*.
%
%   See also: gabrieszbounds, gabwin

%% ---------- Assert correct input.

if nargin<3
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.keyvals.L=[];
definput.keyvals.lt=[0 1];
[flags,kv,L]=ltfatarghelper({'L'},definput,varargin);

%% ------ step 2: Verify a, M and L
if isempty(L)
    if isnumeric(g)
        % Use the window length
        Ls=length(g);
    else
        % Use the smallest possible length
        Ls=1;
    end;

    % ----- step 2b : Verify a, M and get L from the window length ----------
    L=dgtlength(Ls,a,M,kv.lt);

else

    % ----- step 2a : Verify a, M and get L

    Luser=dgtlength(L,a,M,kv.lt);
    if Luser~=L
        error(['%s: Incorrect transform length L=%i specified. Next valid length ' ...
               'is L=%i. See the help of DGTLENGTH for the requirements.'],...
              upper(mfilename),L,Luser);
    end;

end;

%% ----- step 3 : Determine the window

[g,info]=gabwin(g,a,M,L,kv.lt,'callfun',upper(mfilename));

if L<info.gl
  error('%s: Window is too long.',upper(mfilename));
end;

%% ----- actual computation ------------

g=fir2long(g,L);
R=size(g,2);

if kv.lt(2)==1
    % Rectangular case
    % Get the factorization of the window.
    gf=comp_wfac(g,a,M);

    % Compute all eigenvalues.
    lambdas=comp_gfeigs(gf,L,a,M);
    s=size(lambdas,1);

else

    % Convert to multi-window
    mwin=comp_nonsepwin2multi(g,a,M,kv.lt,L);

    % Get the factorization of the window.
    gf=comp_wfac(mwin,a*kv.lt(2),M);

    % Compute all eigenvalues.
    lambdas=comp_gfeigs(gf,L,a*kv.lt(2),M);
    s=size(lambdas,1);

end;

% Min and max eigenvalue.
if a>M*R
    % This can is not a frame, so A is identically 0.
    AF=0;
else
    AF=lambdas(1);
end;

BF=lambdas(s);

if nargout<2
    % Avoid the potential warning about division by zero.
    if AF==0
        AF=Inf;
    else
        AF=BF/AF;
    end;
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%  Figure 9.47      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
%% fig9.47.m
%% Simulink Simulation for Position Control
clf;
[t,x,y]=sim('fig9_46');
plot(t',x');
grid;
xlabel('Time (sec)');
ylabel('Amplitude');
title('Fig. 9.47: Position Control');
text(5,0.95,'Output');
text(3.5,0.3,'x_2');


text:
function M = rmPlotGUI_getModel(view, roi, allTimePoints, preserveCoords)
% retrieve the model data for the specified retinotopy model, in a struct M.
%
%  M = rmPlotGUI_getModel(view, roi, allTimePoints, [preserveCoords])
%
% The M structure will reside as the GUI's UserData, and contain a compact
% description of both the model and the current state of the GUI. I broke
% off the accessor function, in case anyone wants to try to do command-line
% analyses.
%
% ras, 11/2008; broken off into its own function.
if notDefined('preserveCoords'), preserveCoords = 1;                end

% load file with data
rmFile = viewGet(view, 'rmFile');
if isempty(rmFile),
    fprintf('[%s]:No file selected\n', mfilename);
    return;
end;

% load model
try
    model = viewGet(view, 'rmModel');
catch
    model = [];
end;
if isempty(model),
    load(viewGet(view, 'rmFile'), 'model');
    view  = viewSet(view, 'rmModel', model);
    model = viewGet(view, 'rmModel');
end;

% load params
try
    params = viewGet(view, 'rmParams');
catch
    params = [];
end;
params = rmRecomputeParams(view, params, allTimePoints);

modelNames = viewGet(view, 'rmModelNames');

% initialize the M struct with obvious fields
M.roi = roi;
M.model    = model;
M.modelNum = 1;
M.params = params;
M.modelNames = modelNames;
M.dataType = getDataTypeName(view);
M.viewType = view.viewType;

M.prevVoxel = 1;  % this keeps track of the last selected voxel, to detect changes

% get time series and roi-coords
[M.tSeries, M.coords, M.params] = rmLoadTSeries(view, params, roi, preserveCoords);

% detrend
% get/make trends
verbose = false;
trends  = rmMakeTrends(params);
%if isfield(M.params.analysis,'allnuisance')
%    trends = [trends M.params.analysis.allnuisance];
%end

% recompute
b = pinv(trends)*M.tSeries;
M.tSeries = M.tSeries - trends*b;
nt = size(trends,2);

% make x-axis for time series plot
nFramesInd = [M.params.stim(:).nFrames]./[M.params.stim(:).nUniqueRep];
TR = [M.params.stim(:).framePeriod];
nScans = length(M.params.stim);

x = (0:nFramesInd(1)-1)'.*TR(1);
for n=2:nScans,
    timend = x(end);
    % For scans 2:end, we need to start the indexing at 1, not 0, otherwise
    % the first time point of scan 2 and the last time point of scan 1 will
    % have the same x value. We do not want this.
    % x = [x; ((0:nFramesInd(n)-1)'.*TR(n))+timend];
    x = [x; ((1:nFramesInd(n))'.*TR(n))+timend];

end;
% this does not give the correct separator locations between scans
% M.sepx = cumsum(nFramesInd-1) .* TR;
% rather we want this:
M.sepx = (cumsum(nFramesInd)-1) .* TR;
M.x    = x;

return;

%% if we got here, we haven't changed the voxel:
% we can adjust the rfParams, *if* the option is selected...

% test whether the 'adjust PRF' option is selected
if isequal( get(M.ui.movePRF, 'Checked'), 'on' )
	% either manually move the pRF according to slider values...
	mrvSliderSet(M.ui.moveX, 'Visible', 'on');
    mrvSliderSet(M.ui.moveY, 'Visible', 'on');
	mrvSliderSet(M.ui.moveSigma, 'Visible', 'on');
    rfParams(1) = get(M.ui.moveX.sliderHandle, 'Value');
    rfParams(2) = get(M.ui.moveY.sliderHandle, 'Value');
	rfParams(3) = get(M.ui.moveSigma.sliderHandle, 'Value');
	rfParams(5) = rfParams(3);   % circular case

else
	% ... or use stored values  (and keep the sliders hidden)
    mrvSliderSet(M.ui.moveX, 'Visible', 'off');
    mrvSliderSet(M.ui.moveY, 'Visible', 'off');
    mrvSliderSet(M.ui.moveSigma, 'Visible', 'off');

end

return

author: vistalab

repo: vistasoft