text:
function [sorted_vertices, ...
          h_fes, h_bnd, h_fill, h_vert, h_int, h_max, g_labels] = ...
    plot_feasible(A, b, c, lower_b, upper_b, varargin)
% Plot the feasible region of a linear program
%
% file:      	plot_feasible.m, (c) Matthew Roughan, 2009
% created: 	Mon Jun 22 2009
% author:  	Matthew Roughan
% email:   	matthew.roughan@adelaide.edu.au
%
% Plot the feasible region of the 2D linear program
%           maximize            f = c'*x
%           subject to       A x <= b
% on the region bounded by lower_b and upper_b. Note that the inputs
%     A  is a mx2 matrix, where there are m constraints
%     b  is a mx1 column vector, where there are m constraints
%     c  is a 2x1 column vector
%     lower_b, upper_b are 2x1 column vectors
%
% If the problem includes non-negativity bounds, i.e.,
%               x >= 0
% which is quite common, then these need to be explicitly included into
% the constraints:  A x <= b, not implicitly through the bounds, which are mainly there to
% allow us to make a nice plot.
%
% The feasible region is filled with lines at a setable angle and separation, much as fill
% would do with a solid colour, though options allow you to fill with a color, or both.
%
% This code is primarily (for me) useful for teaching about linear programming, and only
% works for 2D problems (for the obvious reason that problems with more than two variables
% are hard to plot).
%
% Note that if you want to cross-hatch just call this routine twice, with a change in angle
% of 90 degrees.
%
% There are a bunch of optional inputs in pairs, as for the plot command
%    'linecolor'         color of feasible region boundary lines (default black)
%    'backgroundcolor'   background color of feasible region (default no background)
%    'linewidth'         width of lines (default 1)
%    'linestyle'         style of the boundary lines
%    'filllinestyle'     style of the fill lines
%    'linesep'           separation of cross-hatch lines (default 1)
%                        if this is set to -1, don't put any fill lines in
%    'lineangle'         angle (in degrees) of cross-hatch lines
%                           the default is to put them along lines of equal values of the
%                           objective function
%    'hold'              hold_n=1 means keep previous plot
%                        hold_n=0 (default) means clear previous figure
%                        use this to plot multiple feasible regions on same plot, or to do
%                        crosshatching.
%
%    'plot_vertices'     if defined indicates that we should plot the vertices of the feasible region
%                        the value will determine the symbol to use to plot the vertices,
%                        e.g., 'rx', or 'o+'
%    'label_vertices'    if not 0, then label the vertices with their position (label_vertices=1)
%                        or value (label_vertices=1) or both (label_vertices=3)
%    'label_vertices_size'  size of vertex labels, default = 12
%    'label_vertices_color' size of vertex labels, default = 'k'
%    'label_vertices_prec'  precision of vertex labels (default = 2 decimal places)
%
%    'plot_intersections' if defined indicates that we should plot the intersection points of boundaries
%                        the value will determine the symbol to use to plot the vertices,
%                        e.g., 'rx', or 'o+'
%                        Note that vertices are also intersection points, and will be double
%                        plotted if used with the plot_vertices option
%
%    'plot_max'          if not 0, then plot the location of the maximum, using the defined symbol
%
%    'extend_boundaries' if defined, plot each of the constraint lines past the feasible region
%                        using the linewidth and color defined above. The style of the line can
%                        be given as the value, e.g., ':'
%
% OUTPUTS:
%     sorted_vertices = a 2x(n+1) vector of the vertices of the feasible region
%                        note the last one is a repeat of the first to aid in plotting the region
%
%
%     h_fes    = handles to boundary lines of the linear program
%     h_bnd    = handles to the boundary lines
%     h_fill   = handles to the fill lines
%     h_vert   = handles to the vertices (if plotted)
%     h_int    = handles to the intersection points (if plotted)
%     h_max    = a handle to the plotted maximum point
%     g_labels = handles to vertex labels if used
%
% version 0.1, July 22nd 2009
%
% TODO
%     speckled fill
%     testing for boundedness and autosizing the figure
%     plotting bounds with cross-hatch on feasible side
%     at present the code uses the optimization toolkit (through the linprog) function
%       which we use to solve the linear program, but as this is a 2D problem we should be able to
%       remove this dependence with a little work.
%     currently the intersections with upper and lower bounds for plotting are also included
%       in intersection set, and plotted
%     automate ability to put text description of the problem in the window
%

% check sizes of inputs
sA = size(A);
sb = size(b);
if (sA(1) ~= sb(1))
  error('incorrect input sizes\n');
end
if (sA(2) ~= 2)
  error('can only plot for 2 variables\n');
end
n = sA(1); % number of constraints
if (size(upper_b,1) ~=2 | size(upper_b,2) ~= 1)
  error('upper_b is the wrong size\n');
end
if (size(lower_b,1) ~=2 | size(lower_b,2) ~= 1)
  error('lower_b is the wrong size\n');
end

% set default output values
h_fes    = [];
h_bnd    = [];
h_fill   = [];
h_vert   = [];
h_int    = [];
h_max    = [];

% parse the variable input arguments
linecolor = 'k';
backgroundcolor = 0;
linewidth = 1;
linestyle = '-';
filllinestyle = '-';
linesep = 1;
% default line angle is so that lines will be iso-objective
if (abs(c(1)) < 1.0e-12)
  lineangle = 90; % degrees
else
  lineangle = atand(c(2)/c(1));
end
hold_n = 0;
plot_vertices = 0;
label_vertices = 0;
label_vertices_size = 12;
label_vertices_color = 'k';
label_vertices_prec = 2;
plot_intersections = 0;
extend_boundaries = 0;
plot_max = 0;
if (length(varargin) > 0)
  % process variable arguments
  for k = 1:2:length(varargin)
    if (ischar(varargin{k}))
      argument = char(varargin{k});
      value = varargin{k+1};
      switch argument
       case {'linecolor','lc'}
	linecolor = value;
       case {'backgroundcolor','bgc'}
	backgroundcolor = value;
       case {'linewidth','lw'}
	linewidth = value;
       case {'linestyle','ls'}
	 linestyle = value;
       case {'filllinestyle','fls'}
	 filllinestyle = value;
       case {'linesep','ls'}
	linesep = value;
       case {'lineangle','la'}
	lineangle = value;
       case {'hold'}
	hold_n = value;
       case {'plot_vertices', 'pv'}
	plot_vertices = value;
       case {'label_vertices', 'lv'}
	label_vertices = value;
       case {'label_vertices_size', 'lvs'}
	label_vertices_size = value;
       case {'label_vertices_color', 'lvc'}
	label_vertices_color = value;
       case {'label_vertices_prec', 'lvp'}
	label_vertices_prec = value;
       case {'plot_intersections', 'pi'}
	plot_intersections = value;
       case {'extend_boundaries', 'eb'}
	extend_boundaries = value;
       case {'plot_max', 'pm'}
	 plot_max = value;
       otherwise
	error('incorrect input parameters');
      end
    end
  end
end

% set up the plot region
if (hold_n==0)
  hold off
  plot(lower_b(1), lower_b(2));
end
diff_x = (upper_b(1)-lower_b(1))/20;
diff_y = (upper_b(2)-lower_b(2))/20;
set(gca, 'xlim', [lower_b(1)-diff_x upper_b(1)+diff_x]);
set(gca, 'ylim', [lower_b(2)-diff_y upper_b(2)+diff_y]);
hold on
epsilon = 1.0e-14;

% plot the boundary line lines of the linear program across the plot region
if ~(extend_boundaries == 0)
  for i=1:n
    % find end-points of the lines along the boundaries
    constraint = A(i,:);
    k = b(i);
    if (abs(constraint(1)) < epsilon)
      % fprintf('(abs(constraint(1) < epsilon): i = %d\n', i);
      if (lower_b(1)<=k/constraint(2) & k/constraint(2)<=upper_b(1))
	x = [lower_b(1) upper_b(1)];
	y = [1 1] * k/constraint(2);
      else
	x = [];
	y = [];
      end
    elseif (abs(constraint(2)) < epsilon)
      % fprintf('(abs(constraint(2) < epsilon): i = %d\n', i);
      if (lower_b(2)<=k/constraint(1) & k/constraint(1)<=upper_b(2))
	x = [1 1] * k/constraint(1);
	y = [lower_b(2) upper_b(2)];
      else
	x = [];
	y = [];
      end
    else
      % compute all four intersections, and decide which is feasible
      % horizontal intercept points
      y_1 = [lower_b(2) upper_b(2)];
      x_1 = (k - constraint(2)*y_1) / constraint(1);
      i_1 = find(lower_b(1)<x_1 & x_1<upper_b(1));
      % vertical intercept points
      x_2 = [lower_b(1) upper_b(1)];
      y_2 = (k - constraint(1)*x_2) / constraint(2);
      i_2 = find(lower_b(2)<y_2 & y_2<upper_b(2));

      x = [x_1(i_1) x_2(i_2)];
      y = [y_1(i_1) y_2(i_2)];
    end
    % plot the line
    h_fes(i) = plot(x, y, 'color', linecolor, 'linestyle', extend_boundaries);
  end
end

% find the vertices of the lines, plus the borders
Aex = [A; -eye(2); eye(2)];
bex = [b; -lower_b; upper_b];
m = n + 4;
intersections = [];
vertices = [];
for i=1:m-1
  for j=i+1:m
    M = [Aex(i,:); Aex(j,:)];
    bm = bex([i,j]);
    x = M \ bm;
    intersections = [intersections, x];
    % test for feasibility
    if (A*x <= b)
      vertices = [vertices, x];
    end
  end
end

% find and fill the polygon defined by the vertices
k = convhull(vertices(1,:), vertices(2,:));
sorted_vertices(1,:) = vertices(1,k)';
sorted_vertices(2,:) = vertices(2,k)';
if (backgroundcolor == 0)
  h_bnd = plot(sorted_vertices(1,:), sorted_vertices(2,:),  ...
	       'color', linecolor, 'linewidth', linewidth, 'linestyle', linestyle);
else
  h_bnd = fill(sorted_vertices(1,:), sorted_vertices(2,:), backgroundcolor, ...
	       'edgecolor', linecolor, 'linewidth', linewidth, 'linestyle', linestyle);
end

% plot the intersection points and vertices
if ~(plot_intersections == 0)
  h_int = plot(intersections(1,:), intersections(2, :), plot_intersections);
end

if ~(plot_vertices == 0)
  h_vert = plot(sorted_vertices(1,:), sorted_vertices(2, :), plot_vertices);

  if (label_vertices == 1)
    % write the position of the vertex
    format_str = sprintf('(%%.%.0ff,%%.%.0ff)', label_vertices_prec, label_vertices_prec);
    for i=1:size(sorted_vertices, 2)
      temp = 10^label_vertices_prec;
      x = round(temp*sorted_vertices(1,i))/temp;
      y = round(temp*sorted_vertices(2,i))/temp;
      g_labels(i) = text(x+diff_x/2, y, sprintf(format_str, x, y), ...
			 'fontsize', label_vertices_size, 'color', label_vertices_color);
    end
  elseif (label_vertices == 2)
    % write the value of the vertex
    format_str = sprintf('f = %%.%.0ff', label_vertices_prec);
    for i=1:size(sorted_vertices, 2)
      temp = 10^label_vertices_prec;
      x = sorted_vertices(1,i);
      y = sorted_vertices(2,i);
      f = round(temp* (c' * [x; y]))/temp;
      g_labels(i) = text(x+diff_x/2, y, sprintf(format_str, f), ...
			 'fontsize', label_vertices_size, 'color', label_vertices_color);
    end
  elseif (label_vertices == 3)
    % write the position and value
    % write the position of the vertex
    format_str = sprintf('(%%.%.0ff,%%.%.0ff), f=%%.%.0ff', ...
			 label_vertices_prec, label_vertices_prec, label_vertices_prec);
    for i=1:size(sorted_vertices, 2)
      temp = 10^label_vertices_prec;
      x = sorted_vertices(1,i);
      y = sorted_vertices(2,i);
      f = round(temp* (c' * [x; y]))/temp;
      x = round(temp*x)/temp;
      y = round(temp*y)/temp;
      g_labels(i) = text(x+diff_x/2, y, sprintf(format_str, x, y, f), ...
			 'fontsize', label_vertices_size, 'color', label_vertices_color);
    end
  end

end

% plot the max if needed
if ~(plot_max == 0)
  x_max = linprog(-c, A, b, [], [], lower_b, upper_b);
  h_max = plot(x_max(1), x_max(2), plot_max);
end

% do the fill lines
if (linesep > 0)
  %   find min and max via linprog to see start and end points of lines
  %      add in implicit boundary constraints when do it
  c_line = [cosd(lineangle), sind(lineangle)];
  x_min = linprog(c_line, A, b, [], [], lower_b, upper_b);
  % plot(x_min(1), x_min(2), 'b*');

  x_max = linprog(-c_line, A, b, [], [], lower_b, upper_b);
  % plot(x_max(1), x_max(2), 'b*');

  distance = sqrt( sum((x_max-x_min).^2)  );

  %   draw lines clipping them at the boundaries
  for i=0:linesep:distance
    counter = ceil(i/linesep) + 1;

    % sigma = x_min + i*(x_max-x_min)/distance;
    sigma = x_min + i*c_line';
    % plot(sigma(1), sigma(2), '+', 'color', linecolor);

    % for each line, consider its intersepts with each possible
    % boundary, and check whether they are feasible
    edge_vertices = [];
    for j=1:m
      M = [Aex(j,:); c_line];
      d = [bex(j); sum(c_line.*sigma')];
      x = M \ d;
      % test for feasibility
      if (A*x <= b)
	edge_vertices = [edge_vertices, x];
      end
    end
    edge_vertices = round(100*edge_vertices)'/100;
    edge_vertices = unique(edge_vertices, 'rows')';
    if (size(edge_vertices,2) >= 2)
      h_fill(counter) = plot(edge_vertices(1,:), edge_vertices(2,:), ...
			     'color', linecolor, 'linestyle', filllinestyle, 'linewidth', linewidth);
    end
  end
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% split1.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function z = split1(x1,x2,f1,f2)
% Input: points x1 and x2, x1 < x2, and corresponding function values f1
%        and f2
% splits the interval [x1,x2] according to the golden section rule
% the part containing the better point gets the larger fraction of the
% interval

function z = split1(x1,x2,f1,f2)
if f1 <= f2
  z = x1 + 0.5*(-1 + sqrt(5))*(x2 - x1);
else
  z = x1 + 0.5*(3 - sqrt(5))*(x2 - x1);
end

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
%  Imbed a file menu to any figure. If file menu exist, it will append
%  to the existing file menu. This file menu includes: Copy to clipboard,
%  print, save, close etc.
%
%  Usage: rri_file_menu(fig);
%
%         rri_file_menu(fig,0) means no 'Close' menu.
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
%--------------------------------------------------------------------

function rri_file_menu(action, varargin)

   if isnumeric(action)
      fig = action;
      action = 'init';
   end

   %  clear the message line,
   %
   h = findobj(gcf,'Tag','MessageLine');
   set(h,'String','');

   if ~strcmp(action, 'init')
      set(gcbf, 'InvertHardcopy','off');
%      set(gcbf, 'PaperPositionMode','auto');
   end

   switch action
      case {'init'}
         if nargin > 1
            init(fig, 1);		% no 'close' menu
         else
            init(fig, 0);
         end
      case {'print_fig'}
         printdlg(gcbf);
      case {'copy_fig'}
         copy_fig;
      case {'export_fig'}
         export_fig;
   end

   return					% rri_file_menu

%------------------------------------------------
%
%  Create (or append) File menu
%
function init(fig, no_close)

   %  search for file menu
   %
   h_file = [];
   menuitems = findobj(fig, 'type', 'uimenu');

   for i=1:length(menuitems)
      filelabel = get(menuitems(i),'label');

      if strcmpi(strrep(filelabel, '&', ''), 'file')
         h_file = menuitems(i);
         break;
      end
   end

   set(fig, 'menubar', 'none');

   if isempty(h_file)
      if isempty(menuitems)
         h_file = uimenu('parent', fig, 'label', 'File');
      else
         h_file = uimenu('parent', fig, 'label', 'Copy Figure');
      end

      h1 = uimenu('parent', h_file, ...
         'callback','rri_file_menu(''copy_fig'');', ...
         'label','Copy to Clipboard');
   else
      h1 = uimenu('parent', h_file, ...
         'callback','rri_file_menu(''copy_fig'');', ...
         'separator','on', ...
         'label','Copy to Clipboard');
   end

   h2 = uimenu(h_file, ...
      'callback','pagesetupdlg(gcbf);', ...
      'label','Page Setup...');

   h2 = uimenu(h_file, ...
      'callback','printpreview(gcbf);', ...
      'label','Print Preview...');

   h2 = uimenu('parent', h_file, ...
      'callback','printdlg(gcbf);', ...
      'label','Print Figure ...');

   h2 = uimenu('parent', h_file, ...
      'callback','rri_file_menu(''export_fig'');', ...
      'label','Save Figure ...');

   arch = computer;
   if ~strcmpi(arch(1:2),'PC')
      set(h1, 'enable', 'off');
   end

   if ~no_close
      h1 = uimenu('parent', h_file, ...
         'callback','close(gcbf);', ...
         'separator','on', ...
         'label','Close');
   end

   return;					% init

%------------------------------------------------
%
%  Copy to clipboard
%
function copy_fig

   arch = computer;
   if(~strcmpi(arch(1:2),'PC'))
      error('copy to clipboard can only be used under MS Windows');
      return;
   end

   print -noui -dbitmap;

   return					% copy_fig

%------------------------------------------------
%
%  Save as an image file
%
function export_fig

   curr = pwd;
   if isempty(curr)
      curr = filesep;
   end

   [selected_file, selected_path] = rri_select_file(curr,'Save As');

   if isempty(selected_file) | isempty(selected_path)
      return;
   end

   filename = [selected_path selected_file];

   if(exist(filename,'file')==2)		% file exist

      dlg_title = 'Confirm File Overwrite';
      msg = ['File ',filename,' exist. Are you sure you want to overwrite it?'];
      response = questdlg(msg,dlg_title,'Yes','No','Yes');

      if(strcmp(response,'No'))
         return;
      end

   end

   old_pointer = get(gcbf,'pointer');
   set(gcbf,'pointer','watch');

   try
      saveas(gcbf,filename);
   catch
      msg = 'ERROR: Cannot save file';
      set(findobj(gcf,'Tag','MessageLine'),'String',msg);
   end

   set(gcbf,'pointer',old_pointer);

   return;					% export_fig

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [inside] = bounding_mesh(pos, pnt, tri);

% BOUNDING_MESH determines if a point is inside/outside a triangle mesh
% whereby the bounding triangle mesh should be closed.
%
% [inside] = bounding_mesh(pos, pnt, tri)
%
% where
%   pos     position of point of interest (can be 1x3 or Nx3)
%   pnt     bounding mesh vertices
%   tri     bounding mesh triangles
%
% See also SOLID_ANGLE

% Copyright (C) 2003, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: bounding_mesh.m 2885 2011-02-16 09:41:58Z roboos $

global fb;
if isempty(fb)
  fb = 0;
end

npos = size(pos, 1);
npnt = size(pnt, 1);
ntri = size(tri, 1);

% determine a cube that encompases the boundary triangulation
bound_min = min(pnt);
bound_max = max(pnt);

% determine a sphere that is completely inside the boundary triangulation
bound_org = mean(pnt);
bound_rad = sqrt(min(sum((pnt - repmat(bound_org, size(pnt,1), 1)).^2, 2)));

inside = zeros(npos, 1);
for i=1:npos
  if fb
    fprintf('%6.2f%%', 100*i/npos);
  end
  if any(pos(i,:)<bound_min) || any(pos(i,:)>bound_max)
    % the point is outside the bounding cube
    inside(i) = 0;
    if fb, fprintf(' outside the bounding cube\n'); end
  elseif sqrt(sum((pos(i,:)-bound_org).^2, 2))<bound_rad
    % the point is inside the interior sphere
    inside(i) = 1;
    if fb, fprintf(' inside the interior sphere\n'); end
  else
    % the point is inside the bounding cube but outside the interior sphere
    % compute the total solid angle of the surface, which is zero for a point outside
    % the triangulation and 4*pi or -4*pi for a point inside (depending on the triangle
    % orientation)
    tmp = pnt - repmat(pos(i,:), npnt, 1);
    solang = solid_angle(tmp, tri);
    if any(isnan(solang))
      inside(i) = nan;
    elseif (abs(sum(solang))-2*pi)<0
      % total solid angle is (approximately) zero
      inside(i) = 0;
    elseif (abs(sum(solang))-2*pi)>0
      % total solid angle is (approximately) plus or minus 4*pi
      inside(i) = 1;
    end
    if fb, fprintf(' solid angle\n'); end
  end
end

author: PatternRecognition

repo: OpenBMI


text:
clear all; close all; clc;

DATADIR = 'WHERE\IS\DATA';
SAVEDIR = 'WHERE\TO\SAVE';
%% SSVEP
SSVEPDATA = 'EEG_SSVEP.mat';
STRUCTINFO = {'EEG_SSVEP_train', 'EEG_SSVEP_test'};
SESSIONS = {'session1', 'session2'};
TOTAL_SUBJECTS = 54;
%% INITIALIZATION
FS = 100;
FREQBAND = [.5 40];
SEGTIME = [0 4000];

params = {
    'Channels', {'Oz'};
    'Class', {'up'; 'left'; 'right'; 'down'};
    'FFTPlot', 'on';
    };
%%
dataLength = (TOTAL_SUBJECTS * 2 * length(STRUCTINFO));
conSMT = cell(1, dataLength);
idx = 1;
%%
fprintf('Visualization tutorial\n')
for sessNum = 1:length(SESSIONS)
    session = SESSIONS{sessNum};
    for subNum = 1:TOTAL_SUBJECTS
        subject = sprintf('s%d',subNum);
        data = importdata(fullfile(DATADIR,session,subject,SSVEPDATA));
        for info = 1:length(STRUCTINFO)
            CNT = data.(STRUCTINFO{info});
            CNT = rmfield(CNT, 'smt');
            CNT = prep_resample(CNT, FS, {'Nr', 0});
            CNT = prep_filter(CNT, {'frequency', FREQBAND});
            SMT = prep_segmentation(CNT, {'interval', SEGTIME});
            SMT = prep_average(SMT);

            conSMT{idx} = SMT;

            fprintf('%06.2f%%...', idx/dataLength*100);
            if mod(idx, 20) == 0
                fprintf('\n');
            end
            idx = idx + 1;
        end
    end
end
fprintf('\n');

avSMT = grandAverage_prototye(conSMT);

vis_plotController(avSMT, [], params);
mkdir(SAVEDIR);
saveas(gcf, fullfile(SAVEDIR, 'SSVEP_GA'), 'fig');


text:
function [ x, seed ] = f_alpha ( n, q_d, alpha, seed )

%*****************************************************************************80
%
%% F_ALPHA generates a 1/F^ALPHA noise sequence.
%
%  Discussion:
%
%    Thanks to Miro Stoyanov for pointing out that the second half of
%    the data returned by the inverse Fourier transform should be
%    discarded, 24 August 2010.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2010
%
%  Author:
%
%    Original C version by Todd Walter.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jeremy Kasdin,
%    Discrete Simulation of Colored Noise and Stochastic Processes
%    and 1/f^a Power Law Noise Generation,
%    Proceedings of the IEEE,
%    Volume 83, Number 5, 1995, pages 802-827.
%
%  Parameters:
%
%    Input, integer N, the number of samples to generate.
%
%    Input, real Q_D, the variance of the noise.
%
%    Input, real ALPHA, the exponent for the noise.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X(N), a sequence sampled with the given power law.
%
%    Output, integer SEED, a seed for the random number generator.
%

%  Set the deviation of the noise.
%
  q_d = sqrt ( q_d );
%
%  Generate the coefficients Hk.
%
  hfa = zeros ( 2 * n, 1 );
  hfa(1) = 1.0;
  for i = 2 : n
    hfa(i) = hfa(i-1) * ( 0.5 * alpha + ( i - 2 ) ) / ( i - 1 );
  end
  hfa(n+1:2*n) = 0.0;
%
%  Fill Wk with white noise.
%
  wfa = zeros ( 2 * n, 1 );
  for i = 1 : n
    [ wfa(i), seed ] = r8_normal_01 ( seed );
  end
  wfa(1:n) = wfa(1:n) * q_d;
  wfa(n+1:2*n) = 0.0;
%
%  Perform the discrete Fourier transforms of Hk and Wk.
%
  [ h_azero, h_a, h_b ] = r8vec_sftf ( 2 * n, hfa );

  [ w_azero, w_a, w_b ] = r8vec_sftf ( 2 * n, wfa );
%
%  Multiply the two complex vectors.
%
  w_azero = w_azero * h_azero;

  for i = 1 : n
    wr = w_a(i);
    wi = w_b(i);
    w_a(i) = wr * h_a(i) - wi * h_b(i);
    w_b(i) = wi * h_a(i) + wr * h_b(i);
  end
%
%  This scaling is introduced only to match the behavior
%  of the Numerical Recipes code...
%
  w_azero = w_azero * 2 * n;

  w_a(1:n-1) = w_a(1:n-1) * n;
  w_b(1:n-1) = w_b(1:n-1) * n;

  w_a(n) = w_a(n) * 2 * n;
  w_b(n) = w_b(n) * 2 * n;
%
%  Take the inverse Fourier transform of the result.
%
  x = r8vec_sftb ( 2 * n, w_azero, w_a, w_b );
%
%  Discard the second half of the inverse Fourier transform.
%
  x = x(1:n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%{
    Copyright (c) 2014, Philipp Krähenbühl
    All rights reserved.

    Redistribution and use in source and binary forms, with or without
    modification, are permitted provided that the following conditions are met:
        * Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.
        * Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.
        * Neither the name of the Stanford University nor the
        names of its contributors may be used to endorse or promote products
        derived from this software without specific prior written permission.

    THIS SOFTWARE IS PROVIDED BY Philipp Krähenbühl ''AS IS'' AND ANY
    EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
    WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
    DISCLAIMED. IN NO EVENT SHALL Philipp Krähenbühl BE LIABLE FOR ANY
    DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
    (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
	 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
	 ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
	 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%}
% Usage example:
% I = imread('peppers.png');
% os = OverSegmentation( I );
% p = Proposal('max_iou', 0.8,...
%              'unary', 130, 5, 'seedUnary()', 'backgroundUnary({0,15})',...
%              'unary', 0, 5, 'zeroUnary()', 'backgroundUnary({0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15})'
%              );
% props = p.propose( os );

classdef GopProposal
    properties (SetAccess = private, GetAccess = private)
        c_p
    end

    methods
        function obj = GopProposal( varargin )
            obj.c_p = gop_mex( 'newProposal', varargin{:} );
        end
        function r = propose(this,os)
            r = gop_mex( 'Proposal_propose', this.c_p, os.c_s );
        end
        function sobj = saveobj(this)
            error( 'You cannot load/save a Proposal object!' );
        end
        function loadobj(this, sobj)
            error( 'You cannot load/save a Proposal object!' );
        end
        function delete(this)
            gop_mex( 'freeProposal', this.c_p );
        end
    end
end

author: Cloud-CV

repo: object-proposals


text:
function [failure_histograms] = analyze_failures(experiment, trackers, sequences)
% analyze_failures Perform failure frequency analysis
%
% This function performs failure frequency analysis for a set of trackers
% on a set of sequences and experiments.
%
% Input:
% - experiment (cell): A cell array of valid experiment structures.
% - trackers (cell): A cell array of valid tracker descriptor structures.
% - sequences (cell): A cell array of valid sequence descriptor structures.
%
% Output:
% - failure_histograms (cell): A cell array (one element for each experiment) of cell arrays (one for each sequence) of double matrices that contain per-frame failure frequencies for all trackers.
%

    repeat = experiment.parameters.repetitions;

    print_text('Failure analysis for experiment %s ...', experiment.name);

    experiment_sequences = convert_sequences(sequences, experiment.converter);

    failure_histograms = cell(1, numel(experiment_sequences));

    for s = 1:length(experiment_sequences)

        print_indent(1);

        failure_histogram = zeros(numel(trackers), experiment_sequences{s}.length);

        print_text('Processing sequence %s ...', experiment_sequences{s}.name);

        for t = 1:length(trackers)

            print_indent(1);

            result_directory = fullfile(trackers{t}.directory, experiment.name, experiment_sequences{s}.name);

            for j = 1:repeat

                result_file = fullfile(result_directory, sprintf('%s_%03d.txt', experiment_sequences{s}.name, j));

                try
                    trajectory = read_trajectory(result_file);
                catch
                    continue;
                end;

                if length(trajectory) < experiment_sequences{s}.length
                    trajectory{end+1:experiment_sequences{s}.length} = 0;
                end;

                [~, failures] = estimate_failures(trajectory, experiment_sequences{s});
                failures = failures(failures <= experiment_sequences{s}.length);
                failure_histogram(t, failures) = failure_histogram(t, failures) + 1;

            end;

            failure_histograms{s} = failure_histogram;

            print_indent(-1);

        end;

        print_indent(-1);

    end;

author: votchallenge

repo: toolkit-legacy


text:
function rxnLink = connectedRxnsByEFM(model, conComp, rxnInLoops, printLevel)
% Find reactions lying in internal cycles that are connected by any EFMs.
% Used for minimizing the number of constraints for the loopless requirement
% when running loopless FVA using localized loopless constraints. This
% function requires EFMtool (`CalculateFluxModes.m`) to work.
%
% USAGE:
%    rxnLink = connectedRxnsByEFM(model, conComp, rxnInLoops)
%
% INPUTS:
%    model:      COBRA model
%    conComp:    reactions connected in the minimal nullspace for internal cycles, computed by `connectedRxnsInNullSpace`
%    rxnInLoops: n-by-2 logical matrix where n = # of rxns in the model
%                rxnInLoops(k, 1) = true => forward direction of the k-th reaction in internal cycles
%                rxnInLoops(k, 2) = true => reverse direction of the k-th reaction in internal cycles
%                Returned by `findMinNull.m`
%
% OPTIONAL INPUT:
%    printLevel: true to show messages when the linkage matrix cannot be computed
%
% OUTPUT:
%    rxnLink:    n-by-n matrix. rxnLink(i, j) = 1 => reactions i and j are connected
%                by an EFM representing an elementary internal cycle.

if nargin < 4
    printLevel = 0;
end
% the path to EFMtool
efmToolpath = which('CalculateFluxModes.m');
if isempty(efmToolpath)
    rxnLink = [];
    if printLevel
        fprintf('EFMtool not in Matlab path. Unable to calculate EFMs.\n')
    end
    return
end
efmToolpath = strsplit(efmToolpath, filesep);
efmToolpath = strjoin(efmToolpath(1: end - 1), filesep);
p = pwd;
cd(efmToolpath)
% EFMtool call options
options = CreateFluxModeOpts('sign-only', true, 'level', 'WARNING');

rxnLink = sparse(size(model.S, 2), size(model.S, 2));
for jC = 1:max(conComp)
    % for each connected component, find the EFM matrix
    try
        S = model.S(:, conComp == jC);
        S = S(any(S, 2), :);
        % revert the stoichiometries for reactions that are in cycles only in the reverse direction
        S(:, rxnInLoops(conComp == jC, 1) & ~rxnInLoops(conComp == jC, 2)) = -S(:, rxnInLoops(conComp == jC, 1) & ~rxnInLoops(conComp == jC, 2));
        rev = all(rxnInLoops(conComp == jC, :), 2);
        efms = CalculateFluxModes(full(S), double(rev), options);
        % calling Java too rapidly may have problems in tests
        pause(1e-4)
        efms = efms.efms;
        rxnJC = find(conComp == jC);
        for j = 1:numel(rxnJC)
            rxnLink(rxnJC(j), rxnJC) = any(efms(:, efms(j, :) ~= 0), 2)';
        end
    catch msg
        if printLevel
            fprintf('Error encountered during calculation of EFMs:\n%s', getReport(msg))
        end
        rxnLink = [];
        return
    end
end
cd(p)

end

author: opencobra

repo: cobratoolbox


text:
% test for ti quicunx

name = 'turbulence';
name = 'lena';
n = 256;
M = load_image(name);
M = rescale( crop(M,n) );
Jmax = log2(n)-1;
Jmin = Jmax-5;

% boundary conditions
options.bound = 'per';
options.bound = 'sym';
% vanishing moments
vm = 6;
options.primal_vm = vm;
options.dual_vm = vm;

% transform
disp('Computing forward transform');
MW = perform_quicunx_wavelet_transform_ti(M,Jmin,options);
disp('Computing backward transform');
M1 = perform_quicunx_wavelet_transform_ti(MW,Jmin,options);
disp(['Error=' num2str(norm(M-M1,'fro')/norm(M,'fro')) ' (should be 0).']);

rep = 'results/quincunx/';
if not(exist(rep))
    mkdir(rep);
end

% display dual wavelets
clf;
for j=7:12
    MW = MW*0;
    MW(end/2,end/2,j) = 1;
    M1 = perform_quicunx_wavelet_transform_ti(MW,Jmin,options);
    imageplot(M1, '', 4,3,j);
    % save
    if j>6
        warning off;
%        imwrite(rescale(M1), [rep 'quincunx-wavelet-' num2string_fixeddigit(j,2) '.png'], 'png');
        warning on;
        clf;
        surf(M1);
        shading interp; colormap jet(256)
        view(-20,50); axis tight; axis off;
        camlight;
        saveas(gcf, [rep 'quincunx-wavelet-vm' num2str(vm) '-' num2string_fixeddigit(j,2) '.png'], 'png');
    end
end

author: gpeyre

repo: matlab-toolboxes


text:
function H = hurst_estimate(sequence,method,isplot,opt)
%
% 'hurst_estimate' estimate the hurst parameter of a given sequence with
%     an appointed method. The algorithms of the methods can be found in
%     Murad's Taqqu, Vadim Teverovsky and Walter Willinger's paper
%     "Estimators for long-range dependence: an empirical study" or other
%     related papers.
% Inputs:
%     sequence: the input sequence for estimate
%     method: the name of a function which used to estimate the hurst
%             parameter of the sequence,e.g.
%               'aggvar': use aggvar function to estimate.
%               'RS': use RS function to estimate.
%               'per': use per function to estimate.
%     isplot: whether display the plot. without a plot if isplot equal to 0
%     opt: a optional parameter for some methods
% Outputs:
%     H: the estimated hurst coeffeient of the input sequence
% Examples:
%     H = hurst_estimate('peng',sequence,1);
%

%  Author: Chu Chen
%  Version 1.0,  03/10/2008
%  chen-chu@163.com
%

if nargin == 2
    isplot = 0;
    H = feval(method,sequence,isplot);
elseif nargin == 3
    H = feval(method,sequence,isplot);
elseif nargin == 4
    H = feval(method,sequence,isplot,opt);
end

author: Sable

repo: mcbench-benchmarks


text:
function value = identity_condition ( n )

%*****************************************************************************80
%
%% IDENTITY_CONDITION returns the L1 condition of the IDENTITY matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real VALUE, the L1 condition.
%
  a_norm = 1.0;
  b_norm = 1.0;
  value = a_norm * b_norm;

  return
end


text:
classdef GraphicalAnalysisIndepVarEnum < matlab.mixin.SetGet
    %GraphicalAnalysisIndepVarEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        Time('Time')
        MissionElapsedTime('Mission Elapsed Time');
        TrueAnomaly('True Anomaly');
        Longitude('Longitude')
        Altitude('Altitude')
    end

    properties
        name(1,:) char
    end

    methods
        function obj = GraphicalAnalysisIndepVarEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function [listBoxStr, m] = getListboxStr()
            m = enumeration('GraphicalAnalysisIndepVarEnum');
            listBoxStr = {m.name};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [J, grad] = cofiCostFunc(params, Y, R, num_users, num_movies, ...
                                  num_features, lambda)
%COFICOSTFUNC Collaborative filtering cost function
%   [J, grad] = COFICOSTFUNC(params, Y, R, num_users, num_movies, ...
%   num_features, lambda) returns the cost and gradient for the
%   collaborative filtering problem.
%

% Unfold the U and W matrices from params
X = reshape(params(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(params(num_movies*num_features+1:end), ...
                num_users, num_features);

% You need to return the following values correctly
J = 0;
X_grad = zeros(size(X));
Theta_grad = zeros(size(Theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost function and gradient for collaborative
%               filtering. Concretely, you should first implement the cost
%               function (without regularization) and make sure it is
%               matches our costs. After that, you should implement the
%               gradient and use the checkCostFunction routine to check
%               that the gradient is correct. Finally, you should implement
%               regularization.
%
% Notes: X - num_movies  x num_features matrix of movie features
%        Theta - num_users  x num_features matrix of user features
%        Y - num_movies x num_users matrix of user ratings of movies
%        R - num_movies x num_users matrix, where R(i, j) = 1 if the
%            i-th movie was rated by the j-th user
%
% You should set the following variables correctly:
%
%        X_grad - num_movies x num_features matrix, containing the
%                 partial derivatives w.r.t. to each element of X
%        Theta_grad - num_users x num_features matrix, containing the
%                     partial derivatives w.r.t. to each element of Theta
%

J_temp = (X * Theta' - Y).^2;
J = sum(sum(J_temp(R == 1)))/2 + lambda/2 .* sum(sum(Theta.^2)) + lambda/2 .* sum(sum(X.^2));

X_grad = ((X * Theta' - Y) .* R) * Theta + lambda.*X;
Theta_grad = ((X * Theta' - Y) .* R)' * X + lambda.*Theta;

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

grad = [X_grad(:); Theta_grad(:)];

end

author: Borye

path: github-repos/MATLAB/Borye-machine-learning-coursera-1/machine-learning-coursera-1-033fdc2e6da393eeb1179a09aafe92362021effb/Week 9 Assignments/Anomaly Detection and Recommender Systems/mlclass-ex8/cofiCostFunc.m


text:
%% Create a TSE sequence and export for execution
%
% The |Sequence| class provides functionality to create magnetic
% resonance sequences (MRI or NMR) from basic building blocks.
%
% This provides an implementation of the open file format for MR sequences
% described here: http://pulseq.github.io/specification.pdf
%
% This example performs the following steps:
%
% # Create slice selective RF pulse for imaging.
% # Create readout gradient and phase encode strategy.
% # Loop through phase encoding and generate sequence blocks.
% # Write the sequence to an open file format suitable for execution on a
% scanner.
%
%   Juergen Hennig <juergen.hennig@uniklinik-freiburg.de>
%   Maxim Zaitsev  <maxim.zaitsev@uniklinik-freiburg.de>

%% Instantiation and gradient limits
% The system gradient limits can be specified in various units _mT/m_,
% _Hz/cm_, or _Hz/m_. However the limits will be stored internally in units
% of _Hz/m_ for amplitude and _Hz/m/s_ for slew. Unspecificied hardware
% parameters will be assigned default values.

dG=250e-6;
system = mr.opts('MaxGrad', 30, 'GradUnit', 'mT/m', ...
    'MaxSlew', 170, 'SlewUnit', 'T/m/s', 'rfRingdownTime', 100e-6, ...
    'rfDeadTime', 100e-6, 'adcDeadTime', 10e-6);

%%
% A new sequence object is created by calling the class constructor.
seq=mr.Sequence(system);

%% Sequence events
% Some sequence parameters are defined using standard MATLAB variables
fov=256e-3;
Nx=128; Ny=128; necho=16; Nslices=1;
rflip=180;
if (numel(rflip)==1), rflip=rflip+zeros([1 necho]); end
sliceThickness=5e-3;
TE=12e-3; TR=2000e-3;
TEeff=60e-3;
k0=round(TEeff/TE);
PEtype='linear';

samplingTime= 6.4e-3;
readoutTime = samplingTime + 2*system.adcDeadTime;
tEx=2.5e-3;
tExwd=tEx+system.rfRingdownTime+system.rfDeadTime;
tRef=2e-3;
tRefwd=tRef+system.rfRingdownTime+system.rfDeadTime;
tSp=0.5*(TE-readoutTime-tRefwd);
tSpex=0.5*(TE-tExwd-tRefwd);
fspR=1.0;
fspS=0.5;

rfex_phase=pi/2; % MZ: we need to maintain these as variables because we will overwrtite phase offsets for multiple slice positions
rfref_phase=0;

%%
%%% Base gradients
%%% Slice selection
% Key concepts in the sequence description are *blocks* and *events*.
% Blocks describe a group of events that are executed simultaneously. This
% hierarchical structure means that one event can be used in multiple
% blocks, a common occurrence in MR sequences, particularly in imaging
% sequences.
%
% First, the slice selective RF pulses (and corresponding slice gradient)
% are generated using the |makeSincPulse| function.
% Gradients are recalculated such that their flattime covers the pulse plus
% the rfdead- and rfringdown- times.
%
flipex=90*pi/180;
[rfex, gz] = mr.makeSincPulse(flipex,system,'Duration',tEx,...
    'SliceThickness',sliceThickness,'apodization',0.5,'timeBwProduct',4,'PhaseOffset',rfex_phase);
GSex = mr.makeTrapezoid('z',system,'amplitude',gz.amplitude,'FlatTime',tExwd,'riseTime',dG);
% plotPulse(rfex,GSex);

flipref=rflip(1)*pi/180;
[rfref, gz] = mr.makeSincPulse(flipref,system,'Duration',tRef,...
    'SliceThickness',sliceThickness,'apodization',0.5,'timeBwProduct',4,'PhaseOffset',rfref_phase,'use','refocusing');
GSref = mr.makeTrapezoid('z',system,'amplitude',GSex.amplitude,'FlatTime',tRefwd,'riseTime',dG);
% plotPulse(rfref,GSref);

AGSex=GSex.area/2;
GSspr = mr.makeTrapezoid('z',system,'area',AGSex*(1+fspS),'duration',tSp,'riseTime',dG);
GSspex = mr.makeTrapezoid('z',system,'area',AGSex*fspS,'duration',tSpex,'riseTime',dG);

%%
%%% Readout gradient
% To define the remaining encoding gradients we need to calculate the
% $k$-space sampling. The Fourier relationship
%
% $$\Delta k = \frac{1}{FOV}$$
%
% Therefore the area of the readout gradient is $n\Delta k$.
deltak=1/fov;
kWidth = Nx*deltak;

GRacq = mr.makeTrapezoid('x',system,'FlatArea',kWidth,'FlatTime',readoutTime,'riseTime',dG);
adc = mr.makeAdc(Nx,'Duration',samplingTime, 'Delay', system.adcDeadTime);%,'Delay',GRacq.riseTime);
GRspr = mr.makeTrapezoid('x',system,'area',GRacq.area*fspR,'duration',tSp,'riseTime',dG);
GRspex = mr.makeTrapezoid('x',system,'area',GRacq.area*(1+fspR),'duration',tSpex,'riseTime',dG);

AGRspr=GRspr.area;%GRacq.area/2*fspR;
AGRpreph = GRacq.area/2+AGRspr;%GRacq.area*(1+fspR)/2;
GRpreph = mr.makeTrapezoid('x',system,'Area',AGRpreph,'duration',tSpex,'riseTime',dG);

%%
%%% Phase encoding
% To move the $k$-space trajectory away from 0 prior to the readout a
% prephasing gradient must be used. Furthermore rephasing of the slice
% select gradient is required.

%[PEorder,Ny] = myTSE_PEorder(Ny,necho,k0,PEtype);
nex=floor(Ny/necho);
pe_steps=(1:(necho*nex))-0.5*necho*nex-1;
if 0==mod(necho,2)
    pe_steps=circshift(pe_steps,[0,-round(nex/2)]); % for odd number of echoes we have to apply a shift to avoid a contrast jump at k=0
end
PEorder=reshape(pe_steps,[nex,necho])';
phaseAreas = PEorder*deltak;

%% split gradients and recombine into blocks
% lets start with slice selection....
GS1times=[0 GSex.riseTime];
GS1amp=[0 GSex.amplitude];
GS1 = mr.makeExtendedTrapezoid('z','times',GS1times,'amplitudes',GS1amp);

GS2times=[0 GSex.flatTime];
GS2amp=[GSex.amplitude GSex.amplitude];
GS2 = mr.makeExtendedTrapezoid('z','times',GS2times,'amplitudes',GS2amp);

GS3times=[0 GSspex.riseTime GSspex.riseTime+GSspex.flatTime GSspex.riseTime+GSspex.flatTime+GSspex.fallTime];
GS3amp=[GSex.amplitude GSspex.amplitude GSspex.amplitude GSref.amplitude];
GS3 = mr.makeExtendedTrapezoid('z','times',GS3times,'amplitudes',GS3amp);

GS4times=[0 GSref.flatTime];
GS4amp=[GSref.amplitude GSref.amplitude];
GS4 = mr.makeExtendedTrapezoid('z','times',GS4times,'amplitudes',GS4amp);

GS5times=[0 GSspr.riseTime GSspr.riseTime+GSspr.flatTime GSspr.riseTime+GSspr.flatTime+GSspr.fallTime];
GS5amp=[GSref.amplitude GSspr.amplitude GSspr.amplitude 0];
GS5 = mr.makeExtendedTrapezoid('z','times',GS5times,'amplitudes',GS5amp);

GS7times=[0 GSspr.riseTime GSspr.riseTime+GSspr.flatTime GSspr.riseTime+GSspr.flatTime+GSspr.fallTime];
GS7amp=[0 GSspr.amplitude GSspr.amplitude GSref.amplitude];
GS7 = mr.makeExtendedTrapezoid('z','times',GS7times,'amplitudes',GS7amp);

% and now the readout gradient....

GR3=GRpreph;%GRspex;

GR5times=[0 GRspr.riseTime GRspr.riseTime+GRspr.flatTime GRspr.riseTime+GRspr.flatTime+GRspr.fallTime];
GR5amp=[0 GRspr.amplitude GRspr.amplitude GRacq.amplitude];
GR5 = mr.makeExtendedTrapezoid('x','times',GR5times,'amplitudes',GR5amp);

GR6times=[0 readoutTime];
GR6amp=[GRacq.amplitude GRacq.amplitude];
GR6 = mr.makeExtendedTrapezoid('x','times',GR6times,'amplitudes',GR6amp);

GR7times=[0 GRspr.riseTime GRspr.riseTime+GRspr.flatTime GRspr.riseTime+GRspr.flatTime+GRspr.fallTime];
GR7amp=[GRacq.amplitude GRspr.amplitude GRspr.amplitude 0];
GR7 = mr.makeExtendedTrapezoid('x','times',GR7times,'amplitudes',GR7amp);

% and filltimes
%tex=GS1.t(end)+GS2.t(end)+GS3.t(end);
%tref=GS4.t(end)+GS5.t(end)+GS7.t(end)+readoutTime;
%tend=GS4.t(end)+GS5.t(end);
tex=mr.calcDuration(GS1)+mr.calcDuration(GS2)+mr.calcDuration(GS3);
tref=mr.calcDuration(GS4)+mr.calcDuration(GS5)+mr.calcDuration(GS7)+readoutTime;
tend=mr.calcDuration(GS4)+mr.calcDuration(GS5);

tETrain=tex+necho*tref+tend;
TRfill=(TR-Nslices*tETrain)/Nslices;
% round to gradient raster
TRfill=system.gradRasterTime * round(TRfill / system.gradRasterTime);
if TRfill<0, TRfill=1e-3;
    disp(strcat('Warning!!! TR too short, adapted to include all slices to : ',num2str(1000*Nslices*(tETrain+TRfill)),' ms'));
else
    disp(strcat('TRfill : ',num2str(1000*TRfill),' ms'));
end
delayTR = mr.makeDelay(TRfill);

%% Define sequence blocks
% Next, the blocks are put together to form the sequence
for kex=0:nex % MZ: we start at 0 to have one dummy
    for s=1:Nslices
        rfex.freqOffset=GSex.amplitude*sliceThickness*(s-1-(Nslices-1)/2);
        rfref.freqOffset=GSref.amplitude*sliceThickness*(s-1-(Nslices-1)/2);
        rfex.phaseOffset=rfex_phase-2*pi*rfex.freqOffset*mr.calcRfCenter(rfex); % align the phase for off-center slices
        rfref.phaseOffset=rfref_phase-2*pi*rfref.freqOffset*mr.calcRfCenter(rfref); % dito

        seq.addBlock(GS1);
        seq.addBlock(GS2,rfex);
        seq.addBlock(GS3,GR3);
        %GS4.first=GS4f;
        %GS4.first=GS3.last;
        for kech=1:necho,
            if (kex>0)
                phaseArea=phaseAreas(kech,kex);
            else
                phaseArea=0;
            end
            GPpre = mr.makeTrapezoid('y',system,'Area',phaseArea,'Duration',tSp,'riseTime',dG);
            GPrew = mr.makeTrapezoid('y',system,'Area',-phaseArea,'Duration',tSp,'riseTime',dG);
            seq.addBlock(GS4,rfref);
            seq.addBlock(GS5,GR5,GPpre);
            if (kex>0)
                seq.addBlock(GR6,adc);
            else
                seq.addBlock(GR6);
            end
            seq.addBlock(GS7,GR7,GPrew);
            %GS4.first=GS7.last;
        end
        seq.addBlock(GS4);
        seq.addBlock(GS5);
        seq.addBlock(delayTR);
    end
end

%% check whether the timing of the sequence is correct
[ok, error_report]=seq.checkTiming;

if (ok)
    fprintf('Timing check passed successfully\n');
else
    fprintf('Timing check failed! Error listing follows:\n');
    fprintf([error_report{:}]);
    fprintf('\n');
end

%% k-space trajectory calculation
[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing] = seq.calculateKspacePP();

% plot k-spaces
figure; plot(t_ktraj,ktraj'); title('k-space components as functions of time'); % plot the entire k-space trajectory
figure; plot(ktraj(1,:),ktraj(2,:),'b',...
             ktraj_adc(1,:),ktraj_adc(2,:),'r.'); % a 2D plot
axis('equal'); % enforce aspect ratio for the correct trajectory display
title('2D k-space');

%% Write to file

% The sequence is written to file in compressed form according to the file
% format specification using the |write| method.
seq.write('tse.seq')

%%
% Display the first few lines of the output file
% s=fileread('myTSE.seq');
% disp(s(1:300))
seq.plot();

% seq.install('siemens');

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
function test_suite = test_meeg_baseline_correct
% tests for cosmo_meeg_baseline_correct
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_meeg_baseline_correct_ft_comparison_methods()
    methods={'relative','absolute','relchange'};
    references={'manual','ft','ds'};

    combis=cosmo_cartprod({methods,references});
    ncombis=size(combis,1);

    for k=1:ncombis
        combi=combis(k,:);
        method=combi{1};
        reference=combi{2};

        % random interval
        start=rand(1)*-.25;
        dur=.05+rand(1)*.2;
        interval=start+[0 dur];

        helper_test_meeg_baseline_correct_comparison(interval,...
                                                        method,reference);
    end

function helper_test_meeg_baseline_correct_comparison(interval,method,...
                                                    reference)

    switch reference
        case 'ft'
            if cosmo_skip_test_if_no_external('fieldtrip')
                return;
            end

            if cosmo_wtf('is_octave')
                cosmo_notify_test_skipped(['ft_freqbaseline is not '...
                                                'compatible with Octave']);
            end

        case 'ds'

        case 'manual'

        otherwise
            assert(false,'this should not happen');
    end

    ds=cosmo_synthetic_dataset('type','timefreq','size','big','seed',0);
    ds.samples=randn(size(ds.samples));

    chan_to_select=cosmo_randperm(max(ds.fa.chan),ceil(rand()*3+1));
    freq_to_select=cosmo_randperm(max(ds.fa.freq),ceil(rand()*3+1));

    m=cosmo_match(ds.fa.chan,chan_to_select) & ...
            cosmo_match(ds.fa.freq,freq_to_select);
    ds=cosmo_slice(ds,m,2);
    ds=cosmo_dim_prune(ds);

    y=cosmo_meeg_baseline_correct(ds,interval,method);

    % (unsupported in octave)
    switch reference
        case 'ft'
            ft=cosmo_map2meeg(ds);

            opt=struct();
            opt.baseline=interval;
            opt.baselinetype=method;

            ft_bl=ft_freqbaseline(opt,ft);
            x=cosmo_meeg_dataset(ft_bl);

        case 'ds'
            msk=cosmo_dim_match(ds,'time',...
                        @(t) t>=min(interval) & t<=max(interval));
            ds_ref=cosmo_slice(ds,msk,2);
            x=cosmo_meeg_baseline_correct(ds,ds_ref,method);

        case 'manual'
            msk=cosmo_dim_match(ds,'time',...
                        @(t) t>=min(interval) & t<=max(interval));
            ds_ref=cosmo_slice(ds,msk,2);

            x=ds;

            for chan=1:max(ds.fa.chan)
                for freq=1:max(ds.fa.freq)
                    msk=ds.fa.chan==chan & ds.fa.freq==freq;
                    s=ds.samples(:,msk);

                    ref_msk=ds_ref.fa.chan==chan & ds_ref.fa.freq==freq;
                    r=mean(ds_ref.samples(:,ref_msk),2);

                    switch method
                        case 'absolute'
                            v=bsxfun(@minus,s,r);

                        case 'relative'
                            v=bsxfun(@rdivide,s,r);

                        case 'relchange'
                            v=bsxfun(@rdivide,bsxfun(@minus,s,r),r);

                        otherwise
                            error('not supported: %s', method);
                    end

                   x.samples(:,msk)=v;
                end
            end

        otherwise
            assert(false)

    end

    x_unq=cosmo_index_unique({x.fa.time,x.fa.chan});
    y_unq=cosmo_index_unique({y.fa.time,y.fa.chan});

    n=numel(x_unq);
    max_n_to_choose=10;
    n_to_choose=min(max_n_to_choose,n);
    rp=cosmo_randperm(n,n_to_choose);
    assert(numel(x_unq)==numel(y_unq));
    for j=1:n_to_choose;
        idx=rp(j);
        x_sel=cosmo_slice(x,x_unq{idx},2,false);
        y_sel=cosmo_slice(y,y_unq{idx},2,false);

        p=x_sel.samples;
        q=y_sel.samples;

        desc=sprintf('method %s, reference %s',method,reference);
        assertElementsAlmostEqual(p,q,...
                        'relative',1e-6,desc);
        assertEqual(x_sel.fa,y_sel.fa);

    end

    assertEqual(x.a.fdim,y.a.fdim);

function test_meeg_baseline_correct_regression()
    interval=[-.15 -.04];
    expected_samples={[ 2.7634  -2.6366   2.3689 -0.36892 -0.43043;...
                        -0.49007  -1.0933  0.79082   1.2092   0.1044 ],...
                      [ -0.63503   1.3096 -0.49296  0.49296  0.51511;...
                        -2.5308  -3.5553 -0.35529  0.35529  -1.5211 ],...
                      [ 1.7634  -3.6366   1.3689  -1.3689  -1.4304
                        -1.4901  -2.0933 -0.20918  0.20918  -0.8956 ]};
    [ds,ds_ref]=get_test_dataset(interval);

    methods={'relative','absolute','relchange'};

    ds_feature_msk=ds.fa.chan==3 & ds.fa.freq==2;

    for k=1:numel(methods)
        method=methods{k};

        for j=1:2
            if j==1
                ref=interval;
            else
                ref=ds_ref;
            end

            ds_bl=cosmo_meeg_baseline_correct(ds,ref,method);
            ds_bl_msk=ds_bl.fa.chan==3 & ds_bl.fa.freq==2;
            d=cosmo_slice(ds_bl,ds_bl_msk,2);

            assertElementsAlmostEqual(d.samples,expected_samples{k},...
                                                'absolute',1e-4);
            d_fa=cosmo_slice(ds.fa,ds_feature_msk,2,'struct');
            assertEqual(d.fa,d_fa);
        end
    end

function [ds,ds_ref]=get_test_dataset(interval)
    ds=cosmo_synthetic_dataset('type','timefreq','size','big',...
                                    'senstype','neuromag306_planar',...
                                    'nchunks',1);
    ds.sa=struct();
    ds.sa.rpt=(1:size(ds.samples,1))';
    msk=ds.fa.chan<=4 & ds.fa.freq<=2;
    ds=cosmo_slice(ds,msk,2);
    ds=cosmo_dim_prune(ds);

    matcher=@(x) interval(1) <= x & x <= interval(2);
    ds_ref=cosmo_slice(ds,cosmo_dim_match(ds,'time',matcher),2);
    ds_ref=cosmo_dim_prune(ds_ref);

function test_meeg_baseline_correct_nonmatching_sa
    ds=cosmo_synthetic_dataset('size','big','ntargets',8,...
                            'nchunks',1,'type','timelock');
    nsamples=size(ds.samples,1);

    while true
        rp=cosmo_randperm(nsamples);
        if ~isequal(rp,1:nsamples)
            break;
        end
    end
    ds_ref=cosmo_slice(ds,rp);
    assertExceptionThrown(@()cosmo_meeg_baseline_correct(...
                            ds,ds_ref,'relative'),'');

function test_meeg_baseline_correct_nonmatching_fa
    ds_big=cosmo_synthetic_dataset('size','big','ntargets',8,...
                            'nchunks',1,'type','timefreq');
    ds=cosmo_slice(ds_big,ds_big.fa.chan<=2 & ds_big.fa.freq<=3,2);
    ds_ref=cosmo_slice(ds_big,ds_big.fa.chan<=3 & ds_big.fa.freq<=2,2);

    assertExceptionThrown(@()cosmo_meeg_baseline_correct(...
                                        ds,ds_ref,'relative'),'');

function test_meeg_baseline_correct_illegal_inputs
    bc=@cosmo_meeg_baseline_correct;
    aet=@assertExceptionThrown;
    ds=cosmo_synthetic_dataset('type','timefreq','size','tiny');

    if cosmo_wtf('is_matlab')
        v=cosmo_wtf('version');
        is_prior_to_2012b=str2num(v(1))<=7;

        if is_prior_to_2012b
            id_missing_arg='MATLAB:inputArgUndefined';
        else
            id_missing_arg='MATLAB:minrhs';
        end
    else
        id_missing_arg='Octave:undefined-function';
    end

    aet(@()bc(ds,ds),id_missing_arg);
    aet(@()bc(ds,ds,'foo'),'');

    aet(@()bc(ds,cosmo_slice(ds,1),'relative'),'');
    aet(@()bc(ds,cosmo_slice(ds,1,2),'relative'),'');

    % test slicing
    bc(ds,cosmo_slice(ds,[1 2 3 4 5 6],1),'relative');
    aet(@()bc(ds,cosmo_slice(ds,[1 2 4 6 4 3],1),'relative'),'');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
%kkblend 'Output = (Input 1 * Constant) + (Input 2 (1 - Constant))'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kblend.pane file
%
% Parameters:
% InputFile: i1 'Input 1 ', required: 'first input data object'
% InputFile: i2 'Input 2 ', required: 'second input file object'
% OutputFile: o 'Output', required: 'resulting output data object'
% InputFile: iratio 'Blending file', optional: 'Input object for providing blending ratios'
%
% Example: o = kkblend({i1, i2, iratio}, {'i1','';'i2','';'o','';'iratio',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kblend - Output = (Input 1 * Constant) + (Input 2 (1 - Constant))
%
%  DESCRIPTION
% The "Blend" operator (kblend) forms a blend of two objects by
% adding X*(object1)+(1-X)*(object2).
%
% Instead of specifying a constant for the blending ratio, an object can be
% supplied for the blending values. The value segment of the ratio object
% is used to determine the proportion in which the two objects are blended.
% It must take on values between 0 and 1, otherwise an error is generated.
%
%
%  "Data Type" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_2input
%
%  "Map Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/map_2input
%
%  "Validity Mask" 5
% If there is a validity mask associated with any of the input objects, the
% resulting output object will have a mask.  The output object mask is the
% logical AND of the input object masks.  If only one input object has a mask,
% then the other input object's mask is assumed to be "all valid" (mask
% value = 1).  If a mask is resized for the operation, the mask padding value
% will be 1.
%
%  "Input Objects of Different Sizes" 5
% .cI $DATAMANIP/repos/shared/man/sections/resize_2input
%
%  "Explicit Location and Time Data" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_2input
%
%  "Failure Modes - Two Input Objects"
% .cI $DATAMANIP/repos/shared/man/sections/fail_2input
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% The blend operator has not been updated to be fully polymorphic.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkblend(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkblend(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'o', '__output';'iratio', '__input'};
maxval={0,0,0,1};
minval={0,0,0,1};
istoggle=[0,0,0,1];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','OutputFile','InputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kblend"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function Rob = map2rob(Rob)
% MAP2ROB  Update Rob structure from the Map information.
%   ROB = MAP2ROB(ROB) updates the structure ROB to reflect the information
%   contained in the golbal map Map.
%
%   See also UPDATEFRAME.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

global Map

% normalize quaternion - mean and cov
[Map.x(Rob.frame.r(4:7)), NQ_q] = normvec(Map.x(Rob.frame.r(4:7)),1);
Map.P(Rob.frame.r(4:7),Map.used) = NQ_q*Map.P(Rob.frame.r(4:7),Map.used);
Map.P(Map.used,Rob.frame.r(4:7)) = Map.P(Map.used,Rob.frame.r(4:7))*NQ_q';

% means
Rob.state.x = Map.x(Rob.state.r);
Rob.frame.x = Map.x(Rob.frame.r);

Rob.frame   = updateFrame(Rob.frame);

% covariances
% Rob.state.P = Map.P(Rob.state.r,Rob.state.r);
% Rob.frame.P = Map.P(Rob.frame.r,Rob.frame.r);

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [hdr] = read_biosig_header(filename)

% READ_BIOSIG_HEADER reads header from EEG file using the BIOSIG
% toolbox and returns it in the FCDC framework standard format
%
% Use as
%  [hdr] = read_biosig_header(filename)
%
% The following data formats are supported: EDF, BKR, CNT, BDF, GDF,
% see for full documentation http://biosig.sourceforge.net/
%
% See also READ_BIOSIG_DATA

% Copyright (C) 2004, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_biosig_header.m 945 2010-04-21 17:41:20Z roboos $

% open the file, read the header and close it again
biosig = sopen(filename,'r');
sclose(biosig);

% the BIOSIG header is defined in full detail below
% the FCDC header should at least contain the following fields
%   hdr.Fs           sampling frequency
%   hdr.nChans       number of channels
%   hdr.nSamples     number of samples per trial
%   hdr.nSamplesPre  number of pre-trigger samples in each trial
%   hdr.nTrials      number of trials
%   hdr.label        cell-array with labels of each channel

if length(biosig.SampleRate)>1 & any(diff(biosig.SampleRate))
  error('channels with different sampling rates are not supported');
else
  hdr.Fs          = biosig.SampleRate(1);
end

if length(biosig.SPR)>1 & any(diff(biosig.SPR))
  error('channels with different number of samples are not supported');
else
  hdr.nSamples  = biosig.SPR(1);
end

hdr.nChans      = biosig.NS;
hdr.nTrials     = biosig.NRec;
hdr.nSamplesPre = 0;      % this one is not in the biosig header
hdr.label       = {};     % start with empty labels and fill them below

if isfield(biosig, 'Label')
  hdr.label = biosig.Label;
end

if length(hdr.label)~=hdr.nChans
  % make default channel labels
  hdr.label = {};
  for i=1:hdr.nChans
    hdr.label{i} = sprintf('%d', i);
  end
end

% I prefer to have them as column vector
hdr.label = hdr.label(:);

% also remember the biosig header details
hdr.orig = biosig;
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the BIOSIG header always contains these elements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HDR.TYPE                string          type of data format
% HDR.VERSION             string          (depends on data format)
% HDR.T0                  float[1..6]     [yyyy mm dd hh MM ss.cc] see HELP CLOCK
% HDR.NS                  integer         number of channels
% HDR.SampleRate          integer         sampling frequency in [Hz]
% HDR.NRec                integer         number of records or blocks; 1 for continous data
% HDR.SPR                 integer         samples per record
% HDR.Dur                 float           Duration (in [s]) of minimal block length
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the BIOSIG header optinally contains the following elements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HDR.Filter.LowPass      float           [Hz]
% HDR.Filter.HighPass     float           [Hz]
% HDR.Filter.Notch        int8            0=Off, 1=ON
% HDR.PreFilt             string          filter setting
% HDR.Label               char-array      z.B. '+C3a - C3p  '
% HDR.PhysDim             string          physical dimension e.g. 'uV'
% HDR.PhysMax             float           physical maximum
% HDR.DigMax              integer         digital maximum
% HDR.PhysMin             float           physical minimum
% HDR.DigMin              integer         digital minimum
% HDR.FLAG.TRIGGERED      int             0=no, 1=yes
% HDR.FLAG.REFERENCE      string          COM, CAR: common average reference; LOC,LAR local average ref; LAP Laplacian derivation, WGT weighted average
% HDR.Classlabel          int             0: left, 1: right, etc.
% HDR.ID.Doctor                           Identification of doctor
% HDR.ID.Hospital                         Identification of Hospital
% HDR.Patient.Name                        Name of Patient
% HDR.Patient.Age                         Age of Patient
% HDR.Patient.Sex                         Patient Gender
% HDR.Patient.Handedness                  Patient Handedness
% HDR.Patient.Medication                  Medication
% HDR.Patient.Classification              Classification of Patient


text:
% -------------------------------------------------------------------------------------------------
function [bb] = get_minmax_BB(region)
%get_minmax_BB computes minmax bbox from the rotated one (REGION)
% -------------------------------------------------------------------------------------------------

bb = [min(region(:,1:2:end),[],2),...
    min(region(:,2:2:end),[],2),...
    max(region(:,1:2:end),[],2),...
    max(region(:,2:2:end),[],2)]-1;

end

author: foolwood

repo: DCFNet

save_path: github-repos/MATLAB/foolwood-DCFNet


text:
% The observation function:
%
%               function [F] = g_CaBBI(Xt,Phi,I_inp,inG)
%
% this function maps [Ca2+] kinetics to fluorescence observations (F)
% through a non-linear saturating function

function [F] = g_CaBBI(Xt,Phi,I_inp,inG)

ind = inG.ind;
Ca = Xt(ind);                        % [Ca2+] kinetics
d_F =  Phi(1);                       % offset parameter
Kd =  200;                           % dissociation constant
scale = inG.k_F0*exp(Phi(2));        % scale paraemeter
F =  scale*(Ca/(Ca + Kd)) + d_F;     % nonlinear mapping from [Ca2+] to F

end

author: MBB-team

repo: VBA-toolbox


text:
function p = fixOperatorProperties(p)

for i = 1:length(p.evalMap)
    % y = f(x)
    y = p.evalVariables(i);
    x = p.evalMap{i}.variableIndex;
    xL = p.lb(x);
    xU = p.ub(x);

    properties = p.evalMap{i}.properties;
    % Can convexity be fixed for xL <= x <= xL?
    if any(xL<xU) && (isa(properties.convexity,'function_handle') || ~isempty(properties.inflection) || strcmp(properties.convexity,'none'))
        vexity = 'none';
        % User-supplied code to derive convexity
        if isa(properties.convexity,'function_handle')
            vexity = properties.convexity(xL,xU);
        end
        % No such code, or code does not work as expected? Check inflection
        if isequal(vexity,'none') && ~isempty(properties.inflection)
            vexity = DeriveVexityFromInflection(properties,xL,xU);
        end
        if ~isequal(vexity,'none')
            % We have managed to fix convexity
            properties.convexity = vexity;
            % Since convexity is known, we can append a convex hull method
            if isempty(properties.convexhull) && ~isempty(properties.derivative)
                if isequal(vexity,'convex')
                    %f0 = @(x)real(eval([p.evalMap{i}.fcn '(x)']));
                    f0 = p.evalMap{i}.properties.function;
                    f = @(xL,xU)createConvexHullMethodConvex(xL,xU,f0,properties.derivative);
                elseif isequal(vexity,'concave')
                    %f0 = @(x)real(eval([p.evalMap{i}.fcn '(x)']));
                    f0 = p.evalMap{i}.properties.function;
                    f = @(xL,xU)createConvexHullMethodConcave(xL,xU,f0,properties.derivative);
                else
                    f = [];
                end
                properties.convexhull = f;
            end
        end
    end
    if any(xL<xU)
        if isa(properties.monotonicity,'function_handle')
            mono = properties.monotonicity(xL,xU);
            if ~isequal(mono,'none')
                properties.monotonicity = mono;
            end
        elseif ~isempty(properties.stationary) && (isequal(properties.shape,'bell-shape') || isequal(properties.shape,'v-shape'))
            mono = DeriveMonotonicityFromShape(properties,xL,xU);
            if ~isequal(mono,'none')
                properties.monotonicity = mono;
            end
        end
    end
    % Save possibly updated properties
    p.evalMap{i}.properties = properties;

    if isempty(p.evalMap{i}.properties.f_upper)
        if isequal(p.evalMap{i}.properties.convexity,'convex')
            %f = makefunction(p.evalMap{i}.fcn);
            f = p.evalMap{i}.properties.function;
            f_upper = @(z,xL,xU)(f(xL) + (z-xL)*(f(xU)-f(xL))/(xU-xL));
            df_upper = @(z,xL,xU)((f(xU)-f(xL))/(xU-xL));
            p.evalMap{i}.properties.f_upper = f_upper;
            p.evalMap{i}.properties.df_upper = df_upper;
        end
    end
    if isempty(p.evalMap{i}.properties.f_lower)
         if isequal(p.evalMap{i}.properties.convexity,'concave')
            %f = makefunction(p.evalMap{i}.fcn);
            f = p.evalMap{i}.properties.function;
            f_lower = @(z,xL,xU)(f(xL) + (z-xL)*(f(xU)-f(xL))/(xU-xL));
            df_lower = @(z,xL,xU)((f(xU)-f(xL))/(xU-xL));
            p.evalMap{i}.properties.f_lower = f_lower;
            p.evalMap{i}.properties.df_lower = df_lower;
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function[variables] = replaceVariables(obj, variables)
% [variables] = replaceVariables(obj, variables)
%
% Replace all default variables (x\d+) with free variables.
%
% Copyright by Holger Caesar, 2015

for i = 1 : numel(variables),
    oldVariable = variables{i};

    if regexp(oldVariable, 'x\d+'),
        freeVariable = getFreeVariable(obj);
        variables{i} = freeVariable;
    end;
end;

author: nightrome

repo: matconvnet-calvin


text:
function W = nullspaceLUSOLapply2Modes(mode, m, n, V, nullS)
% Computes the matrix vector product with the operator nullspace
% function handle of the form `y = pdMat(mode, m, n, x)`
%
% USAGE:
%
%    W = nullspaceLUSOLapply2Modes(mode, m, n, V, nullS)
%
% INPUTS:
%    mode:     :math:`mode = 1` returns :math:`W = Z V`, `mode = 2` returns :math:`W = Z^T V`
%    m:        first dimension of the matrix
%    n:        second dimension of the matrix
%    V:        one of the components of the multiplication
%    nullS:    structure `nullS` from the function `nullspaceLUSOLform(S)`;
%              where `m x n` sparse matrix `S` (:math:`m < n`).
%
% OUTPUT:
%    W:        Matrix vector product
%
% .. 16 May 2008: (MAS) First version of nullspaceLUSOLapply.m. See nullspaceLUSOLtest.m for testing.
% .. 13 Mar 2009: (RF) Tried to add matrix vector product of the transpose for use with pdco

if mode==1
    %returns   W = Z*V (mode=1)

    % Second, if V is an (n-m) x k sparse matrix (k >= 1),
    %        W = nullspaceLUSOLapply(nullS,V);
    % computes an n x k sparse matrix W from V such that S*W = 0.
    %
    % This is an operator form of finding an n x (n-m) matrix Z
    % such that S*Z = 0 and then computing W = Z*V.
    % The aim is to obtain W without forming Z explicitly.

    % 16 May 2008: (MAS) First version of nullspaceLUSOLapply.m.
    %              See nullspaceLUSOLtest.m for testing.

    Cinv    = nullS.Cinv; % Column scales
    L       = nullS.L;    % Strictly triangular L
    p       = nullS.p;    % Column permutation for S
    rankS   = nullS.rank; % rank(S)

    [n ,n ] = size(L);
    [mV,nV] = size(V);

    B       = [sparse(rankS,nV)
                V        ];
    Z       = (L')\B;
    W       = Z;
    W(p,:)  = Z;
    W       = Cinv*W;
else
    %returns  W = Z'*V (mode=2).

    % This is an operator form of finding an n x (n-m) matrix Z
    % such that S*Z = 0 and then computing W = Z'*V.
    % The aim is to obtain W without forming Z explicitly.

    % When S is scaled as S = R*A*C, we have Z = Cinv*P'*inv(L')[0;I]

    % This was the maths I tried to follow -RF
    %     W   =   A'*V
    %         =   (Cinv*P'*inv(L')[0;I])'*V
    %         =   [0,I]'*inv(L)*P*Cinv'*V

    Cinv    = nullS.Cinv; % Column scales
    L       = nullS.L;    % Strictly triangular L
    p       = nullS.p;    % Column permutation for S
    rankS   = nullS.rank; % rank(S)

    [n ,n ] = size(L);

    V = Cinv'*V;

    % I am not sure how to use these permutations
    % e.g. if P*A = B then is it B(p,:)=A or B=A(p,:)
    V(p,:) = V;
%     V(fliplr(p),:) = V;

    W = L\V;

    W = W(rankS+1:n,:);

end

return


text:
function sgmga_aniso_balance_test ( alpha_max, dim_num, level_weight )

%*****************************************************************************80
%
%% SGMGA_ANISO_BALANCE_TEST calls SGMGA_ANISO_BALANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_ANISO_BALANCE_TEST\n' );
  fprintf ( 1, '  ALPHA_MAX = %f\n', alpha_max );
  fprintf ( 1, '  Input weight sum: %f\n', sum ( level_weight(1:dim_num) ) );
  for dim = 1 : dim_num
    fprintf ( 1, '  %12f', level_weight(dim) );
  end
  fprintf ( 1, '\n' );

  level_weight2 = sgmga_aniso_balance ( alpha_max, dim_num, level_weight );

  fprintf ( 1, '  Output weight sum: %f\n', sum ( level_weight2(1:dim_num) ) );
  for dim = 1 : dim_num
    fprintf ( 1, '  %12f', level_weight2(dim) );
  end
  fprintf ( 1, '\n' );

  return
end


text:
function [Z,L,obj,err,iter] = latlrr(X,lambda,opts)

% Solve the Latent Low-Rank Representation by M-ADMM
%
% min_{Z,L,E} ||Z||_*+||L||_*+lambda*loss(E),
% s.t., XZ+LX-X=E.
% loss(E) = ||E||_1 or 0.5*||E||_F^2 or ||E||_{2,1}
% ---------------------------------------------
% Input:
%       X       -    d*n matrix
%       lambda  -    >0, parameter
%       opts    -    Structure value in Matlab. The fields are
%           opts.loss       -   'l1' (default): loss(E) = ||E||_1
%                               'l2': loss(E) = 0.5*||E||_F^2
%                               'l21': loss(E) = ||E||_{2,1}
%           opts.tol        -   termination tolerance
%           opts.max_iter   -   maximum number of iterations
%           opts.mu         -   stepsize for dual variable updating in ADMM
%           opts.max_mu     -   maximum stepsize
%           opts.rho        -   rho>=1, ratio used to increase mu
%           opts.DEBUG      -   0 or 1
%
% Output:
%       Z       -    n*n matrix
%       L       -    d*d matrix
%       E       -    d*n matrix
%       obj     -    objective function value
%       err     -    residual
%       iter    -    number of iterations
%
% version 1.0 - 19/06/2016
%
% Written by Canyi Lu (canyilu@gmail.com)
%

tol = 1e-8;
max_iter = 500;
rho = 1.1;
mu = 1e-4;
max_mu = 1e10;
DEBUG = 0;
loss = 'l1';

if ~exist('opts', 'var')
    opts = [];
end
if isfield(opts, 'loss');        loss = opts.loss;            end
if isfield(opts, 'tol');         tol = opts.tol;              end
if isfield(opts, 'max_iter');    max_iter = opts.max_iter;    end
if isfield(opts, 'rho');         rho = opts.rho;              end
if isfield(opts, 'mu');          mu = opts.mu;                end
if isfield(opts, 'max_mu');      max_mu = opts.max_mu;        end
if isfield(opts, 'DEBUG');       DEBUG = opts.DEBUG;          end

eta1 = 1.02*2*norm(X,2)^2; % for Z
eta2 = eta1; % for L
eta3 = 1.02*2; % for E

[d,n] = size(X);
E = zeros(d,n);
Z = zeros(n,n);
L = zeros(d,d);
Y = E;

XtX = X'*X;
XXt = X*X';

iter = 0;
for iter = 1 : max_iter
    Lk = L;
    Ek = E;
    Zk = Z;
    % first super block {Z}
    [Z,nuclearnormZ] = prox_nuclear(Zk-(X'*(Y/mu+L*X-X-E)+XtX*Z)/eta1,1/(mu*eta1));
    % second super block {L,E}
    temp = Lk-((Y/mu+X*Z-Ek)*X'+Lk*XXt-XXt)/eta2;
    [L,nuclearnormL] = prox_nuclear(temp,1/(mu*eta2));
    if strcmp(loss,'l1')
        E = prox_l1(Ek+(Y/mu+X*Z+Lk*X-X-Ek)/eta3,lambda/(mu*eta3));
    elseif strcmp(loss,'l21')
        E = prox_l21(Ek+(Y/mu+X*Z+Lk*X-X-Ek)/eta3,lambda/(mu*eta3));
    elseif strcmp(loss,'l2')
        E = (Y+mu*(X*Z+Lk*X-X+(eta3-1)*Ek))/(lambda+mu*eta3);
    else
        error('not supported loss function');
    end

    dY = X*Z+L*X-X-E;
    chgL = max(max(abs(Lk-L)));
    chgE = max(max(abs(Ek-E)));
    chgZ = max(max(abs(Zk-Z)));
    chg = max([chgL chgE chgZ max(abs(dY(:)))]);
    if DEBUG
        if iter == 1 || mod(iter, 10) == 0
            obj = nuclearnormZ+nuclearnormL+lambda*comp_loss(E,loss);
            err = norm(dY,'fro')^2;
            disp(['iter ' num2str(iter) ', mu=' num2str(mu) ...
                    ', obj=' num2str(obj) ', err=' num2str(err)]);
        end
    end

    if chg < tol
        break;
    end
    Y = Y + mu*dY;
    mu = min(rho*mu,max_mu);
end
obj = nuclearnormZ+nuclearnormZ+lambda*comp_loss(E,loss);
err = norm(dY,'fro')^2;

function out = comp_loss(E,loss)

switch loss
    case 'l1'
        out = norm(E(:),1);
    case 'l21'
        out = 0;
        for i = 1 : size(E,2)
            out = out + norm(E(:,i));
        end
    case 'l2'
        out = 0.5*norm(E,'fro')^2;
end

author: canyilu

repo: LibADMM-toolbox


text:
function [ meta ] = read_sicd01_meta( filename )
%READ_SICD01_META Read metadata from Sensor Independent Complex Data (SICD)
%file, version 0.1
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

% Open file
fid=fopen(filename,'r','b');
if(fid<0)
    error('READ_SICD01_META:InvalidFile','Unable to open file.');
end
% Read in non-xml header data
meta.format=fgetl(fid);
if(~strcmp(meta.format,'SICD-E/0.1'))
    error('READ_SICD01_META:UnrecognizedVersion','Unrecognized version of SICD format.  Only version 0.1 supported.');
end
for i=1:10
    [fieldname,value]=strtok(fgetl(fid),':=');
    fieldname=regexprep(fieldname,'-','_'); % Make sure field is valid Matlab variable name
    meta.header.(fieldname)=value(3:end);
    doubval=str2double(meta.header.(fieldname));
    if(~isnan(doubval)) % If value is numeric, store as such
        meta.header.(fieldname)=doubval;
    end
end
meta.header_length_in_bytes=ftell(fid);
fread(fid,2,'uchar'); % Blank line (/12/10) in between header and xml
% Read XML
SICD_xml_string = fread(fid,meta.header.xml_content_length,'uchar=>char',0)';
meta.xmldata = xml2simplestruct( xmlread( java.io.StringBufferInputStream( SICD_xml_string ) ) );
fclose(fid);

end

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

author: ngageoint


text:
%{
This file shows an example of using the full power of TFOCS. If you
are not familiar with basic usage of TFOCS, please see other demos
first.

We make an example that uses:
    -several variables (2), and the variables are matrices, not just vectors
    -several constraints (4)
    -affine operators (linear plus offset), and two of the operators
        are matrix --> matrix operators
    -debug mode

For small complicated examples like this, TFOCS is not necessarily
faster than software like CVX, since it takes more iterations than an
interior point method and there is some overhead in the TFOCS software
since the software is meant for flexibility rather than absolute speed.

However, if you take any complicated problem and scale the size by a factor
of 100, then CVX won't be able to handle it at all. TFOCS will do
just fine (and it might even be less than 100x as slow, since now the overhead
is not significant).

The problem we will solve:

min_{X1, X2} smooth1(X1)+smooth2(X2) +  sum_{i=1}^4  g_i( A1_i*X1 + A2_i*X2 + B_i )

meaning that we have 2 variables (X1 and X2), both of which are matrices,
and each variable has its own smooth function.
For non-smooth and/or constraint terms (indexed by "i"), we have 4 functions
( i = 1,2,3,4) g_i, and each has it's own affine operator in X1 and X2.
So each of the 4 affine operators has three parts: the portion linear in X1,
the portion linear in X2, and the constant offset "B_i".

To be explicit, the smooth functions are linear and quadratic, resp.:

smooth1(X1) = dot( s1, X1 ) + 3.4                  ("s1" is a matrix the same size as X1)
smooth2(X2) = dot( X2, X2 ) + dot( s2, X2 ) + 4.5  ("s2" is a matrix the same size as X2)

and the non-smooth/constraint terms are:

g_1(z)  = indicator set of the positive orthant of R^10
g_2(z)  = ||z||_2 (usual Euclidean norm) in R^15
g_3(z)  = ||z||_1 (l1 norm) in R^{10 x 20 }. This views a 10 x 20 matrix as a 200 x 1 vector.
g_4(z)  = indicator set of positive orthant of R^{20 x 22 }, i.e. each element must be >= 0 coordinate-wise

The affine operators are picked arbitrary, but some of them (#3 and #4) have a range that
is a set of matrices, rather than a set of vectors, in order to make
this more interesting.

%}

fileName = fullfile('reference_solutions','complicatedProblem1');
randn('state',29324);
rand('state',9332);
% rng(3481); % this only works on new versions of Matlab

% -- Variables --
% Two sets of variables, X1 (a matrix of size n1 x n2 ) and X2 (matrix of N1 x N2 )
n1 = 10; n2 = 20;
N1 = 12; N2 = 18;
X1 = zeros( n1, n2 );
X2 = zeros( N1, N2 );

% -- the smooth terms --
%   Note: the inner product used is the matrix inner product
%       that induces the Frobenius norm.
s1      = randn(n1,n2);
smooth1 = smooth_linear( s1, 3.4 );
s2      = randn(N1,N2);
smooth2 = smooth_quad( 1, s2, 4.5 );
% and the same thing, but in a format that CVX likes:
smoothF = @(X1,X2) vec(X1)'*vec(s1) + 3.4 + vec(X2)'*vec(s2) + 4.5 + ...
    sum_square(vec(X2))/2;

% -- some proximal terms --

nProx   = 4;
prox1   = proj_Rplus;
prox2   = proj_l2;          % primal is norm( , 2)
prox3   = proj_linf(1);     % this can take matrix varaibles...
prox4   = proj_Rplus;       % this can take matrix variables...
% Sizes of proxes (i.e. sizes of dual variables, i.e. size of range of linear terms)
d1      = [ 10, 1  ];
d2      = [ 15, 1  ];
d3      = [ n1, n2 ];
d4      = [ 20, 22 ];

% -- and linear terms --

% for prox1:
const1  = randn(d1);
temp1   = randn( prod(d1), n1*n2);
A1_X1   = linop_compose( linop_matrix(temp1), linop_vec([n1,n2]) );
temp2   = randn( prod(d1), N1*N2);
A1_X2   = linop_compose( linop_matrix(temp2), linop_vec([N1,N2]) );

A1      = @(X1,X2) temp1*vec(X1) + temp2*vec(X2) + const1;

% for prox2: (matrix variable)
const2  = randn(d2);
temp1   = randn( prod(d2), n1*n2);
A2_X1   = linop_compose( linop_matrix(temp1), linop_vec([n1,n2]) );
temp2   = randn( prod(d2), N1*N2);
A2_X2   = linop_compose( linop_matrix(temp2), linop_vec([N1,N2]) );

A2      = @(X1,X2) temp1*vec(X1) + temp2*vec(X2) + const2;

% for prox3:
const3  = 0;
A3_X1   = 63.4;     % this represents abstract scaling, i.e. any size input
A3_X2   = 0;        % this reprsents the zero linear operator

A3      = @(X1,X2) A3_X1*X1;

% for prox4: (matrix variable)
const4  = randn(d4);
temp1   = randn( prod(d4), n1*n2);
A4_X1   = linop_compose( linop_matrix(temp1), linop_vec([n1,n2]) );
rs      = linop_adjoint( linop_vec(d4) );
A4_X1   = linop_compose( rs, A4_X1 );
mat1    = @(x) reshape( x, d4(1), d4(2) );

temp2   = randn( prod(d4), N1*N2);
A4_X2   = linop_compose( linop_matrix(temp2), linop_vec([N1,N2]) );
rs      = linop_adjoint( linop_vec(d4) );
A4_X2   = linop_compose( rs, A4_X2 );

A4      = @(X1,X2) mat1( temp1*vec(X1) + temp2*vec(X2) ) + const4;

% -- set the smoothing parameter --
mu = 1;

if exist([fileName,'.mat'],'file')
    load(fileName); % contains X_CVX
    fprintf('Loaded problem from %s\n', fileName );
else
    % Get reference solution in CVX
    % First, get a solution to the smoothed version:
    cvx_begin
        variables X1(n1,n2) X2(N1,N2)
        % it's important that we use norm(vec(...),1) and not just norm(...,1),
        %   otherwise CVX interprets this with the wrong implicit inner product
        minimize(     smoothF( X1, X2 ) + ...
            norm( A2(X1,X2), 2 ) + norm( vec(A3(X1,X2)) , 1 ) + ...
            mu*( sum_square(vec(X1)) + sum_square(vec(X2)) )/2       )
        subject to
            A1(X1,X2) >= 0  % constraint is dual of prox1
            A4(X1,X2) >= 0  % constraint is dual of prox2
    cvx_end
    X_CVX_smoothed{1} = X1;
    X_CVX_smoothed{2} = X2;

    % Second, get a solution to the unsmoothed version
    cvx_begin
        variables X1(n1,n2) X2(N1,N2)
        minimize(     smoothF( X1, X2 ) + ...
            norm( A2(X1,X2), 2 ) + norm( vec(A3(X1,X2)) , 1 ) )
        subject to
            A1(X1,X2) >= 0
            A4(X1,X2) >= 0
    cvx_end
    X_CVX{1} = X1;
    X_CVX{2} = X2;

    save(fileName,'X_CVX', 'X_CVX_smoothed');
    fprintf('Saved data to file %s\n', fileName);
end

% Verify constraint are satisfied
fprintf('Is A1(x1,x2) >= 0? The min element is %g\n', min(    A1(X_CVX{1},X_CVX{2} )) )
fprintf('Is A2(x1,x2) >= 0? The min element is %g\n', min(min(A4(X_CVX{1},X_CVX{2}))) )
%% before running TFOCS, scale the problem perhaps?
% This is one way to see how big the norms are:
% nrm11    = linop_test( A1_X1 ); % 15.7
% nrm12    = linop_test( A1_X2 ); % 16.6
%
% nrm21    = linop_test( A2_X1 ); % 17.8
% nrm22    = linop_test( A2_X2 ); % 18.6
%
% nrm31    = linop_test( A3_X1 ); % 63.4
% nrm32    = linop_test( A3_X2 ); % 0
%
% nrm41    = linop_test( A4_X1 ); % 35.6
% nrm42    = linop_test( A4_X2 ); % 34.6

%% now, run TFOCS. First, try it without continuation

x0 = [];
z0 = [];
opts    = struct('continuation',false,'maxits',1500,'debug',true); % using 'debug' mode to print out useful information
opts.printEvery     = 50;

% Pick a scaling factor "s" (optional: set to "1" to have no effect)
s = .5; % helps a little bit
% (note that if we scale prox3 by "s", then we multiply the corresponding
%  affine part by "s", rather than divide them by "s", since prox3
%  is really for the dual and not the primal).

mu  = 1;
opts.errFcn{1} = @(f,d,p) norm( p{1}-X_CVX_smoothed{1}); % compare to smoothed reference solution
opts.errFcn{2} = @(f,d,p) norm( p{2}-X_CVX_smoothed{2});

[xAll,outParam,optsOut] = tfocs_SCD( {smooth1,smooth2}, ...
    { A1_X1, A1_X2, const1; A2_X1, A2_X2, const2; ...
    A3_X1*s, A3_X2*s, const3*s; A4_X1, A4_X2, const4 }, ...
    {prox1,prox2,prox_scale(prox3,s),prox4},...
    mu, x0, z0, opts );

mnConstraint1    = min(min( A1(xAll{1},xAll{2} ) ) );
mnConstraint2    = min(min( A4(xAll{1},xAll{2} ) ) );
fprintf('First constraint violated by:   %g\n', mnConstraint1);
fprintf('Second constraint violated by:  %g\n', mnConstraint2 );

% Check that we are within acceptable limits
er = sqrt(norm( xAll{1} - X_CVX_smoothed{1} )^2 + norm( xAll{2} - X_CVX_smoothed{2} )^2 )/...
    sqrt( norm(X_CVX_smoothed{1})^2 + norm(X_CVX_smoothed{2})^2 ); % should be about .006

if er > 0.04 || mnConstraint1 < -.01 || mnConstraint2 < -.01
    error('Failed the test');
else
    disp('This test successfully passed');
end

%% run TFOCS with continuation
x0 = [];
z0 = [];
% type "tfocs" at the command to see possible options
opts    = struct('continuation',true,'maxits',1500);
opts.printEvery     = 50;
opts.tol            = 1e-4;
opts.stopCrit       = 4;
opts.printStopCrit  = true;

% type "continuation" at the command to see possible options
contOpts    = struct( 'maxIts', 8 , 'muDecrement', 0.8 );
% ask for increased accuracy on the final solve
contOpts.finalTol = 1e-5;

s = .5; % helps a little bit
% (note that if we scale prox3 by "s", then we multiply the corresponding
%  affine part by "s", rather than divide them by "s", since prox3
%  is really for the dual and not the primal).

mu  = 10;
opts.errFcn{1} = @(f,d,p) norm( p{1}-X_CVX{1}); % compare to unsmoothed reference solution
opts.errFcn{2} = @(f,d,p) norm( p{2}-X_CVX{2});

[xAll,outParam,optsOut] = tfocs_SCD( {smooth1,smooth2}, ...
    { A1_X1, A1_X2, const1; A2_X1, A2_X2, const2; ...
    A3_X1*s, A3_X2*s, const3*s; A4_X1, A4_X2, const4 }, ...
    {prox1,prox2,prox_scale(prox3,s),prox4},...
    mu, x0, z0, opts, contOpts );

mnConstraint1    = min(min( A1(xAll{1},xAll{2} ) ) );
mnConstraint2    = min(min( A4(xAll{1},xAll{2} ) ) );
fprintf('First constraint violated by:   %g\n', mnConstraint1);
fprintf('Second constraint violated by:  %g\n', mnConstraint2 );

% Check that we are within acceptable limits
er = sqrt(norm( xAll{1} - X_CVX{1} )^2 + norm( xAll{2} - X_CVX{2} )^2 )/...
    sqrt( norm(X_CVX{1})^2 + norm(X_CVX{2})^2 ); % should be about .006

if er > 0.4 || mnConstraint1 < -.01 || mnConstraint2 < -.01
    error('Failed the test');
else
    disp('This test successfully passed');
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [c] = redblue(m)
%
% redblue  	Shades of red to white to blue colormap.
%		REDBLUE(M) returns an M-by-3 matrix containing a
%		"redblue" colormap.
%               CAVEAT this actually only works properly for even length
%               colormaps.  Some fiddling will get it to work generally.
%
%		See also HSV, COLORMAP, RGBPLOT.

if nargin<1
  [m,n] = size(colormap);
end;

m0=floor(m*0.0);
m1=floor(m*0.20);
m2=floor(m*0.20);
m3=floor(m/2)-m2-m1;

b_= [ 0.4*(0:m1-1)'/(m1-1)+0.6; ones(m2+m3,1)];
g_=  [zeros(m1,1); (0:m2-1)'/(m2-1);ones(m3,1)];
r_=  [zeros(m1,1); zeros(m2,1); (0:m3-1)'/(m3-1)];

r=[r_; flipud(b_)];
g=[g_; flipud(g_)];
b=[b_; flipud(r_)];

c=[r g b];


text:
% Example of use for function 'power_flow_solver'.
% Written by Dr. Yoash Levron, October 2012.
%
% All documentation may be obtained by typing:
% 'HELP power_flow_solver' at the matlab command prompt.
%
% This script defines a simple power network
% and uses the 'power_flow_solver' function to
% compute the power flow in it.
%
% Associated file: 'power_flow_solver.m'

%%%%%%% Define the network  %%%%%%%
% The network composes three buses.
% bus 1 - (slack bus) a voltage source, E=1.0
% bus 2 - a generator and a load.
% bus 3 - a load.
%
% Define the line impedances:
z12 = 0.05 + j*0.1;  % impedance (inductive) connecting bus 1 to 2.
z13 = 0.02 + j*0.05; % impedance (inductive) connecting bus 1 to 3.
z23 = j*0.05; % impedance (inductive) connecting bus 2 to 3.
z11 = -j*100;  % shunt impedance (capacitive) at bus 1
z22 = inf;  % no shunt impedance at bus 2
z33 = -j*40;  % shunt impedance (capacitive) at bus 3

% A matrix of the network impedances
Zmat = [z11 z12 z13;
        z12 z22 z23;
        z13 z23 z33];

% compute the admittance matrix (admittamce = 1/impedance):
Ymat = 1./Zmat;

% Define voltage at bus 1
E1_phase = 0; % degree (reference phase angle).
E1_mag = 1.0; % voltage magnitude
E1 = E1_mag*exp(j*E1_phase*(pi/180));

% Define generators and loads.
% Pg&Qg - generators.   Pl&Ql - loads.
Pg1=0;  Qg1=0; Pl1=0; Ql1=0; % no power source on the slack bus.
Pg2=2;  Qg2=0;  Pl2=1; Ql2=1;  % bus 2. generator and load.
Pg3=0;  Qg3=0;  Pl3=3; Ql3=0.5; % bus 3. a load.

% sum generators and load to form united power sources.
% (generators taken positive, loads taken negative).
P1 = Pg1 - Pl1;
P2 = Pg2 - Pl2;
P3 = Pg3 - Pl3;
Q1 = Qg1 - Ql1;
Q2 = Qg2 - Ql2;
Q3 = Qg3 - Ql3;

% Define power vectors:
Pbus = [P1 P2 P3];
Qbus = [Q1 Q2 Q3];

% run function to solve the power flow
[Ebus, Ibus, Imat, iter] = ...
    power_flow_solver(Ymat, Pbus, Qbus, E1);

% display results
clc;
disp('---  SUMMARY OF RESULTS  ---');
disp('number of iterations until convergence:');
disp(iter);
disp('voltage magnitudes');
disp(abs(Ebus.'));
disp('voltage phase angles (degree)');
disp(phase(Ebus.')*180/pi);
disp('current supplied by each source (amplitude)');
disp(abs(Ibus.'));
disp('Active power supplied by each source');
disp(real(Ebus.*conj(Ibus)));
disp('Reactive power supplied by each source');
disp(imag(Ebus.*conj(Ibus)));
disp('current in bus 1 --> bus 2');
disp(Imat(1,2));
disp('current in bus 1 --> bus 3');
disp(Imat(1,3));
disp('current in bus 2 --> bus 3');
disp(Imat(2,3));
disp('shunt currents in each bus (magnitude)');
disp(abs(diag(Imat).'));


text:
function y = permute( x, order )

%   Disciplined convex/geometric programming information for PERMUTE:
%      PERMUTE(A,ORDER) imposes no convexity restrictions on A. ORDER
%      must be constant.

%
% Determine the permutation
%

s = x.size_;
ndxs = 1 : prod( s );
ndx2 = permute( reshape( ndxs, s ), order );

%
% Permute the data
%

b = x.basis_;
try
    b = x.basis_( :, ndx2 );
catch %#ok
    ndxs( ndx2( : ).' ) = ndxs;
    [ r, c, v ] = find( b );
    b = sparse( r, ndxs( c ), v, size( b, 1 ), size( b, 2 ) );
    clear r c v
end
y = cvx( size( ndx2 ), b );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX