text:
classdef Par < handle
  % Par  TIC-TOC class for parfor loops
  %   This is a class for timing PARFOR loops. Once the measurement is
  %   made, it can be plotted to observe the various overheads that may
  %   exist in parallel for loops. It also shows the utilization of each
  %   worker.
  %
  % Par Properties:
  %  Worker       - Process ID of the worker
  %  ItStart      - Start time of iteration
  %  ItStop       - Stop time of iteration
  %
  % Par Methods:
  %  Par          - Construct Par objects and start the measurement
  %  stop         - Stop the measurement
  %  par2struct   - Convert the results to a flat structure
  %  plot         - Create a custom plot of the results
  %  horzcat      - Horizontally concatenate multiple PARFOR measurements
  %  StartTime    - Return the start time of the measurement
  %  StopTime     - Return the stop time of the measurement
  %  tic (Static) - Record start time of parfor loop iteration
  %  toc (Static) - Record stop time of parfor loop iteration
  %
  % This class should be used in the following way:
  %
  %       p = Par(n);                   (1)
  %       parfor id = 1:n
  %           Par.tic;                  (2)
  %
  %           <usual computations>
  %
  %           p(id) = Par.toc;          (3)
  %       end
  %       stop(p);                      (4)
  %
  %       plot(p);                      (5)
  %
  % 1. Construct a Par object, with the number of iterations as the input.
  %    This constructs the object.
  % 2. Call Par.tic just inside the PARFOR loop. This records the start
  %    time of the iteration.
  % 3. Call Par.toc just before the end of the PARFOR loop. Store the
  %    output to the appropriate index of the Par object. This is necessary
  %    for PARFOR to recognize that the variable is sliced.
  % 4. Stop the measurement. This records the final end time.
  % 5. Visualize.
  %
  % There may be some overhead in adding the Par construct. Because of
  % this, the numbers you get may not accurately portray the true timing,
  % especially for short computations.
  %
  % Doc page in Help browser
  %   doc Par

  % Copyright 2010 The MathWorks, Inc.

  properties (SetAccess = protected)
    % WORKER  Process ID of the worker running this iteration
    %   This is a read-only property and is only set using Par.tic and
    %   Par.toc
    Worker = NaN

    % ITSTART  Start time (from the start of object creation) of the iteration
    %   This is a read-only property and is only set using Par.tic and
    %   Par.toc
    ItStart = NaN

    % ITSTOP  Stop time (from the start of object creation) of the iteration
    %   This is a read-only property and is only set using Par.tic and
    %   Par.toc
    ItStop = NaN
  end

  properties (Hidden, GetAccess = protected)
    MainStartTime = []
    MainStopTime  = []
  end

  methods

    function obj = Par(n, varargin)
      % PAR Contructor for the Par object
      %
      % obj = Par(N)
      %   Creates an array of loop timer object for PARFOR loops. N is a
      %   positive integer equal to the length of iterations in the PARFOR
      %   loop.

      % Special case used internally
      %   check first since we want this to happen immediately.
      if nargin == 3
        obj.Worker  = n;
        obj.ItStart = varargin{1};
        obj.ItStop  = varargin{2};
        return;
      end

      switch nargin
        case 0
          return;

        case 1
          validateattributes(n, {'numeric'}, {'integer', 'scalar', 'positive'});
          obj(n) = Par();

          % Record the current time
          getnow(1);

          % Start the initial tic on this machine and on the workers
          tictoc(1);
          spmd
            tictoc(1);
          end

        otherwise
          error('Par:Par:InvalidNumberOfInputArguments', ...
            'You must call the constructor with one input argument: Par(n)');
      end

    end

    function newobj = horzcat(varargin)
      % HORZCAT Horizontally concatenate multiple PARFOR measurements.
      %
      % horzcat(p1, p2, ...) or [p1, p2, ...]
      %   Concatenates multiple Par objects. This allows you to combine
      %   results from multiple PARFOR blocks that are run sequentially.
      %   P1, P2, ... must be either Par objects or an empty array ([]).

      varargin(cellfun(@isempty, varargin)) = [];

      if ~all(cellfun('isclass', varargin, 'Par'))
        error('Par:horzcat:InvalidInputType', ...
          'Input arguments must [] or class type Par');
      end

      % Sort Par objects by MainStartTime
      [sortedGlobalTimes, idx] = sort(cellfun(@(x) min([x.MainStartTime]), varargin));
      % Calculate time to offset each Par objects
      offset = (sortedGlobalTimes - sortedGlobalTimes(1))*24*3600;

      % Reorder in ascending MainStartTime
      varargin = varargin(idx);

      % Create a new Par object with the appropriate number of elements
      num = sum(cellfun(@length, varargin));
      newobj = Par(num);

      % Concatenate one by one, offsetting the times by the offset values
      iCounter = 0;
      for iPar = 1:length(varargin)
        idx = (1:length(varargin{iPar})) + iCounter;
        os = offset(iPar);

        [newobj(idx).Worker]        = deal(varargin{iPar}.Worker);
        [newobj(idx).MainStartTime] = deal(varargin{iPar}.MainStartTime);
        [newobj(idx).MainStopTime]  = deal(varargin{iPar}.MainStopTime);

        newobj(idx(end)).MainStopTime = newobj(idx(end)).MainStopTime + os;

        newval = num2cell([varargin{iPar}.ItStart] + os);
        [newobj(idx).ItStart] = deal(newval{:});
        newval = num2cell([varargin{iPar}.ItStop] + os);
        [newobj(idx).ItStop] = deal(newval{:});

        iCounter = iCounter + length(varargin{iPar});
      end

    end

    function stop(obj)
      % STOP Stop the measurement.
      %
      % obj.stop() or stop(obj)
      %   Stops the measurement process. This should be called right after
      %   the parfor loop.

      if ~isempty(obj(end).MainStopTime)
        error('Par:stop:MultipleStopsNotAllowed', ...
          'You cannot stop more than once.');
      end

      % Assign only to the last element to minimize memory usage.
      obj(end).MainStartTime = getnow(2);
      obj(end).MainStopTime  = tictoc(2);                       %#ok<NASGU>

    end

    function val = StartTime(obj)
      % STARTTIME Return the end time of the measurement.
      %
      % val = obj.StartTime() or val = StartTime(obj)
      %   Returns the start time relative to the initial creation of the
      %   object. Usually, this is zero. If this object was a concatenation
      %   of multiple Par objects, then this will return multiple start
      %   times for each of the objects, relative to the initial start
      %   time.

      % The start time is stored in the hidden property MainStartTime
      gst = [obj.MainStartTime];
      % Convert to seconds (MainStartTime is a serial date)
      val = (gst - gst(1)) * 24 * 3600;

    end

    function val = StopTime(obj)
      % STOPTIME Return the end time of the measurement.
      %
      % val = obj.StopTime() or val = StopTime(obj)
      %   Returns the time passed from the creation of the object to
      %   stopping of the object.

      % The stop time is stored in the hidden property MainStopTime
      val = [obj.MainStopTime];

    end

    function parInfo = par2struct(obj)
      % PAR2STRUCT Convert the results to a flat structure.
      %
      % p = obj.par2struct() or p = par2struct(obj)
      %   Converts the results into a flat structure with the following
      %   fields:
      %       Start   : start times of each PARFOR block
      %       Stop    : end times of each PARFOR block
      %       Worker  : array of process id's for each iteration
      %       ItStart : array of start times for each iteration
      %       ItStop  : array of stop times for each iteration

      parInfo.Start   = obj.StartTime;
      parInfo.Stop    = obj.StopTime;
      parInfo.Worker  = [obj.Worker];
      parInfo.ItStart = [obj.ItStart];
      parInfo.ItStop  = [obj.ItStop];

    end

    function plot(obj, varargin)
      % PLOT Create a custom plot of the results.
      %
      % obj.plot() or plot(obj)
      %   Creates a custom plot of the results. The plot consists of 3
      %   parts:
      %     1) a line graph showing start and stop times of each iteration
      %     on each worker.
      %     2) a stacked bar graph showing the efficiency of each worker.
      %     3) a plot showing actual duration of each iteration.
      %
      % obj.plot(obj2) or plot(obj, obj2)
      %   Creates a custom plot that also compares the results of obj and
      %   obj2. obj2 must be of class Par and must be a result from a
      %   serial computation or a single-worker computation.

      % Error checking
      addSerial = 0;

      if nargin == 3

        obj2 = varargin{1};
        newMap = varargin{2};
        addSerial = 1;
        if ~isa(newMap, 'function_handle')
          error('Par:plot:invalidclass', ...
            'The third argument must be a function handle');
        end

      end

      if nargin == 2

        if isa(varargin{1}, 'function_handle')
          newMap = varargin{1};
        else
          obj2 = varargin{1};
          addSerial = 1;
          if ~(isa(obj2, 'Par'))
            error('Par:plot:invalidclass', ...
              'The second argument must also be of class Par or a function handle');
          end

          if length(unique([obj2.Worker])) > 1
            error('Par:plot:notserial', ...
              'The second argument must be the results from a serial (or 1 worker) calculation');
          end
        end

      end

      if any(isnan(obj.StopTime)) || (addSerial == 1 && isnan(obj2.StopTime))
        error('Par:plot:unterminated', ...
          'You need to stop() the object before you can plot()');
      end

      if any(isnan([obj.Worker])) || (addSerial == 1 && any(isnan([obj.Worker])))
        error('Par:plot:uncompleted', ...
          'You need to have timed all iterations before you can plot()');
      end

      % Call plotting subfunction
      if nargin == 1
        parPlot(addSerial, obj.par2struct)
      elseif nargin == 2
        if addSerial == 1
          parPlot(addSerial, obj.par2struct, obj2.par2struct)
        else
          parPlot(addSerial, obj.par2struct, newMap)
        end
      else
        parPlot(addSerial, obj.par2struct, obj2.par2struct, newMap)
      end

    end

  end

  methods (Static)

    function val = tic()
      % TIC Record start time of parfor loop iteration.
      %
      % Par.tic()
      %   Records the starting time within the loop. This should be called
      %   immediately inside the parfor loop.
      %
      % See also Par.toc

      persistent start_time

      if nargout
        if isempty(start_time)
          val = NaN;
        else
          val = start_time;
        end
      else
        start_time = tictoc(2);
      end

    end

    function obj = toc()
      % TOC Record stop time of parfor loop iteration.
      %
      % obj = Par.toc()
      %   Records the ending time within the loop. This should be called
      %   just before the end statement of the parfor loop.
      %
      % See also Par.tic

      if nargout

        workerNo = get(getCurrentTask,'ID');

        if isempty(workerNo)
          workerNo  = 0;
        end

        obj = Par(workerNo, Par.tic(), tictoc(2));

      end

    end

  end

  % These are overridden handle class methods. It's hidden so that it
  % doesn't show up as one of the standard methods in the DOC.
  methods (Hidden)

    function varargout = vertcat(varargin) %#ok<STOUT>

      error('Par:vertcat:VerticalConcatenationNotAllowed', ...
        ['Vertical concatenation is not allowed for a Par object.\n', ...
        'Consider horizontal concatenation.']);

    end

    function out = addlistener(varargin)
      out = addlistener@handle(varargin{:});
    end

    function delete(varargin)
      delete@handle(varargin{:});
    end

    function out = eq(varargin)
      out = eq@handle(varargin{:});
    end

    function findobj(varargin)
      findobj@handle(varargin{:});
    end

    function out = findprop(varargin)
      out = findprop@handle(varargin{:});
    end

    function out = ge(varargin)
      out = ge@handle(varargin{:});
    end

    function out = gt(varargin)
      out = gt@handle(varargin{:});
    end

    function out = le(varargin)
      out = le@handle(varargin{:});
    end

    function out = lt(varargin)
      out = lt@handle(varargin{:});
    end

    function out = ne(varargin)
      out = ne@handle(varargin{:});
    end

    function notify(varargin)
      notify@handle(varargin{:})
    end

  end

end

function out = getnow(flag)
% getnow
%   Internal function for recording global start time (in serial dates).

persistent nowtime

switch flag
  case 1
    nowtime = now;
  case 2
    if isempty(nowtime)
      error('Par:getnow:InvalidInitialization', ...
        ['The Par object may not have been initialized correctly.\n', ...
        'Be sure to initialize Par object with an input argument: p = Par(n)']);
    else
      out = nowtime;
    end
end

end

function out = tictoc(flag)
% tictoc
%   Internal function for recording start and end times.

persistent t1

switch flag
  case 1
    t1 = tic;
  case 2
    if isempty(t1)
      error('Par:tictoc:InvalidInitialization', ...
        ['The Par object may not have been initialized correctly.\n', ...
        'Be sure to initialize Par object with an input argument: p = Par(n)']);
    else
      out = toc(t1);
    end
  case 3
    if isempty(t1)
      error('Par:tictoc:InvalidInitialization', ...
        ['The Par object may not have been initialized correctly.\n', ...
        'Be sure to initialize Par object with an input argument: p = Par(n)']);
    else
      out = t1;
    end
end

end

function  parPlot(addSerial, parInfo, varargin)

cMap = @winter;

if nargin == 4
  parInfoS = varargin{1};
  cMap = varargin{2};
end

if nargin == 3
  if addSerial
    parInfoS = varargin{1};
  else
    cMap = varargin{1};
  end
end

% Determing unique workers
pid = parInfo.Worker;
allPIDs = unique(parInfo.Worker);
numWorkers = length(allPIDs);

% Preallocation
startTimes = cell(1, numWorkers);
endTimes = cell(1, numWorkers);
itTimes = cell(1, numWorkers);
totalTimes = nan(1, numWorkers);
overHead = nan(1, numWorkers);
overPerct = nan(1, numWorkers);
avgLoopTime = nan(1, numWorkers);

% Extracting start times, end times, and iteration times
if addSerial
  startTimesS = parInfoS.ItStart;
  endTimesS = parInfoS.ItStop;
  itTimesS = parInfoS.ItStop - parInfoS.ItStart;
  speedUp = parInfoS.Stop/parInfo.Stop;
end

totalParTime = sum(parInfo.Stop - parInfo.Start);

for ii = 1:numWorkers
  myIndex = (pid == allPIDs(ii));

  startTimes{ii} = parInfo.ItStart(myIndex);
  endTimes{ii} = parInfo.ItStop(myIndex);
  itTimes{ii} = parInfo.ItStop(myIndex) - parInfo.ItStart(myIndex);
  totalTimes(ii) = sum(endTimes{ii}-startTimes{ii});

  overHead(ii) = totalParTime-totalTimes(ii);
  overPerct(ii) = (overHead(ii)/totalParTime)*100;
  avgLoopTime(ii) = totalTimes(ii)/length(endTimes{ii});
end

% Plotting Data

myColors = cMap(numWorkers+addSerial); % setting up color scheme

hFig = figure;
subplot(2, 2, 1:2, 'Parent', hFig);
box on;

if addSerial
  plotdata = [startTimesS;endTimesS;nan(1, length(startTimesS))];
  plotdata2 = startTimesS;
  plotdata3 = endTimesS;

  line(plotdata(:), zeros(1, numel(plotdata)), ...
    'color', myColors(1,:));
  hold on
  line(plotdata2(:), zeros(1, numel(plotdata2)), ...
    'marker','.','MarkerSize',5, 'LineStyle', 'none',...
    'color', myColors(1,:));
  hold on
  line(plotdata3(:), zeros(1, numel(plotdata3)), ...
    'marker','X','MarkerSize',5, 'LineStyle', 'none',...
    'color', myColors(1,:));

  speedStr=['Speed Up = ', num2str(speedUp)];
  text(parInfo.Stop + (0.005*parInfo.Stop),numWorkers,speedStr)
end

for ii = 1:numWorkers
  plotdata = [startTimes{ii};endTimes{ii};nan(1, length(startTimes{ii}))];
  plotdata2 = startTimes{ii};
  plotdata3 = endTimes{ii};

  line(plotdata(:), ii*ones(1, numel(plotdata)), ...
    'color', myColors(ii+addSerial,:));
  hold on
  line(plotdata2(:), ii*ones(1, numel(plotdata2)), ...
    'marker','.','MarkerSize',5, 'LineStyle', 'none', ...
    'color',myColors(ii+addSerial,:));
  hold on
  line(plotdata3(:), ii*ones(1, numel(plotdata3)), ...
    'marker','X','MarkerSize',5, 'LineStyle', 'none', ...
    'color', myColors(ii+addSerial,:));

end

% Plotting total start and stop times

line([parInfo.Start; parInfo.Start], ...
  repmat([0-addSerial; (numWorkers+1)], 1, length(parInfo.Start)),...
  'color', [.5 .5 .5]);
line([parInfo.Stop; parInfo.Stop], ...
  repmat([0-addSerial; (numWorkers+1)], 1, length(parInfo.Stop)),...
  'color', [.5 .5 .5]);
line([0 0], [0-addSerial, (numWorkers+1)], 'color', 'k');
ylim([0-addSerial (numWorkers+1)])
set(gca, 'YTick', (1-addSerial):numWorkers);
set(gca, 'TickDir', 'out');

if addSerial == 1
  % Setting 0 label to Serial for clarity
  ylabels = get(gca, 'YTickLabel') ;
  newlabels = char(zeros(length(ylabels),6));
  newlabels(:,1) = ylabels;
  newlabels(1,:) = 'Serial';
  set(gca, 'YTickLabel', newlabels) ;
end

title('Iterations in Time By Worker','fontweight','b')
xlabel('Time (s)')
ylabel('Worker')

% Plotting Efficency
subplot(2,2,3)

if numWorkers > 1

  % Giving bars different colors
  for ii = 1:numWorkers
    pp = nan(numWorkers, 2);
    pp(ii, :) = [100-overPerct(ii), overPerct(ii)];
    hB = barh(1:numWorkers, pp, 'stack');
    set(hB(1),'FaceColor',myColors(addSerial+ii,:));
    set(hB(2),'FaceColor',[0.5 0.5 0.5]);
    hold on;
  end

else
  hB = barh( [1 2] ,[100-overPerct, NaN; overPerct, NaN]','Stack');
  set(hB(1),'FaceColor',myColors(1,:));
  set(hB(2),'FaceColor',[0.5 0.5 0.5])
end

ylim([0, numWorkers+1])
set(gca, 'YTick', 1:numWorkers);

title('Worker Utilization','fontweight','b')
xlabel('% of Total Time on Worker')
ylabel('Worker')

% Plotting Stats on Iteration Times
% Adding menu to change between linear and semilog axis
hMenu = uicontextmenu();
hAx = subplot(2,2,4, 'UIContextMenu', hMenu);
uimenu('Parent', hMenu, 'Label', 'Linear Scale', 'Checked', 'on', ...
  'Callback', {@setScaleFcn, hAx});
uimenu('Parent', hMenu, 'Label', 'Log Scale', 'Checked', 'off', ...
  'Callback', {@setScaleFcn, hAx});

if addSerial
  counter = (1:length(parInfo.ItStop));
  plot(counter,itTimesS,'color', ...
    myColors(1,:),'Marker','.','LineStyle','none')
  hold on
end

for ii = 1:numWorkers
  myIndex = (pid == allPIDs(ii));
  counter = (1:length(parInfo.ItStop));
  plot(counter(myIndex),itTimes{ii},'color', ...
    myColors(ii+addSerial,:),'Marker','*','LineStyle','none')
  hold on
end

title('Time Of Individual Iterations','fontweight','b')
xlabel('Iteration Number')
ylabel('Iteration Time (s)')

set(hFig, 'HandleVisibility', 'off');

% if numWorkers + addSerial > 1
%     % Add legend if more than one point
%     legendTxt = cell(addSerial+numWorkers,1);
%
%     for ii = 1-addSerial:numWorkers
%         legendTxt{ii+addSerial,1} = num2str(ii);
%     end
%
%     legend(legendTxt)
%
% end
end

function setScaleFcn(hObject, ~, hAx)

set(get(get(hObject, 'Parent'), 'Children'), 'Checked', 'off')
set(hObject, 'Checked', 'on');
switch get(hObject, 'Label')
  case 'Log Scale'
    set(hAx, 'YScale', 'log');
  case 'Linear Scale'
    set(hAx, 'YScale', 'linear');
end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function pts = getPoints(dens,ind)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% getPoints(P,ind)
%  returns the [Nd x Np] points of the nonparametric density estimate P
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2003 Alexander Ihler; distributable under GPL -- see README.txt

  if (nargin < 2) ind=1:dens.N; end;
  pts = zeros(dens.D,dens.N);
  pts(:,double(dens.perm(dens.N + (1:dens.N)))+1) = dens.centers(:,dens.N + (1:dens.N));
  pts = pts(:,ind);
%  pts = dens.centers(:,dens.N + ind);

author: ShapeNet

repo: RenderForCNN


text:
function [varargout]=bss_energy_ratios(F_s_target,F_e_interf,varargin)

% compute energy ratios corresponding to SDR/SIR/SNR/SAR given a decomposition of an estimated source into target/interference/noise/artifacts over frames.
%
% Usage:
%
%    [SDR,SIR,SAR]    =bss_energy_ratios(F_s_target,F_e_interf,F_e_artif)
%    [SDR,SIR,SNR,SAR]=bss_energy_ratios(F_s_target,F_e_interf,F_e_noise,F_e_artif)
%
% Input:
%   - F_s_target: n_frames x T matrix containing the frames of the target source contribution,
%   - F_e_interf: n_frames x T matrix containing the frames of the interferences contribution,
%   - F_e_noise: n_frames x T matrix containing the frames of the noise contribution (if any),
%   - F_e_artif: n_frames x T matrix containing the frames of the artifacts contribution.
%
% Ouput:
%   - SDR: n_frames x 1 vector contaning the Source to Distortion Ratios per frame,
%   - SIR: n_frames x 1 vector contaning the Source to Interferences Ratios per frame,
%   - SNR: n_frames x 1 vector contaning the Signal to Noise Ratios (if noise) per frame,
%   - SAR: n_frames x 1 vector contaning the Signal to Artifacts Ratios per frame.
%
% Developers:  - Cedric Fevotte (cf269@cam.ac.uk) - Emmanuel Vincent
% (vincent@ircam.fr) - Remi Gribonval (remi.gribonval@irisa.fr)

switch nargin
    case 3
        F_e_artif=varargin{1};
        % SDR
        F_e_total=F_e_interf+F_e_artif;
        varargout{1}= sum(F_s_target.^2,2)./sum(F_e_total.^2,2);
        % SIR
        varargout{2}=sum(F_s_target.^2,2)./sum(F_e_interf.^2,2);
        % SAR
        varargout{3}=sum((F_s_target+F_e_interf).^2,2)./sum(F_e_artif.^2,2);

    case 4
        F_e_noise=varargin{1};
        F_e_artif=varargin{2};
        % SDR
        F_e_total=F_e_interf+F_e_noise+F_e_artif;
        varargout{1}=sum(F_s_target.^2,2)./sum(F_e_total.^2,2);
        % SIR
        varargout{2}=sum(F_s_target.^2,2)./sum(F_e_interf.^2,2);
        % SNR
        varargout{3}=sum((F_s_target+F_e_interf).^2,2)./sum(F_e_noise.^2,2);
        % SAR
        varargout{4}=sum((F_s_target+F_e_interf+F_e_noise).^2,2)./sum(F_e_artif.^2,2);
end

author: posenhuang

repo: deeplearningsourceseparation


text:
function [ know, x ] = p06_sol ( n )

%*****************************************************************************80
%
%% P06_SOL returns the solution for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 1;

  x(1:n) = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function obj = updateBar3(obj, surfaceIndex)

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(surfaceIndex).AssociatedAxis);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,axIndex);

%-SURFACE DATA STRUCTURE- %
bar_data = get(obj.State.Plot(surfaceIndex).Handle);
figure_data = get(obj.State.Figure.Handle);

%-AXIS STRUCTURE-%
axis_data = get(ancestor(bar_data.Parent,'axes'));

%-GET SCENE-%
eval(['scene = obj.layout.scene' num2str(xsource) ';']);

%-------------------------------------------------------------------------%

%-associate scene-%
obj.data{surfaceIndex}.scene = sprintf('scene%d', xsource);

%-------------------------------------------------------------------------%

%-surface type-%
obj.data{surfaceIndex}.type = 'mesh3d';

%-------------------------------------------------------------------------%

%-FORMAT DATA-%
xdata = bar_data.XData;
ydata = bar_data.YData;
zdata = bar_data.ZData;
cdata = bar_data.CData;

%-parse xedges-%
xedges = xdata(2, 1:2:end);

%-parse yedges-%
yedges = ydata(2:6:end, 2);
yedges = [yedges', mean(diff(yedges(1:2)))];

%-parse values-%
values = [];
for n = 1:6:size(zdata, 1)
    values = [values, diff(zdata(n:n+1, 2))];
end

%-parse offsets-%
offsets = zdata(1:6:end, 2)';

%-------------------------------------------------------------------------%

%-get the values to use plotly's mesh3D-%
bargap = diff(yedges(1:2)) - diff(ydata(2:3));
[X, Y, Z, I, J, K] = get_plotly_mesh3d(xedges, yedges, values, bargap);

%---------------------------------------------------------------------%

%-reformat Z according to offsets-%
m = 1;
lz2 = 0.5*length(Z);

for n = 1:4:lz2
    Z(n:n+3) = Z(n:n+3)+offsets(m);
    Z(n+lz2:n+lz2+3) = Z(n+lz2:n+lz2+3)+offsets(m);
    m = m + 1;
end

%-------------------------------------------------------------------------%

%-set mesh3d data-%
obj.data{surfaceIndex}.x = X;
obj.data{surfaceIndex}.y = Y;
obj.data{surfaceIndex}.z = Z;
obj.data{surfaceIndex}.i = int16(I-1);
obj.data{surfaceIndex}.j = int16(J-1);
obj.data{surfaceIndex}.k = int16(K-1);

%-------------------------------------------------------------------------%

%-coloring-%
cmap = figure_data.Colormap;

if isnumeric(bar_data.FaceColor)

    %-paper_bgcolor-%
    col = 255*bar_data.FaceColor;
    col = sprintf('rgb(%f,%f,%f)', col);

else
    switch bar_data.FaceColor

        case 'none'
            col = 'rgba(0,0,0,0)';

        case {'flat','interp'}

            switch bar_data.CDataMapping

                case 'scaled'
                    capCD = max(min(cdata(1,1),axis_data.CLim(2)),axis_data.CLim(1));
                    scalefactor = (capCD - axis_data.CLim(1))/diff(axis_data.CLim);
                    col =  255*(cmap(1+ floor(scalefactor*(length(cmap)-1)),:));
                case 'direct'
                    col =  255*(cmap(cdata(1,1),:));

            end

            col = sprintf('rgb(%f,%f,%f)', col);

        case 'auto'
            col = 'rgb(0,113.985,188.955)';
    end
end

obj.data{surfaceIndex}.color = col;

%-------------------------------------------------------------------------%

%-some settings-%
obj.data{surfaceIndex}.contour.show = true;
obj.data{surfaceIndex}.contour.width = 6;
obj.data{surfaceIndex}.contour.color='rgb(0,0,0)';
obj.data{surfaceIndex}.flatshading = false;

%-------------------------------------------------------------------------%

%-lighting settings-%
obj.data{surfaceIndex}.lighting.diffuse = 0.8;
obj.data{surfaceIndex}.lighting.ambient = 0.65;
obj.data{surfaceIndex}.lighting.specular = 1.42;
obj.data{surfaceIndex}.lighting.roughness = 0.52;
obj.data{surfaceIndex}.lighting.fresnel = 0.2;
obj.data{surfaceIndex}.lighting.vertexnormalsepsilon = 1e-12;
obj.data{surfaceIndex}.lighting.facenormalsepsilon = 1e-6;

obj.data{surfaceIndex}.lightposition.x = 0;
obj.data{surfaceIndex}.lightposition.y = 0;
obj.data{surfaceIndex}.lightposition.z = 0;

%-------------------------------------------------------------------------%

%-surface name-%
obj.data{surfaceIndex}.name = bar_data.DisplayName;

%-------------------------------------------------------------------------%

%-surface visible-%
obj.data{surfaceIndex}.visible = strcmp(bar_data.Visible,'on');

%-------------------------------------------------------------------------%

leg = get(bar_data.Annotation);
legInfo = get(leg.LegendInformation);

switch legInfo.IconDisplayStyle
    case 'on'
        showleg = true;
    case 'off'
        showleg = false;
end

obj.data{surfaceIndex}.showlegend = showleg;

%-------------------------------------------------------------------------%

%-SETTING SCENE-%

%-------------------------------------------------------------------------%

%-aspect ratio-%
ar = obj.PlotOptions.AspectRatio;

if ~isempty(ar)
    if ischar(ar)
        scene.aspectmode = ar;
    elseif isvector(ar) && length(ar) == 3
        xar = ar(1);
        yar = ar(2);
        zar = ar(3);
    end
else

    %-define as default-%
    xar = max(xedges(:));
    yar = max(yedges(:));
    zar = 0.7*max([xar, yar]);
end

scene.aspectratio.x = xar;
scene.aspectratio.y = yar;
scene.aspectratio.z = zar;

%-------------------------------------------------------------------------%

%-camera eye-%
ey = obj.PlotOptions.CameraEye;

if ~isempty(ey)
    if isvector(ey) && length(ey) == 3
        scene.camera.eye.x = ey(1);
        scene.camera.eye.y = ey(2);
        scene.camera.eye.z = ey(3);
    end
else

    %-define as default-%
    scene.camera.eye.x = xar + 7;
    scene.camera.eye.y = yar - 2;
    scene.camera.eye.z = zar + 0.5;
end

%-------------------------------------------------------------------------%

%-axis configuration-%
scene.xaxis.range = axis_data.XLim(end:-1:1);
scene.yaxis.range = axis_data.YLim;
scene.zaxis.range = axis_data.ZLim;

scene.xaxis.tickvals = axis_data.XTick;
scene.xaxis.ticktext = axis_data.XTickLabel;

scene.yaxis.tickvals = axis_data.YTick;
scene.yaxis.ticktext = axis_data.YTickLabel;

scene.zaxis.tickvals = axis_data.ZTick;
scene.zaxis.ticktext = axis_data.ZTickLabel;

scene.xaxis.zeroline = false;
scene.yaxis.zeroline = false;
scene.zaxis.zeroline = false;

scene.xaxis.showline = true;
scene.yaxis.showline = true;
scene.zaxis.showline = true;

scene.xaxis.tickcolor = 'rgba(0,0,0,1)';
scene.yaxis.tickcolor = 'rgba(0,0,0,1)';
scene.zaxis.tickcolor = 'rgba(0,0,0,1)';

scene.xaxis.ticklabelposition = 'outside';
scene.yaxis.ticklabelposition = 'outside';
scene.zaxis.ticklabelposition = 'outside';

scene.xaxis.title = axis_data.XLabel.String;
scene.yaxis.title = axis_data.YLabel.String;
scene.zaxis.title = axis_data.ZLabel.String;

%-------------------------------------------------------------------------%

%-SET SCENE TO LAYOUT-%
obj.layout = setfield(obj.layout, sprintf('scene%d', xsource), scene);

%-------------------------------------------------------------------------%

end

function bar_ = bar_data(position3d, size_)
    % position3d - 3-list or array of shape (3,) that represents the point of coords (x, y, 0), where a bar is placed
    % size = a 3-tuple whose elements are used to scale a unit cube to get a paralelipipedic bar
    % returns - an array of shape(8,3) representing the 8 vertices of  a bar at position3d

    if nargin < 2
        size_ = [1, 1, 1];
    end

    bar_ = [...
            0, 0, 0; ...
            1, 0, 0; ...
            1, 1, 0; ...
            0, 1, 0; ...
            0, 0, 1; ...
            1, 0, 1; ...
            1, 1, 1; ...
            0, 1, 1 ...
            ]; % the vertices of the unit cube

    for n =1:size(bar_, 1)
        bar_(n,:) = bar_(n,:) .* size_; % scale the cube to get the vertices of a parallelipipedic bar_
    end

    bar_ = bar_ + position3d; %translate each  bar_ on the directio OP, with P=position3d
end

function [vertices, I, J, K] = triangulate_bar_faces(positions, sizes)
    % positions - array of shape (N, 3) that contains all positions in the plane z=0, where a histogram bar is placed
    % sizes -  array of shape (N,3); each row represents the sizes to scale a unit cube to get a bar
    % returns the array of unique vertices, and the lists i, j, k to be used in instantiating the go.Mesh3d class

    if nargin < 2
        sizes = ones(size(positions,1), 3); %[(1,1,1)]*len(positions)
    else
        sizes;
        % if isinstance(sizes, (list, np.ndarray)) and len(sizes) != len(positions):
        %     raise ValueError('Your positions and sizes lists/arrays do not have the same length')
    end

    c = 1;
    for n = 1:size(positions, 1)
        if sizes(n, 3) ~= 0
            all_bars(:,:,c) = bar_data(positions(n,:), sizes(n,:))';
            c = c+1;
        end
    end

    % all_bars = [bar_data(pos, size)  for pos, size in zip(positions, sizes) if size[2]!=0]
    [r, q, p] = size(all_bars);

    % extract unique vertices from the list of all bar vertices
    all_bars = reshape(all_bars, [r, p*q])';
    [vertices, ~, ixr] = unique(all_bars, 'rows');

    %for each bar, derive the sublists of indices i, j, k assocated to its chosen  triangulation
    I = [];
    J = [];
    K = [];

    for k = 0:p-1
        aux = ixr([1+8*k, 1+8*k+2,1+8*k, 1+8*k+5,1+8*k, 1+8*k+7, 1+8*k+5, 1+8*k+2, 1+8*k+3, 1+8*k+6, 1+8*k+7, 1+8*k+5]);
        I = [ I;  aux(:)];
        aux = ixr([1+8*k+1, 1+8*k+3, 1+8*k+4, 1+8*k+1, 1+8*k+3, 1+8*k+4, 1+8*k+1, 1+8*k+6, 1+8*k+7, 1+8*k+2, 1+8*k+4, 1+8*k+6]);
        J = [ J;  aux(:)];
        aux = ixr([1+8*k+2, 1+8*k, 1+8*k+5, 1+8*k, 1+8*k+7, 1+8*k, 1+8*k+2, 1+8*k+5, 1+8*k+6, 1+8*k+3, 1+8*k+5, 1+8*k+7]);
        K = [ K;  aux(:)];
    end

end

function [X, Y, Z, I, J, K] = get_plotly_mesh3d(xedges, yedges, values, bargap)
    % x, y- array-like of shape (n,), defining the x, and y-ccordinates of data set for which we plot a 3d hist

    xsize = xedges(2)-xedges(1)-bargap;
    ysize = yedges(2)-yedges(1)-bargap;
    [xe, ye]= meshgrid(xedges(1:end-1), yedges(1:end-1));
    ze = zeros(size(xe));

    positions = zeros([size(xe), 3]);
    positions(:,:,1) = xe;
    positions(:,:,2) = ye;
    positions(:,:,3) = ze;

    [m, n, p] = size(positions);
    positions = reshape(positions, [m*n, p]);

    h = values'; h = h(:);
    sizes = [];
    for n = 1:length(h)
        sizes = [sizes; ysize, ysize, h(n)];
    end

    [vertices, I, J, K]  = triangulate_bar_faces(positions, sizes);
    X = vertices(:,1);
    Y = vertices(:,2);
    Z = vertices(:,3);

end

author: plotly


text:
function export(t,fname,varargin)
% export tensor to json file
%
% Syntax
%   m = matrix(T)
%   m = matrix(T,'voigt')
%
% Input
%  T - @tensor
%
% Output
%  m - matrix
%
% Options
%  voigt - give a 4 rank tensor in voigt notation, i.e. as a 6 x 6 matrix
%
% See also
%

CS.abg = [t.CS.alpha,t.CS.beta,t.CS.gamma]./degree;
CS.abc = norm([t.CS.aAxis,t.CS.bAxis,t.CS.cAxis]);
CS.pointGroup = t.CS.pointGroup;
CS.mineral = t.CS.mineral;
CS.alignment = t.CS.alignment;
warning off
data = struct(t);
warning on
data.M = matrix(t,'voigt');
data.CS = CS;

savejson(class(t),data,fname);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [x,r,g,info] = spg_bp(A,b,varargin )
%SPG_BP  Solve the basis pursuit (BP) problem
%
%   SPG_BP is designed to solve the basis pursuit problem
%
%   (BP)  minimize  ||X||_1  subject to  AX = B,
%
%   where A is an M-by-N matrix, B is an M-vector, and SIGMA is a
%   nonnegative scalar.  In all cases below, A can be an explicit M-by-N
%   matrix or matrix-like object for which the operations  A*x  and  A'*y
%   are defined (i.e., matrix-vector multiplication with A and its
%   adjoint.)
%
%   Also, A can be a function handle that points to a function with the
%   signature
%
%   v = A(w,mode)   which returns  v = A *w  if mode == 1;
%                                  v = A'*w  if mode == 2.
%
%   X = SPG_BP(A,B) solves the BP problem.
%
%   X = SPG_BP(A,B,OPTIONS) specifies options that are set using
%   SPGSETPARMS.
%
%   [X,R,G,INFO] = SPG_BP(A,B,OPTIONS) additionally returns the
%   residual R = B - A*X (which should be small), the objective gradient G
%   = A'*R, and an INFO structure.  (See SPGL1 for a description of this
%   last output argument.)
%
%   See also spgl1, spgSetParms, spg_bpdn, spg_lasso.

%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
%   http://www.cs.ubc.ca/labs/scl/spgl1
%   $Id: spg_bp.m 1074 2008-08-19 05:24:28Z ewout78 $

if ~exist('b','var') || isempty(b)
    error('Second argument cannot be empty.');
end
if ~exist('A','var') || isempty(A)
    error('First argument cannot be empty.');
end

sigma = 0;
tau = 0;
x0  = [];
[x,r,g,info] = spgl1(A,b,tau,sigma,x0,varargin{:});

author: mpf


text:
% JN Kather 2018

function [lgraph,imageInputSize,networkType] = getAndModifyNet(nnmodel,hyperparam,numOutputClasses)

% load pre-trained network model for transfer learning
switch nnmodel
    case 'vgg19'
        rawnet = vgg19;
        networkType = 'series';
    case 'vgg16'
        rawnet = vgg16;
        networkType = 'series';
    case 'alexnet'
        rawnet = alexnet;
        networkType = 'series';
    case 'inceptionv3'
        rawnet = inceptionv3;
        networkType = 'DAG';
        layersForRemoval = {'predictions', 'predictions_softmax','ClassificationLayer_predictions'};
        layersForReconnection = {'avg_pool','fc'};
    case 'googlenet'
        rawnet = googlenet;
        networkType = 'DAG';
        layersForRemoval = {'loss3-classifier','prob','output'};
        layersForReconnection = {'pool5-drop_7x7_s1','fc'};
    case 'resnet18'
        rawnet = resnet18;
        networkType = 'DAG';
        layersForRemoval = {'fc1000', 'prob','ClassificationLayer_predictions'};
        layersForReconnection = {'pool5','fc'};
    case 'resnet50'
        rawnet = resnet50;
        networkType = 'DAG';
        layersForRemoval = {'ClassificationLayer_fc1000', 'fc1000_softmax','fc1000'};
        layersForReconnection = {'avg_pool','fc'};
    case 'resnet101'
        rawnet = resnet101;
        networkType = 'DAG';
        layersForRemoval = {'fc1000', 'prob','ClassificationLayer_predictions'};
        layersForReconnection = {'pool5','fc'};
    case 'squeezenet'
        rawnet = squeezenet;
        networkType = 'DAG';
        layersForRemoval = {'pool10', 'prob','ClassificationLayer_predictions'};
        layersForReconnection = {'relu_conv10','fc'};
    case 'inceptionresnetv2'
        rawnet = inceptionresnetv2;
        networkType = 'DAG';
        layersForRemoval = {'predictions', 'predictions_softmax','ClassificationLayer_predictions'};
        layersForReconnection = {'avg_pool','fc'};
    otherwise
        error('wrong network model specified');
end

% prune and rewire network
switch networkType
    case 'series' % e.g. alexnet
        lgraph = rawnet.Layers;
        % freeze shallow layers
        freezeIndex = 1:(numel(lgraph)-hyperparam.hotLayers);
        lgraph(freezeIndex) = freezeWeights(lgraph(freezeIndex));
        % overwrite penultimate and last layer
        lgraph(end-2) = fullyConnectedLayer(numOutputClasses,'Name','fc',...
            'WeightLearnRateFactor',hyperparam.learnRateFactor,...
            'BiasLearnRateFactor',hyperparam.learnRateFactor);
        lgraph(end) = classificationLayer;
        imageInputSize = lgraph(1).InputSize(1:2);
    case 'DAG' % e.g. googlenet
        % freeze shallow layers
        lgraph = layerGraph(rawnet); % convert network to layer graph
        layers = lgraph.Layers;      % extract layers
        connections = lgraph.Connections; % exctract connections
        freezeIndex = 1:(numel(layers)-hyperparam.hotLayers);
        layers(freezeIndex) = freezeWeights(layers(freezeIndex));
        lgraph = createLgraphUsingConnections(layers,connections);
        % remove old layers
        lgraph = removeLayers(lgraph,layersForRemoval);
        % add new layers and connect
        newLayers = [
            fullyConnectedLayer(numOutputClasses,'Name','fc',...
            'WeightLearnRateFactor', hyperparam.learnRateFactor,...
            'BiasLearnRateFactor', hyperparam.learnRateFactor),...
            softmaxLayer('Name','softmax'),...
            classificationLayer('Name','classoutput')];
        lgraph = addLayers(lgraph,newLayers);
        lgraph = connectLayers(lgraph,layersForReconnection{1},layersForReconnection{2});
        imageInputSize = lgraph.Layers(1).InputSize(1:2);
    otherwise, error('undefined network type');
end
end

author: jnkather

repo: MSIfromHE

save_path: github-repos/MATLAB/jnkather-MSIfromHE


text:
function RunningOneMicrostructureVademecum

%dSmooth = obtainSettings('SquareMesh4b','Rectangle');
dSmooth = obtainSettings('SmallCircleQ2','SmoothRectangle');

vc = VademecumCellVariablesCalculator(dSmooth);
vc.computeVademecumData()

end

function d = obtainSettings(prefix,freeFemFile)
d = SettingsVademecumCellVariablesCalculator();
d.freeFemFileName = freeFemFile;
d.fileName   = prefix;
d.mxMin = 0.1;
d.mxMax = 0.1;
d.myMin = 0.1;
d.myMax = 0.1;
d.nMx   = 1;
d.nMy   = 1;
d.outPutPath = [];
d.print = true;

%i = 4;
%d.freeFemSettings.hMax  = 10^(-1)/(2^(i - 1));%0.0025;

d.freeFemSettings.hMax  = 0.02;%0.0025;

%d.smoothingExponentSettings.type = 'Optimal';

d.smoothingExponentSettings.type = 'Given';
d.smoothingExponentSettings.q = 2;

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function startPreciseHandGuiding1( t,wightOfTool,COMofTool )
%% About:
% This function is used to start the handguiding functionality on the robot

%% Syntax:
% [ output_args ] = startPreciseHandGuiding1( t,wightOfTool,COMofTool )

%% Arreguments:
% t: is the TCP/IP connection object
% wightOfTool: weight of the tool connected to the flange in newton, unit Newtons.
% COMofTool: coordinates of the center of mass of the tool, descirbed in
% the reference frame of the flange (mm).

% Copy right, Mohammad SAFEEA, 22nd of Oct 2017

com=colVec(COMofTool)/1000; % convert from mm to meter
wot=-wightOfTool; % negtive, wight of tool described in base frame of robot

if(sum(size(com)==[3,1])==2)
else
    fprintf('Error: COMofTool vector is not in the right dimention \n');
    disp('Function termintated');
    return;
end

if norm(com)>0.5
    fprintf('Error: COMofTool vector shall not have a norm bigger than 0.5 meters \n');
    disp('Function termintated');
    return;
end

if(sum(size(wot)==[1,1])==2)
else
    fprintf('Error: wightOfTool variable is not in the right dimention \n');
    disp('Function termintated');
    return;
end

s='preciseHandGuiding1';
s=[s,'_',num2str(wot)];
for i=1:3
	s=[s,'_',num2str(com(i))];
end
theCommand=s;
fprintf(t, theCommand);

fprintf('Precise hand gudining functionality is started \n');
fprintf('To terminate the precise hand guding function, press the green button for more than 5 sec');
fprintf('\n keep pressing until the red light starts to flicker then release your hand, ')
fprintf('\n the hand guiding mode will be terminated ');

    message='';
    w = warning ('off','all');
    readingFlag=false;
    while readingFlag==false
        message=fgets(t);

        if checkAcknowledgment(message)
            readingFlag=true;
        end
    end
     w = warning ('on','all');
    fprintf('Hand gudining was terminated by user at robot side \n');
end

function y=colVec(x)
    if(size(x,2)==1)
        y=x;
    else
        y=x';
    end
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function edge = clipLine(line, box, varargin)
%CLIPLINE Clip a line with a box.
%
%   EDGE = clipLine(LINE, BOX);
%   LINE is a straight line given as a 4 element row vector: [x0 y0 dx dy],
%   with (x0 y0) being a point of the line and (dx dy) a direction vector,
%   BOX is the clipping box, given by its extreme coordinates:
%   [xmin xmax ymin ymax].
%   The result is given as an edge, defined by the coordinates of its 2
%   extreme points: [x1 y1 x2 y2].
%   If line does not intersect the box, [NaN NaN NaN NaN] is returned.
%
%   Function works also if LINE is a N-by-4 array, if BOX is a Nx4 array,
%   or if both LINE and BOX are N-by-4 arrays. In these cases, EDGE is a
%   N-by-4 array.
%
%
%   Example
%   line = [30 40 10 0];
%   box = [0 100 0 100];
%   res = clipLine(line, box)
%   res =
%       0 40 100 40
%
%   See also
%   lines2d, boxes2d, edges2d
%   clipEdge, clipRay, clipLine3d
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2007-08-27, using Matlab 7.4.0.287 (R2007a)
% Copyright 2007-2022 INRA - Cepia Software Platform

% adjust size of two input arguments
nLines = size(line, 1);
nBoxes = size(box, 1);
if nLines == 1 && nBoxes > 1
    line = repmat(line, nBoxes, 1);
elseif nBoxes == 1 && nLines > 1
    box = repmat(box, nLines, 1);
elseif nLines ~= nBoxes
    error('bad sizes for input');
end

% allocate memory
nLines = size(line, 1);
edge   = zeros(nLines, 4);

% main loop on lines
for i = 1:nLines
    % extract limits of the box
    xmin = box(i, 1);
    xmax = box(i, 2);
    ymin = box(i, 3);
    ymax = box(i, 4);

    % use direction vector for box edges similar to direction vector of the
    % line in order to reduce computation errors
    delta = hypot(line(i,3), line(i,4));

	% compute intersection with each edge of the box
    px1 = intersectLines(line(i,:), [xmin ymin delta 0]);   % lower edge
    px2 = intersectLines(line(i,:), [xmax ymin 0 delta]);   % right edge
    py1 = intersectLines(line(i,:), [xmax ymax -delta 0]);  % upper edge
    py2 = intersectLines(line(i,:), [xmin ymax 0 -delta]);  % left edge

    % remove undefined intersections (case of lines parallel to box edges)
    points = [px1 ; px2 ; py1 ; py2];
    points = points(isfinite(points(:,1)), :);

    % sort points according to their position on the line
    pos = linePosition(points, line(i,:));
    [pos, inds] = sort(pos); %#ok<ASGLU>
    points = points(inds, :);

    % create clipped edge by using the two points in the middle
    ind = size(points, 1)/2;
    inter1 = points(ind,:);
    inter2 = points(ind+1,:);
    edge(i, 1:4) = [inter1 inter2];

    % check that middle point of the edge is contained in the box
    midX = mean(edge(i, [1 3]));
    xOk = xmin <= midX && midX <= xmax;
    midY = mean(edge(i, [2 4]));
    yOk = ymin <= midY && midY <= ymax;

    % if one of the bounding condition is not met, set edge to NaN
    if ~(xOk && yOk)
        edge (i,:) = NaN;
    end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%   CLASS Meteo_Network
% =========================================================================
%
% DESCRIPTION
%   Class to store and manage a network of meteo station
%
% EXAMPLE
%   settings = Meteo_Network();
%
% FOR A LIST OF CONSTANTs and METHODS use doc Meteo_Data
%
% REFERENCE

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Andrea Gatti, Giulio Tagliaferro ...
%  Contributors:      Andrea Gatti, Giulio Tagliaferro ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%----------------------------------------------------------------------------------------------
classdef Meteo_Network < handle
    properties
        mds % list of meteo data
        log % handfle of logger
    end
    %% METHOD CREATOR
    % ==================================================================================================================================================
    methods (Static)
        % Concrete implementation.  See Singleton superclass.
        function this = Meteo_Network()
            this.log = Core.getLogger();
        end
    end

    %% METHODS UI
    % ==================================================================================================================================================
    methods (Static, Access = public)
        function ok_go = openGUI()
            ok_go = gui_goGPS;
        end
    end

    methods (Access = public)
        function initSession(this, date_start, date_stop)
            % Load all the metereological data
            %
            % SYNTAX
            % this.initSession(data_start, data_stop)

            this.mds = [];
            % load all meteo file present in current settings
            state = Core.getCurrentSettings();
            fnames = state.getMetFileName(date_start, date_stop);
            this.initFromFileList(fnames)
        end

        function initFromFileList(this, file_name_list)
            % Load metereological data directly from file list
            %
            % INPUT
            %   file_name_list: two level cells, first level station second level days
            %
            % SYNTAX
            %   this.initFromFileList(file_name_list)

            n_stations = numel(file_name_list);
            for  i = 1 : n_stations
                n_sss = numel(file_name_list{i});
                for  s = 1 : n_sss
                    if ~isempty(File_Name_Processor.getFileName(file_name_list{i}{s}))
                        if ~exist(file_name_list{i}{s}, 'file')
                            this.log.addWarning(sprintf('Skipping %s - file not found', file_name_list{i}{s}));
                        else
                            md = Meteo_Data(file_name_list{i}{s});
                            if md.isValid()
                                % md.setMaxBound(0);
                                if length(this.mds) < i
                                    this.mds = [this.mds; md];
                                else
                                    this.mds(i).inject(md);
                                end
                            end
                        end
                    end
                end
            end
        end

        function md = getVMS(this, name, xyz, time)
            % Get Virtual Meteo Station
            %
            % INPUT
            %   name    name of the new Meteo_Data virtual station
            %   xyz     coordinates of the new meteo station
            %   time    time of interpolation
            %
            %
            % OUTPUT
            %   md      virtual Meteo_Data station generated at xyz coordinates
            %
            % SYNTAX
            %   md = this.getVMS(name, xyz, time)
            %
            % EXAMPLE
            %   [x, y, z, amsl] = station(1).getLocation();
            %   md1 = this.getVMS('test', [x y z], station(1).getObsTime)
            md = Meteo_Data.getVMS( name, xyz, time, this.mds);
        end
    end

end

author: goGPS-Project


text:
function  message=read_text(file_name)
      % Copyright 2002-2010 The MathWorks, Inc.
      fid=fopen(file_name,'rt');
      message=char(fread(fid))';
      fclose(fid);

author: Sable

repo: mcbench-benchmarks


text:
function hpaData=parseHPA(fileName, version)
% parseHPA
%   Parses a database dump of the Human Protein Atlas (HPA)
%
%   Input:
%   fileName            comma- or tab-separated database dump of HPA. For details
%                       regarding the format, see
%                       http://www.proteinatlas.org/about/download.
%   version             version of HPA [optional, default=19]
%
%
%   Output:
%   hpaData
%       genes               cell array with the unique gene names. In
%                           version >=18 this is the ensemble name, see
%                           geneNames below for the names in ver >=18
%       geneNames           cell array with the gene names, indexed the
%                           same way as genes.
%       tissues             cell array with the tissue names. The list may not be
%                           unique, as there can be multiple cell types per tissue
%       celltypes           cell array with the cell type names for each tissue
%       levels              cell array with the unique expression levels
%       types               cell array with the unique evidence types
%       reliabilities       cell array with the unique reliability levels
%
%       gene2Level          gene-to-expression level mapping in sparse matrix form.
%                           The value for element i,j is the index in
%                           hpaData.levels of gene i in cell type j
%       gene2Type           gene-to-evidence type mapping in sparse matrix form.
%                           The value for element i,j is the index in
%                           hpaData.types of gene i in cell type j. Doesn't
%                           exist in version >=18.
%       gene2Reliability    gene-to-reliability level mapping in sparse matrix form.
%                           The value for element i,j is the index in
%                           hpaData.reliabilities of gene i in cell type j
%
%   Usage: hpaData=parseHPA(fileName,version)

if nargin<2
    version=19; %Change this and add code for more versions when the current HPA version is increased and the format is changed
end

fileName=char(fileName);
if ~isfile(fileName)
    error('HPA file %s cannot be found',string(fileName));
end

if (version >= 17)
    fid=fopen(fileName,'r');
    hpa=textscan(fid,'%q %q %q %q %q %q','Delimiter','\t');
    fclose(fid);

    %Go through and see if the headers match what was expected
    headers={'Gene' 'Gene name' 'Tissue' 'Cell type' 'Level' 'Reliability'};
    for i=1:numel(headers)
        if ~strcmpi(headers(i),hpa{i}(1))
            EM=['Could not find the header "' headers{i} '". Make sure that the input file matches the format specified at http://www.proteinatlas.org/about/download'];
            dispEM(EM);
        end
        %Remove the header line here
        hpa{i}(1)=[];
    end

    %Get the unique values of each data type
    [hpaData.genes, P, I]=unique(hpa{1});
    hpaData.geneNames=hpa{2}(P); %make this vector use the index as genes
    [~, J, K]=unique(strcat(hpa{3},'€',hpa{4}));
    hpaData.tissues=hpa{3}(J);
    hpaData.celltypes=hpa{4}(J);
    [hpaData.levels, ~, L]=unique(hpa{5});
    [hpaData.reliabilities, ~, N]=unique(hpa{6});

    %Map the data to be sparse matrises instead
    hpaData.gene2Level=sparse(I,K,L,numel(hpaData.genes),numel(hpaData.tissues));
    hpaData.gene2Reliability=sparse(I,K,N,numel(hpaData.genes),numel(hpaData.tissues));
else
    fid=fopen(fileName,'r');
    hpa=textscan(fid,'%q %q %q %q %q %q','Delimiter',',');
    fclose(fid);

    %Go through and see if the headers match what was expected
    headers={'Gene' 'Tissue' 'Cell type' 'Level' 'Expression type' 'Reliability'};
    for i=1:numel(headers)
        if ~strcmpi(headers(i),hpa{i}(1))
            EM=['Could not find the header "' headers{i} '". Make sure that the input file matches the format specified at http://www.proteinatlas.org/about/download'];
            dispEM(EM);
        end
        %Remove the header line here
        hpa{i}(1)=[];
    end

    %Get the unique values of each data type
    [hpaData.genes, ~, I]=unique(hpa{1});
    [~, J, K]=unique(strcat(hpa{2},'€',hpa{3}));
    hpaData.tissues=hpa{2}(J);
    hpaData.celltypes=hpa{3}(J);
    [hpaData.levels, ~, L]=unique(hpa{4});
    [hpaData.types, ~, M]=unique(hpa{5});
    [hpaData.reliabilities, ~, N]=unique(hpa{6});

    %Map the data to be sparse matrises instead
    hpaData.gene2Level=sparse(I,K,L,numel(hpaData.genes),numel(hpaData.tissues));
    hpaData.gene2Type=sparse(I,K,M,numel(hpaData.genes),numel(hpaData.tissues));
    hpaData.gene2Reliability=sparse(I,K,N,numel(hpaData.genes),numel(hpaData.tissues));
end
end

author: SysBioChalmers

repo: RAVEN


text:
function [model, directions] = directionalityStats(model, directions, cumNormProbCutoff, printLevel)
% Build Boolean vectors with reaction directionality statistics
%
% USAGE:
%
%    [model, directions] = directionalityStats(model, directions, cumNormProbCutoff, printLevel)
%
% INPUTS:
%    model:                 a structue of boolean vectors with different directionality
%                           assignments where some vectors contain subsets of others with field:
%
%                             * .directions
%
%    directions:            structure with field:
%
%                             * .forwardProbability
%
%                           qualitatively assigned internal reaction direactions
%
%                             * .forwardRecon
%                             * .reverseRecon
%                             * .reversibleRecon
%                             * .equilibriumRecon
%
%                           quantitatively assigned internal reaction direactions
%                           thermodynamic data is lacking
%
%                             * .forwardThermo
%                             * .reverseThermo
%                             * .reversibleThermo
%                             * .uncertainThermo
%                             * .equilibriumThermo
%
% OPTIONAL INPUT
%    cumNormProbCutoff:     {0.2} cutoff for probablity that reaction is
%                           reversible within this cutoff of 0.5
%    printLevel:            -1 to print out to file,
%                           0 to silent,
%                           1 to print out to command window
%
% OUTPUT:
%    directions:            a structue of boolean vectors with different directionality
%                           assignments where some vectors contain subsets of others
%
%                           qualtiative -> quantiative changed reaction directions
%
%                             * .forward2Forward
%                             * .forward2Reverse
%                             * .forward2Reversible
%                             * .forward2Uncertain
%                             * .reversible2Forward
%                             * .reversible2Reverse
%                             * .reversible2Reversible
%                             * .reversible2Uncertain
%                             * .reverse2Forward
%                             * .reverse2Reverse
%                             * .reverse2Reversible
%                             * .reverse2Uncertain
%                             * .tightened
%
%                           subsets of qualtiatively forward  -> quantiatively reversible
%
%                             * .forward2Reversible_bydGt0
%                             * .forward2Reversible_bydGt0LHS
%                             * .forward2Reversible_bydGt0Mid
%                             * .forward2Reversible_bydGt0RHS
%                             * .forward2Reversible_byConc_zero_fixed_DrG0
%                             * .forward2Reversible_byConc_negative_fixed_DrG0
%                             * .forward2Reversible_byConc_positive_fixed_DrG0
%                             * .forward2Reversible_byConc_negative_uncertain_DrG0
%                             * .forward2Reversible_byConc_positive_uncertain_DrG0
%
% .. Author: - Ronan M.T. Fleming

if ~exist('cumNormProbCutoff','var')
    directions.cumNormProbCutoff=0.2;
else
    directions.cumNormProbCutoff=cumNormProbCutoff;
end
%must be symmetric about 50:50 to be logically consistent
cumNormProbFwdUpper=0.5+directions.cumNormProbCutoff;
cumNormProbFwdLower=0.5-directions.cumNormProbCutoff;

if ~exist('printLevel','var')
    printLevel=0;
end
if ~exist('fileName','var')
    fileName='directionalityStats.txt';
end

DrGtMin=model.DrGtMin;
DrGtMax=model.DrGtMax;
if any(DrGtMin>DrGtMax)
    error('DrGtMin greater than DrGtMax');
end

DrGtNaNBool=(isnan(model.DrGtMax) | isnan(model.DrGtMin)) & model.SIntRxnBool;
if any(DrGtNaNBool)
    warning([int2str(nnz(DrGtNaNBool)) ' DrGt are NaN']);
end

nEqualDrGt=nnz(DrGtMin==DrGtMax & DrGtMin~=0);
if any(nEqualDrGt)
    fprintf('%s\n',[num2str(nEqualDrGt) '/' num2str(length(DrGtMin)) ' reactions with DrGtMin=DrGtMax~=0' ]);
end

nZeroDrGt=nnz(DrGtMin==0 & DrGtMax==0);
if any(nZeroDrGt)
    fprintf('%s\n',[num2str(nZeroDrGt) '/' num2str(length(DrGtMin)) ' reactions with DrGtMin=DrGtMax=0' ]);
end

[~,nRxn]=size(model.S);

% qualitatively assigned directions
forwardRecon=directions.forwardRecon;
reverseRecon=directions.reverseRecon;
reversibleRecon=directions.reversibleRecon;
% quantitatively assigned directions
forwardThermo=directions.forwardThermo;
reverseThermo=directions.reverseThermo;
reversibleThermo=directions.reversibleThermo;
uncertainThermo=directions.uncertainThermo;

%%%%%%CHANGES IN REACTION DIRECTIONS%%%%%%%%%%%
%thermodynamic constraints tightened
tightened=model.lb<model.lb_reconThermo & model.ub_reconThermo<model.ub;

reversible2Forward=reversibleRecon & forwardThermo;
reversible2Reverse=reversibleRecon & reverseThermo;
reversible2Reversible=reversibleRecon & reversibleThermo;
reversible2Uncertain=reversibleRecon & uncertainThermo;

forward2Reverse=forwardRecon & reverseThermo;
forward2Reversible=forwardRecon & reversibleThermo;
forward2Forward=forwardRecon & forwardThermo;
forward2Uncertain=forwardRecon & uncertainThermo;

reverse2Reverse =  reverseRecon & reverseThermo;
reverse2Forward   =  reverseRecon & forwardThermo;
reverse2Reversible =  reverseRecon & reversibleThermo;
reverse2Uncertain  =  reverseRecon & uncertainThermo;

%%%%%%CAUSES TO CHANGES IN REACTION DIRECTIONS%%%%%%%%%%%
% model.DrGtMax = model.DrGt0Max + gasConstant*T*(R'*log(model.concMax) - F'*log(model.concMin));
% model.DrGtMin = model.DrGt0Min + gasConstant*T*(R'*log(model.concMin) - F'*log(model.concMax));
% model.DrGtMaxMeanConc = model.DrGt0Max + gasConstant*T*(R-F)'*log((model.concMax+model.concMin)/2);
% model.DrGtMinMeanConc = model.DrGt0Min + gasConstant*T*(R-F)'*log((model.concMax+model.concMin)/2);

forward2Reversible_bydGt0=forwardRecon & reversibleThermo & model.DrGtMin<0 & model.DrGtMax>0; % dGfGCforward2ReversibleBool_bydGt0

forward2Reversible_byConc_negative_fixed_DrG0 = forwardRecon & reversibleThermo & model.DrGt0Min==model.DrGt0Max & model.DrGtMax<=0; %dGfGCforward2ReversibleBool_byConc_No_dGt0ErrorLHS
forward2Reversible_byConc_positive_fixed_DrG0 = forwardRecon & reversibleThermo & model.DrGt0Min==model.DrGt0Max & model.DrGtMin>0; %dGfGCforward2ReversibleBool_byConc_No_dGt0ErrorRHS

forward2Reversible_bydGt0LHS=forward2Reversible_bydGt0 & directions.forwardProbability>cumNormProbFwdUpper; %dGfGCforward2ReversibleBool_bydGt0LHS
forward2Reversible_bydGt0Mid=forward2Reversible_bydGt0 & directions.forwardProbability>=cumNormProbFwdLower & directions.forwardProbability<=cumNormProbFwdUpper;%dGfGCforward2ReversibleBool_bydGt0Mid
forward2Reversible_bydGt0RHS=forward2Reversible_bydGt0 & directions.forwardProbability<cumNormProbFwdLower; %dGfGCforward2ReversibleBool_bydGt0RHS

forward2Reversible_byConc_negative_uncertain_DrG0 = forwardRecon & reversibleThermo & model.DrGt0Min~=model.DrGt0Max & model.DrGt0Max<0; %dGfGCforward2ReversibleBool_byConcLHS
forward2Reversible_byConc_positive_uncertain_DrG0 = forwardRecon & reversibleThermo & model.DrGt0Min~=model.DrGt0Max & model.DrGt0Min>0; %dGfGCforward2ReversibleBool_byConcRHS

forward2Reversible_byConc_zero_fixed_DrG0 = forwardRecon & reversibleThermo & model.DrGt0Min==0 & model.DrGt0Max==0;%new to v2

if printLevel<0
    fid=fopen(fileName,'w');
else
    fid=1;
end

if printLevel~=0
    fprintf(fid,'%s\n','Qualitative internal reaction directionality:');
    fprintf(fid,'%10s\t%s\n',int2str(nnz(model.SIntRxnBool)),' internal reconstruction reaction directions.');
    fprintf(fid,'%10s\t%s\n',int2str(nnz(forwardRecon)), ' forward reconstruction assignment.');
    fprintf(fid,'%10s\t%s\n',int2str(nnz(reverseRecon)), ' reverse reconstruction assignment.');
    fprintf(fid,'%10s\t%s\n',int2str(nnz(reversibleRecon)), ' reversible reconstruction assignment.');
    fprintf(fid,'\n');

    fprintf(fid,'%s\n','Quantitative internal reaction directionality:');
    fprintf(fid,'%10s\t%s\n',int2str(nnz(model.SIntRxnBool)),' internal reconstruction reaction directions.');
    fprintf(fid,'%10s\t%s\n',int2str(nnz(forwardThermo)+nnz(reverseThermo)+nnz(reversibleThermo)),  ' of which have a thermodynamic assignment.');
    fprintf(fid,'%10s\t%s\n',int2str(nnz(uncertainThermo)),  ' of which have no thermodynamic assignment.');
    fprintf(fid,'%10s\t%s\n',int2str(nnz(forwardThermo)), ' forward thermodynamic only assignment.');
    fprintf(fid,'%10s\t%s\n',int2str(nnz(reverseThermo)), ' reverse thermodynamic only assignment.');
    fprintf(fid,'%10s\t%s\n',int2str(nnz(reversibleThermo)), ' reversible thermodynamic only assignment.');
    fprintf(fid,'\n');

    fprintf(fid,'%s\n','Qualitiative vs Quantitative:');
    fprintf(fid,'%10i\t%s\n',nnz(reversible2Reversible),' Reversible -> Reversible');
    fprintf(fid,'%10i\t%s\n',nnz(reversible2Forward),' Reversible -> Forward');
    fprintf(fid,'%10i\t%s\n',nnz(reversible2Reverse),' Reversible -> Reverse');
    fprintf(fid,'%10i\t%s\n',nnz(reversible2Uncertain),' Reversible -> Uncertain');
    fprintf(fid,'%10i\t%s\n',nnz(forward2Forward),' Forward -> Forward');
    fprintf(fid,'%10i\t%s\n',nnz(forward2Reverse),' Forward -> Reverse');
    fprintf(fid,'%10i\t%s\n',nnz(forward2Reversible),' Forward -> Reversible');
    fprintf(fid,'%10i\t%s\n',nnz(forward2Uncertain),' Forward -> Uncertain');
    fprintf(fid,'%10i\t%s\n',nnz(reverse2Reversible),' Reverse -> Reverse');
    fprintf(fid,'%10i\t%s\n',nnz(reverse2Forward),' Reverse -> Forward');
    fprintf(fid,'%10i\t%s\n',nnz(reverse2Reversible),' Reverse -> Reversible');
    fprintf(fid,'%10i\t%s\n',nnz(reverse2Uncertain),' Reversible -> Uncertain');
    fprintf(fid,'\n');

    fprintf(fid,'%s\n','Breakdown of relaxation of reaction directionality, Qualitiative vs Quantitative:');
    %total number of qualitatively forward reactions that are
    %quantitatively reversible
    fprintf(fid,'%10i\t%s\n',nnz(forward2Reversible),' qualitatively forward reactions that are quantitatively reversible (total).');
    %qualitatively forward reactions that are quantitatively reversible by
    %the range of dGt0
    fprintf(fid,'%10i\t%s%s\n',nnz(forward2Reversible_bydGt0LHS),' of which are quantitatively reversible by range of dGt0. ',['P(\Delta_{r}G^{\primeo}<0) > ' num2str(cumNormProbFwdUpper)]);
    %qualitatively reverse reactions that are quantitatively
    %reversible by concentration alone (with dGt0 error)
    fprintf(fid,'%10i\t%s%s\n',nnz(forward2Reversible_bydGt0Mid),' of which are quantitatively reversible by range of dGt0. ',[num2str(cumNormProbFwdLower) '< P(\Delta_{r}G^{\primeo}<0) < ' num2str(cumNormProbFwdUpper)]);
    %qualitatively reverse reactions that are quantitatively
    %reversible by concentration alone (with dGt0 error)
    fprintf(fid,'%10i\t%s%s\n',nnz(forward2Reversible_bydGt0RHS),' of which are quantitatively reversible by range of dGt0. ',['P(\Delta_{r}G^{\primeo}<0) < ' num2str(cumNormProbFwdLower)]);
    %qualitatively forward reactions that are quantitatively
    %reversible by concentration alone (no dGt0 error)
    fprintf(fid,'%10i\t%s\n',nnz(forward2Reversible_byConc_zero_fixed_DrG0),' of which are quantitatively forward by fixed dGr0t, but reversible by concentration alone (zero fixed DrGt0).');
    %qualitatively reverse reactions that are quantitatively
    %reversible by concentration alone (with dGt0 error)
    fprintf(fid,'%10i\t%s\n',nnz(forward2Reversible_byConc_negative_fixed_DrG0),' of which are quantitatively reverse by dGr0t, but reversible by concentration (negative fixed DrGt0).');
    %qualitatively reverse reactions that are quantitatively
    %reversible by concentration alone (with dGt0 error)
    fprintf(fid,'%10i\t%s\n',nnz(forward2Reversible_byConc_positive_fixed_DrG0),' of which are quantitatively forward by dGr0t, but reversible by concentration (positve fixed DrGt0).');
    %qualitatively reverse reactions that are quantitatively
    %reversible by concentration alone (with dGt0 error)
    fprintf(fid,'%10i\t%s\n',nnz(forward2Reversible_byConc_negative_uncertain_DrG0),' of which are quantitatively reverse by dGr0t, but reversible by concentration (uncertain negative DrGt0).');
    %qualitatively reverse reactions that are quantitatively
    %reversible by concentration alone (with dGt0 error)
    fprintf(fid,'%10i\t%s\n',nnz(forward2Reversible_byConc_positive_uncertain_DrG0),' of which are quantitatively forward by dGr0t, but reversible by concentration (uncertain positive DrGt0).');
end
if printLevel<0
    fclose(fid);
end

%changed directions
directions.forward2Forward=forward2Forward;
directions.forward2Reverse=forward2Reverse;
directions.forward2Reversible=forward2Reversible;
directions.forward2Uncertain=forward2Uncertain;
directions.reversible2Forward=reversible2Forward;
directions.reversible2Reverse=reversible2Reverse;
directions.reversible2Reversible=reversible2Reversible;
directions.reversible2Uncertain=reversible2Uncertain;
directions.reverse2Forward=reverse2Forward;
directions.reverse2Reverse=reverse2Reverse;
directions.reverse2Reversible=reverse2Reversible;
directions.reverse2Uncertain=reverse2Uncertain;

directions.tightened=tightened;

%all forward reversible classes
directions.forward2Reversible_bydGt0=forward2Reversible_bydGt0;
directions.forward2Reversible_bydGt0LHS=forward2Reversible_bydGt0LHS;
directions.forward2Reversible_bydGt0Mid=forward2Reversible_bydGt0Mid;
directions.forward2Reversible_bydGt0RHS=forward2Reversible_bydGt0RHS;

directions.forward2Reversible_byConc_zero_fixed_DrG0=forward2Reversible_byConc_zero_fixed_DrG0;
directions.forward2Reversible_byConc_negative_fixed_DrG0=forward2Reversible_byConc_negative_fixed_DrG0;
directions.forward2Reversible_byConc_positive_fixed_DrG0=forward2Reversible_byConc_positive_fixed_DrG0;
directions.forward2Reversible_byConc_negative_uncertain_DrG0=forward2Reversible_byConc_negative_uncertain_DrG0;
directions.forward2Reversible_byConc_positive_uncertain_DrG0=forward2Reversible_byConc_positive_uncertain_DrG0;

author: opencobra

repo: cobratoolbox


text:
function x = r8mat_utsol ( n, r, b )

%*****************************************************************************80
%
%% R8MAT_UTSOL solves R' * X = B for an upper triangular matrix R.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2013
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real R(N,N), the upper triangular matrix.
%
%    Input, real B(N), the right hand side.
%
%    Output, real X(N), the solution.
%
  x(1:n,1) = b(1:n);

  x(1,1) = x(1,1) / r(1,1);

  for j = 2 : n
    x(j,1) = ( x(j,1) - r(1:j-1,j)' * x(1:j-1,1) ) / r(j,j);
  end

  return
end


text:
function [sol, infos] = gsp_wiener_optimization(G, x0, f, psd, psd_noise, param)
%GSP_WIENER_OPTIMIZATION Solve wiener optimization problem
%   Usage:  sol = gsp_wiener_optimization(G, x0, ffid, psd, psd_noise)
%           sol = gsp_wiener_optimization(G, x0, ffid, psd, psd_noise, param)
%           [sol, infos] = gsp_wiener_optimization(...)
%
%   Input parameters:
%         G          : Graph (GSP structure)
%         x0         : Starting point (column vector)
%         f          : Fidelity term - UNLocBox structure
%         psd        : PSD filter (anonymous function)
%         psd_noise  : PSD filter of the noise or single number
%         param      : Optional optimization parameters
%   Output parameters:
%         sol        : Solution
%         infos      : Convergence informations
%
%   This function solves the following wiener optimization problem:
%
%     .. argmin_x f(x) + || w(L) x ||_2^2
%
%     .. math:: arg\min_x f(x) + \| w(L) x \|_2^2
%
%   Please refer to the reference for more information about this problem.
%   This function requires the UNLocBox to work.
%
%   Please refer to the function gsp_filter_analysis and solvep to know how
%   *param* can be set.
%
%   References: perraudin2016stationary

% Author : Nathanael Perraudin
% Date: 6 January 2016

if nargin<6
    param = struct;
end

if isnumeric(psd_noise)
    if sum(abs(psd_noise(:)))==0
        error('This function cannot solve this problem')
    end
    if sum(abs(psd_noise(:)))<1e-10
        warning('This function can prabaly not solve this case');
    end
    wl = @(x) psd_noise./(psd(x)+eps);
    fprox = @(T) @(x) psd(x)./(psd(x)+2*T*psd_noise + eps);
else
    wl = @(x) psd_noise(x)./(psd(x)+eps);
    fprox = @(T) @(x) psd(x)./(psd(x)+2*T*psd_noise(x) + eps);

end

%fprox = @(x) 1./(wl(x)+1);

% In order to be faster
param.stopping_criterion = 'rel_norm_obj';

% Wiener term
fwiener.prox = @(x,T) gsp_filter_analysis(G,fprox(T),x, param);
fwiener.eval = @(x) 0.5*norm(gsp_filter_analysis(G,wl,x,param),'fro')^2;

% Call the solver
[sol , infos ] = solvep(x0,{f,fwiener},param);

end

repo: gspbox


text:
function nodeSet = coloring(G,nc,qp)
%% COLORING color vertices of a graph
%
% nodeSet = coloring(G,nc) color the vertices of the graph G by nc colors.
% Vertices in the same color are disjointed except the last one. The input
% graph G is represented by the sparse pattern of a symmetric sparse
% matrix. The values of G is not used. The output nodeSet is a cell
% array and nodeSet{k} is the indices of vertices in the k-th color. By
% default, nc = 6 if not specified by the input.
%
% nodeSet = coloring(G,nc,qp) chose the vertices with a larger value qp
% first. For example, qp can be the degree of vertices. Then nodeSet{1}
% will pick up nodes with a local max degree. If no qp, random values are
% used.
%
% Example
%
%     load lakemesh
%     N = size(node,1); NT = size(elem,1);
%     t2p = sparse([1:NT,1:NT,1:NT], elem(1:NT,:), 1, NT, N);
%     G = t2p'*t2p;
%     nodeSet = coloring(G);
%     showcoloring(G,node,nodeSet);
%
% See also: showcoloring, showgraphmatrix
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Set up
N = size(G,1);
N0 = min(floor(N/100),20);  % number of the coarest nodes
if ~exist('nc','var'), nc = 8; end
nodeSet = cell(nc,1);
if ~exist('qp','var')
    qp = rand(N,1);
else
    idx = (qp ~= 0);
    qp(idx) = (1 + 0.1*rand(sum(idx),1)).*qp(idx);
end

%% Coloring
for k = 1:nc-1
    isC = false(N,1);       % C: coarse nodes
    isF = false(N,1);       % F: fine nodes
    isU = (qp>0);           % U: undecided nodes
    if ~any(isU)            % no more vertices left
        break
    end
    for m = 1:round(log10(N))
        % Mark all available nodes
        isS = isU;          % selected all undecided vertices
        S = find(isS);

        % Find marked nodes with local maximum quality
        [i,j] = find(triu(G(S,S),1));   % i,j and i<j: edges of subgraph S
        idx = (qp(S(i)) >= qp(S(j)));   % compare quality of connected vertices
        isS(S(j(idx))) = false;         % remove vertices with smaller qp
        isS(S(i(~idx))) = false;        % if qualities are equal, keep the nodes with smaller index
        isC(isS) = true;                % set selected nodes as coarse nodes

        % Remove coarse nodes and neighboring nodes from undecided set
        [i,j] = find(G(:,isC)); %#ok<*NASGU>
        isF(i) = true;                  % neighbor of C nodes are F nodes
        isU(isF | isC) = false;         % remove current C and F from U
        if sum(isU) <= N0
            break;    % only break the inner loop
        end
    end
    nodeSet{k} = find(isC);
    qp(nodeSet{k}) = 0;
end
nodeSet{k+1} = find(qp>0);  % all left vertices are added to the last color
% nodeSet = nodeSet{1:k+1};

repo: ifem


text:
%function geyser ( )

%*****************************************************************************80
%
%% GEYSER uses MATLAB to make a histogram.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 April 2011
%
%  Author:
%
%    John Burkardt
%

%
%  Load the data.
%
  t = load ( 'geyser.txt' );
%
%  Our data lies in the range 40 to 110.
%  Dividing into bins of width 5, we want 14 bins.
%  But that means 14+1 "edges".
%
  bin_num = 14;
  bin_edges = linspace ( 40.0, 110.0, bin_num + 1 );
%
%  T_BINNED counts data between successive bin edges.
%
  t_binned = histc ( t, bin_edges );
%
%  Discard the last entry, which counts data greater than the last bin edge.
%
  t_binned = t_binned(1:bin_num);
%
%  Set the average value of each bin.
%
  bar_center = linspace ( ( bin_edges(1)     + bin_edges(2)    ) / 2, ...
                          ( bin_edges(end-1) + bin_edges(end) ) / 2, ...
                          bin_num );
%
%  Set the relative width of the bar.
%
  bar_width = 0.90;
%
%  Draw a bar graph.
%
  bar ( bar_center, t_binned, bar_width, 'r' );

  grid on
  title ( 'Time between eruptions of Old Faithful' )
  xlabel ( 'Minutes' )
  ylabel ( 'Frequency' )

%  return
%end


text:
function h = figcopy(fc,fm)
% FIGCOPY - copies one figure into another figure and lets the
% users move and  resize the copied figure. This utility is great
% if you need a figure in a figure and still would like to use
% MATLAB's export function.
%
% h = FIGCOPY the user selects both figures by mouse clicks. Please
% follow the instructions in the command window. First activate the
% figure to be copied. Second activate the "master" figure and
% select the lower left and upper right corner, where the copied
% figure should appear in the "master" figure. Click on the copied
% figure to move or resize it. (Don't worry, the resize markers
% will not be printed!) h is the handle of the copied figure.
%
% h = FIGCOPY(fc) fc specifies the figure number (handle) to be
% copied. The master figure is selected by mouse click.
%
% h = FIGCOPY(fc,fm) specifies the figure number to be copied and
% the figure number of the master figure.

% Dirk Tenne
% CoDE (Control Dynamics and Estimation)
% web:      "http://code.eng.buffalo.edu"
% created:  "long time ago"
% modified: "December 12. 2002" updated the input handling
%           "November 27. 2007" copy only figure axes ignoring the new
%           uicontrols of Figures.
%

% input handling
  % is the figure number given or should it be clicked?
  if nargin < 3
    select_fc = 0;
    select_fm = 0;
    if nargin < 2
      select_fc = 0;
      select_fm = 1;
      if nargin < 1
	select_fc = 1;
	select_fm = 1;
      end
    end
  end

%activate the figure to copy

  if select_fc
    disp('Click on the figure which you''d like to copy!')
    waitforbuttonpress;
    fc = gcf;
    disp(' ')
    disp(['Figure ',num2str(gcf),' copied into the clipboard.'])
  end

  %hcfig = get(fc,'Children');
  % only look for axes handles as new versions also return uicontrols
  hcfig = findobj(1,'Type','axes');

  if select_fm
    disp('Activate master figure! (Try clicking the frame.)')
    waitforbuttonpress;
    fm = gcf;
  end

% specify the position of the figure
  figure(fm);
  disp(['Click to indicate the LOWER LEFT and UPPER RIGHT'])
  disp(['corners, where the copied figure should appear!'])
  [xi,yi] = ginput(2);

%normalizing x and y
   m_axis = axis;
   m_ax_pos = get(gca,'position');
   alphax = m_ax_pos(3); % rescaling the units
   alphay = m_ax_pos(4);
   x_max = diff(m_axis(1:2)); y_max = diff(m_axis(3:4));
   x = (xi-m_axis(1))/x_max*alphax;
   y = (yi-m_axis(3))/y_max*alphay;

   % placing the figure
   pos(1) = m_ax_pos(1) + x(1);
   pos(2) = m_ax_pos(2) + y(1);
   siz(1) = abs(x(1)-x(2));
   siz(2) = abs(y(1)-y(2));
   newhcfig = copyobj(hcfig,fm);
   for j = 1:length(newhcfig)
     set(newhcfig(j),'Position',[pos siz])
   end
   set(newhcfig(1),'ButtonDownFcn','selectmoveresize')
   %disable:set(gco,'ButtonDownFcn','','Selected','off')

% output handling
   if nargout > 0
     h = newhcfig;
   end


text:
classdef PerformanceTests < handle & matlab.unittest.TestCase

    properties (TestParameter)

    end

    methods (Test, TestTags = {'Performance'})

        function test2D(testCase)
            index = 1;
            s.dim    = '2D';
            s.length    = 1;
            s.height = 0.1;
            steps = 0.0001:0.0001:0.1;
            index = 1;
            tic
            % tic should start AFTer assembling the connec
            % maybe store it?
            temps = zeros(size(steps));
            nelems = zeros(size(steps));
%             for step = steps
            step = 0.0001;
                test = PerformanceTest(s);
                sol = test.compute(step);
                nelem = sol.getMesh().nelem;
                disp(nelem)
                temps(index) = toc;
                nelems(index) = nelem;
                index = index+1;
%             end
            disp(temps)
            disp(nelems)
            plot(nelems, temps)
        end

%         function test3D(testCase)
%             index = 1;
%             for i = 0.2:0.1:0.5
%                 gg = @() testCase.example3D(i);
%                 temps(index) = timeit(gg);
%                 connecs(index) = size(testCase.example3D(i),1);
%                 index = index+1;
%             end
%         end

    end

    methods (Access = private)

        function connec = example2D(testCase, step)
            s.dim    = '2D';
            s.length    = 1;
            s.height = 0.1;
            test = PerformanceTest(s);
            sol = test.compute(step);
            connec = test.getConnec();
        end

        function connec = example3D(testCase, step)
            s.dim    = '3D';
            s.length    = 1;
            s.height = 0.1;
            test = PerformanceTest(s);
            sol = test.compute(step);
            connec = test.getConnec();
        end

    end

end


text:
function result = merge_results(num_subjects, solver_name)
    filepath = sprintf('bin/yale_faces_%d_subjects_test_%s.mat', num_subjects, solver_name);
    load(filepath);
    num_trials = numel(trials);
    fprintf('Number of subjects: %d, Number of trials: %d, solver: %s\n', num_subjects,  num_trials, solver_name);
    result = struct;
    trial = trials{1};
    result.M = trial.M;
    result.K = trial.K;
    result.D = trial.D;
    result.elapsed_time_vec = zeros(num_trials, 1);
    result.connectivity_vec = zeros(num_trials, 1);
    result.clustering_error_perc_vec = zeros(num_trials, 1);
    result.clustering_acc_perc_vec = zeros(num_trials, 1);
    result.spr_error_vec = zeros(num_trials, 1);
    result.spr_flag_vec = zeros(num_trials, 1);
    result.spr_perc_vec = zeros(num_trials, 1);
    num_points_vec = zeros(num_trials, 1);
    result.trials = trials;
    for i=1:numel(trials)
        trial = trials{i};
        result.elapsed_time_vec(i) = trial.elapsed_time;
        result.connectivity_vec(i) = trial.connectivity;
        result.clustering_error_perc_vec(i) = trial.clustering_error_perc;
        result.clustering_acc_perc_vec(i) = trial.clustering_acc_perc;
        result.spr_error_vec(i) = trial.spr_error;
        result.spr_flag_vec(i) = trial.spr_flag;
        result.spr_perc_vec(i) = trial.spr_perc;
        num_points_vec(i) = trial.S;
    end
    result.num_points = mean(num_points_vec);
    result.elapsed_time = mean(result.elapsed_time_vec);
    result.connectivity = min(result.connectivity_vec);
    result.clustering_error_perc = mean(result.clustering_error_perc_vec);
    result.clustering_acc_perc = mean(result.clustering_acc_perc_vec);
    result.spr_error = mean(result.spr_error_vec);
    result.spr_flag = mean(result.spr_flag_vec);
    result.spr_perc = mean(result.spr_perc_vec);

    fprintf('\n\n');
    fprintf('Mean Points: ');
    fprintf('%.2f ', result.num_points);
    fprintf('\n\n');
    fprintf('Time: ');
    fprintf('%.2f ', result.elapsed_time);
    fprintf('\n\n');
    fprintf('Connectivity: ');
    fprintf('%.2f ', result.connectivity);
    fprintf('\n\n');
    fprintf('Clustering error(%%): ');
    fprintf('%.2f ', result.clustering_error_perc);
    fprintf('\n\n');
    fprintf('Clustering accuracy (%%): ');
    fprintf('%.2f ', result.clustering_acc_perc);
    fprintf('\n\n');
    fprintf('Subspace preserving representation error: ');
    fprintf('%.4f ', result.spr_error);
    fprintf('\n\n');
    fprintf('Subspace preserving representation flag: ');
    fprintf('%.2f ', result.spr_flag);
    fprintf('\n\n');
    fprintf('Subspace preserving representations(%%): ');
    fprintf('%.2f ', result.spr_perc);
    fprintf('\n\n');
end

author: indigits

repo: sparse-plex


text:
% Demonstrates the projection of x^2 function
clear all; close all; clc;

% The function to be projected.
func = @(x) (x-2).^2;

% epsilon
polarity = 0;
step = 0.01;
epsilon = 0.2;
stop = 2;
x = 0:step:stop;
func_x = func(x);
projection = spx.wavelet.lcs.project_0_right(func, step, stop, epsilon, polarity);
subplot(311);
plot(x, func_x, 'r');
subplot(312);
plot(x, projection, 'b');
residual = func_x - projection;
subplot(313);
plot(x, residual, 'g');


text:
function I = mi_model_gd_vec(x, y, Ym, biascorrect, demeaned)
% MI_MODEL_GD_VEC Vectorized MI calculation between multiple Gaussian variables
%         and a common discrete variable in bits based on ANOVA style model comparison.
%   I = mi_model_gd_vec(x,y,Ym) returns the MI between the (possibly multidimensional)
%   Gaussian variables x and the discrete variable y.
%   size(x) = [Ntrl Nvec Ndim]
%   so each output I(i) = mi_model_gd(squeeze(x(:,i,:)), y, Ym);
%
%   For 1D x this is a lower bound to the mutual information.
%   y should contain integer values in the range [0 Ym-1] (inclusive).
%
%   biascorrect : true / false option (default true) which specifies whether
%   bias correction should be applied to the esimtated MI.
%   demeaned : false / true option (default false) which specifies whether the
%   input data already has zero mean (true if it has been copula-normalized)
%   See also: MI_MODEL_GD, MI_MIXTURE_GD_VEC

% ensure samples first axis for vectors
if isvector(x)
    x = x(:);
end
if ndims(x)>3
    error('mi_model_gd: input arrays should be 3d')
end
if isvector(y)
    y = y(:);
else
    error('mi_model_gd: only univariate discrete variable supported');
end

Ntrl = size(x,1);
Nvec = size(x,2);
Xdim = size(x,3);

if size(y,1) ~= Ntrl
    error('mi_model_gd: number of trials do not match');
end

% default option values
if nargin<4
    biascorrect = true;
end

if nargin<5
    demeaned = false;
end

% unconditional demean
if ~demeaned
    x = bsxfun(@minus,x,sum(x,1)/Ntrl);
end

I = zeros(Nvec,1);

% y = y-1;
% for vi=1:Nvec
%   I(vi) = mi_model_gd(squeeze(x(:,vi,:)),y,Ym,true,true);
% end
%
% return

% one-hot encoding of Y
Yhot = indexed2boolean(y);

% remove class means
[Xcen class_means] = removeclassmeans(x, Yhot);

% allocate memory for class-conditional entropies and covariances
Ntrl_y = sum(Yhot);
Hcond  = zeros(Nvec,Ym);
Cm     = zeros(Nvec,Xdim,Xdim,Ym);

% allocate memory for overall entropy and covariance
Hunc = zeros(Nvec,1);
Cx   = zeros(Nvec,Xdim,Xdim);

%  data is class-demeaned, this needs to be accounted for in the
% unconditional entropies
c = diag(sqrt(Ntrl_y))*class_means.';
c = reshape(c, [Ym Nvec Xdim]);

for vi1=1:Xdim
  % all voxels for this dimension
  x1 = Xcen(:,:,vi1);
  c1 = squeeze(c(:,:,vi1));

  Cx(:,vi1,vi1) = sum(x1.^2)+sum(c1.^2);
  for yi=1:Ym
    tmp = x1(Yhot(:,yi),:);
    Cm(:,vi1,vi1,yi) = sum(tmp.^2);
  end

  for vi2=(vi1+1):Xdim
    x2  = Xcen(:,:,vi2);
    c2  = squeeze(c(:,:,vi2));
    tmp = transpose(sum(x1.*x2) + sum(c1.*c2));

    Cx(:,vi1,vi2) = tmp;
    Cx(:,vi2,vi1) = tmp;
    for yi=1:Ym
      tmp = transpose(sum(x1(Yhot(:,yi),:).*x2(Yhot(:,yi),:)));
      Cm(:,vi1,vi2,yi) = tmp;
      Cm(:,vi2,vi1,yi) = tmp;
    end
  end
end

Cx = Cx / (Ntrl-1);
Cx = vecchol(Cx);
for vi=1:Xdim
  Hunc = Hunc + shiftdim(log(Cx(:,vi,vi)));
end

for yi=1:Ym
  Cm(:,:,:,yi) = Cm(:,:,:,yi) / (Ntrl_y(yi) - 1);
  Cm(:,:,:,yi) = vecchol(Cm(:,:,:,yi));
  for vi=1:Xdim
    Hcond(:,yi) = Hcond(:,yi)+shiftdim(log(Cm(:,vi,vi,yi)));
  end
end

% apply bias corrections
ln2 = log(2);
if biascorrect

  vars = 1:Xdim;

  psiterms_unc = psi((Ntrl - vars)/2) / 2;
  dterm_unc    = (ln2 - log(Ntrl-1)) / 2;
  bias_unc     = Xdim'*dterm_unc + sum(psiterms_unc);

  dterm_cond    = (ln2 - log(Ntrl_y-1)) / 2;
  psiterms_cond = zeros(1,Ym);
  for vi=vars
    idx = (Ntrl_y-vi);
    psiterms_cond = psiterms_cond + psi(idx/2);
  end
  bias_cond = Xdim*dterm_cond + (psiterms_cond/2);

  Hunc  = Hunc  - bias_unc;
  Hcond = Hcond - ones(Nvec,1)*bias_cond;
end

% class weights
w = Ntrl_y ./ Ntrl;

% compute mutual information
I = Hunc - Hcond*w';

% convert to bits
I = I / ln2;

function [Xcen, class_means] = removeclassmeans(X, design)
[Ntrl, Nvec, Ndim] = size(X);
Xcen = X(:,:);
class_means = zeros(Nvec*Ndim,size(design,2));
for k = 1:size(design,2)
  sel = design(:,k);
  tmp = Xcen(sel,:);
  class_means(:,k) = mean(tmp,1);
  Xcen(sel,:) = bsxfun(@minus,tmp,class_means(:,k).');
end
Xcen = reshape(Xcen,[Ntrl Nvec Ndim]);

function Y = indexed2boolean(X)
uX = unique(X);
Y  = false(numel(X),numel(uX));
for k = 1:size(Y,2)
  Y(X==uX(k),k) = true;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function convBpool()
    global config mem;
    curr_layer_idx = config.misc.current_layer;
    mem.deltas{curr_layer_idx+1} = reshape(accumarray(mem.convBpool{curr_layer_idx+2}{1}, mem.deltas{curr_layer_idx+1}(:)), ...
                                           mem.orig_activation_size{curr_layer_idx+1}(1), mem.orig_activation_size{curr_layer_idx+1}(2));
                                            %size(mem.layer_inputs{curr_layer_idx+1}, 1), size(mem.layer_inputs{curr_layer_idx+1}, 2));
    mem.grads{curr_layer_idx+1} = sum(mem.layer_inputs{curr_layer_idx+1} .* mem.deltas{curr_layer_idx+1}, 1)';
    mem.grads{curr_layer_idx+1+config.layer_num} = sum(mem.deltas{curr_layer_idx+1}, 1)';
	mem.delta_act = bsxfun(@times, mem.deltas{curr_layer_idx+1}(mem.convBpool{curr_layer_idx+2}{2})', config.weights{curr_layer_idx+1});
end

author: jimmy-ren


text:
%  Figure 3.30e      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%

clf;
num=1;
a=.5;

zeta =1;
den2=[1/(zeta*a) 1];
den1=[1 2*zeta 1];
den=conv(den1,den2);
t=0:.1:10;
y1=step(num,den,t);

zeta =.7;
den2=[1/(zeta*a) 1];
den1=[1 2*zeta 1];
den=conv(den1,den2);
y2=step(num,den,t);

zeta =.5;
den2=[1/(zeta*a) 1];
den1=[1 2*zeta 1];
den=conv(den1,den2);
y3=step(num,den,t);

axis([1 10 .1 .9])
plot(t,y1,'-',t,y2,'--',t,y3,'-.'),grid
title('Fig. 3.30e Step response with extra pole, \alpha= .5')
xlabel('\omega_n t')
ylabel('y(t)')


text:
function [select] = parameterselection(param, data);

% PARAMETERSELECTION selects the parameters that are present as a volume in the data
% add that have a dimension that is compatible with the specified dimensions of the
% volume, i.e. either as a vector or as a 3D volume.
%
% Use as
%   [select] = parameterselection(param, data)
% where
%   param    cell-array, or single string, can be 'all'
%   data     structure with anatomical or functional data
%   select   returns the selected parameters as a cell-array

% Copyright (C) 2005-2008, 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: parameterselection.m 2885 2011-02-16 09:41:58Z roboos $

if ischar(param)
  param = {param};   % it should be a cell-array
elseif isempty(param)
  param = {};        % even being empty, it should be a cell-array
end

sel = find(strcmp(param, 'all'));
if ~isempty(sel)
  % the old default was a list of all known volume parameters
  % the new default is to try all fields present in the data
  allparam = fieldnames(data);
  % fields can be nested in source.avg
  if isfield(data, 'avg')
    tmp = fieldnames(data.avg);
    for i=1:length(tmp)
      tmp{i} = ['avg.' tmp{i}];
    end
    allparam = cat(1, allparam, tmp);
  end
  % fields can be nested in source.trial
  if isfield(data, 'trial')
    tmp = fieldnames(data.trial);
    for i=1:length(tmp)
      tmp{i} = ['trial.' tmp{i}];
    end
    allparam = cat(1, allparam, tmp);
  end
  param(sel) = [];                          % remove the 'all'
  param      = [param(:)' allparam(:)'];    % add the list of all possible parameters, these will be tested later
else
  % check all specified parameters and give support for some parameters like 'pow' and 'coh'
  % which most often will indicate 'avg.pow' and 'avg.coh'
  for i=1:length(param)
    if ~issubfield(data, param{i}) && issubfield(data, ['avg.' param{i}])
      % replace the parameter xxx by avg.xxx
      param{i} = ['avg.' param{i}];
    end
  end
end

% remove empty fields
param(find(cellfun('isempty', param))) = [];

% ensure that there are no double entries
param = unique(param);

select = {};
for i=1:length(param)
  if issubfield(data, param{i})
    % the field is present, check whether the dimension is correct
    dim = size(getsubfield(data, param{i}));
    if isfield(data, 'dim') && isequal(dim(:), data.dim(:))
      select{end+1} = param{i};
    elseif isfield(data, 'dim') && prod(dim)==prod(data.dim)
      select{end+1} = param{i};
    elseif isfield(data, 'pos') && prod(dim)==size(data.pos, 1)
      select{end+1} = param{i};
    elseif isfield(data, 'dimord') && (isfield(data, 'pos') || isfield(data, 'transform')),
      dimtok = tokenize(data.dimord, '_');
      nels   = 1;
      for k=1:numel(dimtok)
        if strcmp(dimtok{k}, 'rpt') || strcmp(dimtok{k}, 'rpttap')
	  nels = nels*dim(k);
	elseif strcmp(dimtok{k}, 'pos') && isfield(data, 'pos')
	  %for source structure
	  nels = nels*size(data.pos,1);
	elseif strcmp(dimtok{k}, 'pos')
	  %for volume structure FIXME 'pos' in dimord is not OK
	  nels = nels*prod(data.dim(1:3));
	else
          nels = nels*numel(getfield(data, dimtok{k}));
	end
      end
      if nels==prod(dim),
        select{end+1} = param{i};
      end
    end
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function tests = test_spm_cfg_dcm_fmri
% Unit Tests for spm_cfg_dcm_fmri (DCM fMRI spec batch)
%__________________________________________________________________________
% Copyright (C) 2016 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_cfg_dcm_fmri.m 7479 2018-11-09 14:17:33Z peter $

tests = functiontests(localfunctions);

% -------------------------------------------------------------------------
function setup(testCase)
% Delete artefacts before each test

% Delete GCM from the tmp directory
outpath = fullfile(get_data_path(), 'tmp');
if exist(outpath,'file')
    delete(fullfile(outpath,'*.mat'));
end

% Delete DCMs from the GLM directory
glmdir = fullfile(get_data_path(), 'GLM');
dcms   = cellstr(spm_select('FPListRec',glmdir,'^DCM_.*mat$'));
for i = 1:length(dcms)
    spm_unlink(dcms{i});
end

% -------------------------------------------------------------------------
function test_specify_group(testCase)

data_path = get_data_path();

% Output path & directory
outpath = fullfile(data_path, 'tmp');
dir_out = outpath;
if ~exist(dir_out,'file')
    mkdir(dir_out);
end

% Template DCMs
fullm   = fullfile(data_path,'models','DCM_s1_m1.mat');
altm    = {fullfile(data_path,'models','DCM_s1_m2.mat');
           fullfile(data_path,'models','DCM_s1_m3.mat');
           fullfile(data_path,'models','DCM_s1_m4.mat')};

% All subjects' SPMs and VOIs
spms    = cellstr(spm_select('FPListRec',data_path,'SPM.mat'));
voi1    = cellstr(spm_select('FPListRec',data_path,'VOI_R1.mat'));
voi2    = cellstr(spm_select('FPListRec',data_path,'VOI_R2.mat'));

% Run
matlabbatch{1}.spm.dcm.spec.fmri.group.output.dir  = cellstr(dir_out);
matlabbatch{1}.spm.dcm.spec.fmri.group.output.name = 'specify_group_test';
matlabbatch{1}.spm.dcm.spec.fmri.group.template.fulldcm = cellstr(fullm);
matlabbatch{1}.spm.dcm.spec.fmri.group.template.altdcm  = altm;
matlabbatch{1}.spm.dcm.spec.fmri.group.data.spmmats      = spms;
matlabbatch{1}.spm.dcm.spec.fmri.group.data.session      = 1;
matlabbatch{1}.spm.dcm.spec.fmri.group.data.region = {voi1
                                                      voi2}';

out = spm_jobman('run',matlabbatch);

% Load output
GCM = load(out{1}.gcmmat{1});
GCM = GCM.GCM;
GCM = spm_dcm_load(GCM);

% Check each subject's DCM timeseries match their VOIs
for s = 1:size(GCM,1)

    % Load VOI files
    xY1 = load(voi1{s});
    xY2 = load(voi2{s});
    Y_expected = [xY1.Y xY2.Y];

    % Check VOI match actual timeseries
    for m = 1:size(GCM,2)
        Y_actual   = GCM{s,m}.Y.y;
        testCase.assertTrue(all(Y_expected(:)==Y_actual(:)), ...
            sprintf('Timeseries for subject %d model %m wrong',s,m));
    end
end

% -------------------------------------------------------------------------
function data_path = get_data_path()

data_path = fullfile( spm('Dir'), 'tests', ...
    'data', 'fMRI', 'simulated_2region');

author: spm


text:
%SlideDistanceLimit

%w = sqrt(2*g / 3*L);
%cos(th) = 2/3;
%sin(th) = sqrt(1-cos(th)^2);
%
%rg = [L*sin(th); L*cos(th)];
%
% Let L=1, g = 1
%

costh = 2/3;
sinth = sqrt(1-costh^2);
w = sqrt(2/3);

rg = [-sinth; costh];
drg = [-w*costh; -w*sinth];

%0 = rg(2) + drg(2)*t + -0.5*t^2
P = [-0.5,drg(2),rg(2)];
t = max(roots(P));

d = rg(1) + drg(1)*t;

SlipDist = d+1;

%%%%%%%%%%%%%%

syms L g

L=sym(1);
g = sym(1);

costh = sym(2/3);
sinth = sqrt(sym(5/9));
w = sqrt((2*g)/(3*L));

rg = [-L*sinth; L*costh];
drg = -L*w*[costh; sinth];
a = -g/2;
b = drg(2);
c = rg(2);
t = simplify((-b - sqrt(b^2-4*a*c))/(2*a));

d = simplify(rg(1) + drg(1)*t);

SlipDist = simplify(d + L);

pretty(SlipDist)
SlipDist
double(SlipDist)

author: MatthewPeterKelly

repo: dscTutorials


text:
function Set = intersect_elements(Set1,Set2,False1,False2)

% Determines the intersection of Set1 and Set2.

Set = unique_elements([Set1; Set2],False1);
False1(Set1) = true;
False2(Set2) = true;
I = False1(Set)&False2(Set);
Set = Set(I);

author: InverseTampere

repo: TreeQSM


text:
function [out] = infiltration_8(S,Smax,fin)
%infiltration_8

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Infiltration into storage is equal the inflow when current
%               storage is under the maximum storage,
%               and zero when storage reaches maximum capacity
% Constraints:  f <= fin
% @(Inputs):
%               S    - current storage [mm]
%               Smax - maximum storage [mm]
%               fin  - size of incoming flux [mm/d]

out = (S < Smax) .* fin;

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT