text:
classdef TwoBodyImpactPointTime < AbstractConstraint
    %TwoBodyImpactPointTime Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function obj = TwoBodyImpactPointTime(event, lb, ub)
            obj.event = event;
            obj.lb = lb;
            obj.ub = ub;

            obj.id = rand();
        end

        function [lb, ub] = getBounds(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, celBodyData)
            type = obj.getConstraintType();

            switch obj.eventNode
                case ConstraintStateComparisonNodeEnum.FinalState
                    stateLogEntry = stateLog.getLastStateLogForEvent(obj.event);

                case ConstraintStateComparisonNodeEnum.InitialState
                    stateLogEntry = stateLog.getFirstStateLogForEvent(obj.event);

                otherwise
                    error('Unknown event node.');
            end

            value = lvd_TwoBodyImpactPointTasks(stateLogEntry, 'timeToImpact');

            if(isnan(value))
                value = 1E99; %need to throw in a real number
            end

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                switch obj.stateCompNode
                    case ConstraintStateComparisonNodeEnum.FinalState
                        stateLogEntryStateComp = stateLog.getLastStateLogForEvent(obj.stateCompEvent).deepCopy();

                    case ConstraintStateComparisonNodeEnum.InitialState
                        stateLogEntryStateComp = stateLog.getFirstStateLogForEvent(obj.stateCompEvent).deepCopy();

                    otherwise
                        error('Unknown event node.');
                end

                cartElem = stateLogEntryStateComp.getCartesianElementSetRepresentation();
                cartElem = cartElem.convertToFrame(stateLogEntry.centralBody.getBodyCenteredInertialFrame());
                stateLogEntryStateComp.setCartesianElementSet(cartElem);

                valueStateComp = lvd_TwoBodyImpactPointTasks(stateLogEntryStateComp, 'timeToImpact');

                if(isnan(valueStateComp))
                    valueStateComp = 1E99; %need to throw in a real number
                end
            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;
        end

        function tf = canUseSparseOutput(obj)
            tf = true;
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function type = getConstraintType(obj)
            type = 'Two-Body Time To Impact';
        end

%         function name = getName(obj)
%             name = sprintf('%s - Event %i', obj.getConstraintType(), obj.event.getEventNum());
%         end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = 'sec';
            lbLim = -Inf;
            ubLim = Inf;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
%             addConstraintTf = lvd_EditGenericMAConstraintGUI(obj, lvdData);

            output = AppDesignerGUIOutput({false});
            lvd_EditGenericMAConstraintGUI_App(obj, lvdData, output);
            addConstraintTf = output.output{1};
        end
    end

    methods(Static)
        function constraint = getDefaultConstraint(~, ~)
            constraint = TwoBodyImpactPointTime(LaunchVehicleEvent.empty(1,0),0,0);
        end
    end
end

author: Arrowstar

repo: ksptot


text:
%
% output = my_grep(input, pattern)
%
% Similate the function of grep in Unix
% The inputs are
%       input: texts in cell format
%       pattern: the pattern that is searched for
% The output is
%       output: lines in input that contains the pattern
%
% Author: Xiong Xiao,
%         School of Computer Engineering, Nanyang Technological University
% Date Created: 23 Oct, 2008
% Last Modified: 23 Oct, 2008
%
function output = my_grep(input, pattern)

N = length(input);

output{1} = '';
cnt = 1;
for i=1:N
    if length(regexp(input{i},pattern)) > 0
        output{cnt} = input{i};
        cnt = cnt + 1;
    end
end

author: singaxiong

repo: SignalGraph


text:
function varargout = exportfig(varargin)
%EXPORTFIG  Export a figure.
%   EXPORTFIG(H, FILENAME) writes the figure H to FILENAME.  H is
%   a figure handle and FILENAME is a string that specifies the
%   name of the output file.
%
%   EXPORTFIG(H, FILENAME, OPTIONS) writes the figure H to FILENAME
%   with options initially specified by the structure OPTIONS. The
%   field names of OPTIONS must be legal parameters listed below
%   and the field values must be legal values for the corresponding
%   parameter. Default options can be set in releases prior to R12
%   by storing the OPTIONS structure in the root object's appdata
%   with the command
%      setappdata(0,'exportfigdefaults', OPTIONS)
%   and for releases after R12 by setting the preference with the
%   command
%      setpref('exportfig', 'defaults', OPTIONS)
%
%   EXPORTFIG(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies
%   parameters that control various characteristics of the output
%   file. Any parameter value can be the string 'auto' which means
%   the parameter uses the default factory behavior, overriding
%   any other default for the parameter.
%
%   Format Paramter:
%     'Format'  a string
%          specifies the output format. Defaults to 'eps'. For a
%          list of export formats type 'help print'.
%     'Preview' one of the strings 'none', 'tiff'
%          specifies a preview for EPS files. Defaults to 'none'.
%
%   Size Parameters:
%     'Width'   a positive scalar
%          specifies the width in the figure's PaperUnits
%     'Height'  a positive scalar
%          specifies the height in the figure's PaperUnits
%     'Bounds' one of the strings 'tight', 'loose'
%          specifies a tight or loose bounding box. Defaults to 'tight'.
%     'BoundsCode' one of the strings 'internal', 'mcode'
%          specifies whether to use the internal EPS code or
%          general m-code for computing bounding boxes. Defaults to
%          'internal'.
%     'Reference' an axes handle or a string
%          specifies that the width and height parameters
%          are relative to the given axes. If a string is
%          specified then it must evaluate to an axes handle.
%
%     Specifying only one dimension sets the other dimension
%     so that the exported aspect ratio is the same as the
%     figure's or reference axes' current aspect ratio.
%     If neither dimension is specified the size defaults to
%     the width and height from the figure's or reference
%     axes' size. Tight bounding boxes are only computed for
%     2-D views and in that case the computed bounds enclose all
%     text objects.
%
%   Rendering Parameters:
%     'Color'     one of the strings 'bw', 'gray', 'cmyk'
%         'bw'    specifies that lines and text are exported in
%                 black and all other objects in grayscale
%         'gray'  specifies that all objects are exported in grayscale
%         'rgb'   specifies that all objects are exported in color
%                 using the RGB color space
%         'cmyk'  specifies that all objects are exported in color
%                 using the CMYK color space
%     'Renderer'  one of 'painters', 'zbuffer', 'opengl'
%         specifies the renderer to use
%     'Resolution'   a positive scalar
%         specifies the resolution in dots-per-inch.
%     'LockAxes'  one of 0 or 1
%         specifies that all axes limits and ticks should be fixed
%         while exporting.
%
%     The default color setting is 'bw'.
%
%   Font Parameters:
%     'FontMode'     one of the strings 'scaled', 'fixed'
%     'FontSize'     a positive scalar
%          in 'scaled' mode multiplies with the font size of each
%          text object to obtain the exported font size
%          in 'fixed' mode specifies the font size of all text
%          objects in points
%     'DefaultFixedFontSize' a positive scalar
%          in 'fixed' mode specified the default font size in
%          points
%     'FontSizeMin' a positive scalar
%          specifies the minimum font size allowed after scaling
%     'FontSizeMax' a positive scalar
%          specifies the maximum font size allowed after scaling
%     'FontEncoding' one of the strings 'latin1', 'adobe'
%          specifies the character encoding of the font
%     'SeparateText' one of 0 or 1
%          specifies that the text objects are stored in separate
%          file as EPS with the base filename having '_t' appended.
%
%     If FontMode is 'scaled' but FontSize is not specified then a
%     scaling factor is computed from the ratio of the size of the
%     exported figure to the size of the actual figure.
%
%     The default 'FontMode' setting is 'scaled'.
%
%   Line Width Parameters:
%     'LineMode'     one of the strings 'scaled', 'fixed'
%     'LineWidth'    a positive scalar
%     'DefaultFixedLineWidth' a positive scalar
%     'LineWidthMin' a positive scalar
%          specifies the minimum line width allowed after scaling
%     'LineWidthMax' a positive scalar
%          specifies the maximum line width allowed after scaling
%     The semantics of 'Line' parameters are exactly the
%     same as the corresponding 'Font' parameters, except that
%     they apply to line widths instead of font sizes.
%
%   Style Map Parameter:
%     'LineStyleMap'    one of [], 'bw', or a function name or handle
%          specifies how to map line colors to styles. An empty
%          style map means styles are not changed. The style map
%          'bw' is a built-in mapping that maps lines with the same
%          color to the same style and otherwise cycles through the
%          available styles. A user-specified map is a function
%          that takes as input a cell array of line objects and
%          outputs a cell array of line style strings. The default
%          map is [].
%
%   Examples:
%     exportfig(gcf,'fig1.eps','height',3);
%       Exports the current figure to the file named 'fig1.eps' with
%       a height of 3 inches (assuming the figure's PaperUnits is
%       inches) and an aspect ratio the same as the figure's aspect
%       ratio on screen.
%
%     opts = struct('FontMode','fixed','FontSize',10,'height',3);
%     exportfig(gcf, 'fig2.eps', opts, 'height', 5);
%       Exports the current figure to 'fig2.eps' with all
%       text in 10 point fonts and with height 5 inches.
%
%   See also PREVIEWFIG, APPLYTOFIG, RESTOREFIG, PRINT.

%  Copyright 2000-2009 The MathWorks, Inc.

if (nargin < 2)
  error('Too few input arguments');
end

% exportfig(H, filename, [options,] ...)
H = varargin{1};
if ~LocalIsHG(H,'figure')
  error('First argument must be a handle to a figure.');
end
filename = varargin{2};
if ~ischar(filename)
  error('Second argument must be a string.');
end
paramPairs = {varargin{3:end}};
if nargin > 2
  if isstruct(paramPairs{1})
    pcell = LocalToCell(paramPairs{1});
    paramPairs = {pcell{:}, paramPairs{2:end}};
  end
end
verstr = version;
majorver = str2num(verstr(1));
defaults = [];
if majorver > 5
  if ispref('exportfig','defaults')
    defaults = getpref('exportfig','defaults');
  end
elseif exist('getappdata')
  defaults = getappdata(0,'exportfigdefaults');
end
if ~isempty(defaults)
  dcell = LocalToCell(defaults);
  paramPairs = {dcell{:}, paramPairs{:}};
end

% Do some validity checking on param-value pairs
if (rem(length(paramPairs),2) ~= 0)
  error(['Invalid input syntax. Optional parameters and values' ...
	 ' must be in pairs.']);
end

auto.format = 'eps';
auto.preview = 'none';
auto.width = -1;
auto.height = -1;
auto.color = 'bw';
auto.defaultfontsize=10;
auto.fontsize = -1;
auto.fontmode='scaled';
auto.fontmin = 8;
auto.fontmax = 60;
auto.defaultlinewidth = 1.0;
auto.linewidth = -1;
auto.linemode=[];
auto.linemin = 0.5;
auto.linemax = 100;
auto.fontencoding = 'latin1';
auto.renderer = [];
auto.resolution = [];
auto.stylemap = [];
auto.applystyle = 0;
auto.refobj = -1;
auto.bounds = 'tight';
auto.boundscode = 'internal';
explicitbounds = 0;
auto.lockaxes = 1;
auto.separatetext = 0;
opts = auto;

% Process param-value pairs
args = {};
for k = 1:2:length(paramPairs)
  param = lower(paramPairs{k});
  if ~ischar(param)
    error('Optional parameter names must be strings');
  end
  value = paramPairs{k+1};

  switch (param)
   case 'format'
    opts.format = LocalCheckAuto(lower(value),auto.format);
    if strcmp(opts.format,'preview')
      error(['Format ''preview'' no longer supported. Use PREVIEWFIG' ...
	     ' instead.']);
    end
   case 'preview'
    opts.preview = LocalCheckAuto(lower(value),auto.preview);
    if ~strcmp(opts.preview,{'none','tiff'})
      error('Preview must be ''none'' or ''tiff''.');
    end
   case 'width'
    opts.width = LocalToNum(value, auto.width);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.width)
	error('Width must be a numeric scalar > 0');
      end
    end
   case 'height'
    opts.height = LocalToNum(value, auto.height);
    if ~ischar(value) | ~strcmp(value,'auto')
      if(~LocalIsPositiveScalar(opts.height))
	error('Height must be a numeric scalar > 0');
      end
    end
   case 'color'
    opts.color = LocalCheckAuto(lower(value),auto.color);
    if ~strcmp(opts.color,{'bw','gray','rgb','cmyk'})
      error('Color must be ''bw'', ''gray'',''rgb'' or ''cmyk''.');
    end
   case 'fontmode'
    opts.fontmode = LocalCheckAuto(lower(value),auto.fontmode);
    if ~strcmp(opts.fontmode,{'scaled','fixed'})
      error('FontMode must be ''scaled'' or ''fixed''.');
    end
   case 'fontsize'
    opts.fontsize = LocalToNum(value,auto.fontsize);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.fontsize)
	error('FontSize must be a numeric scalar > 0');
      end
    end
   case 'defaultfixedfontsize'
    opts.defaultfontsize = LocalToNum(value,auto.defaultfontsize);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.defaultfontsize)
	error('DefaultFixedFontSize must be a numeric scalar > 0');
      end
    end
   case 'fontsizemin'
    opts.fontmin = LocalToNum(value,auto.fontmin);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.fontmin)
	error('FontSizeMin must be a numeric scalar > 0');
      end
    end
   case 'fontsizemax'
    opts.fontmax = LocalToNum(value,auto.fontmax);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.fontmax)
	error('FontSizeMax must be a numeric scalar > 0');
      end
    end
   case 'fontencoding'
    opts.fontencoding = LocalCheckAuto(lower(value),auto.fontencoding);
    if ~strcmp(opts.fontencoding,{'latin1','adobe'})
      error('FontEncoding must be ''latin1'' or ''adobe''.');
    end
   case 'linemode'
    opts.linemode = LocalCheckAuto(lower(value),auto.linemode);
    if ~strcmp(opts.linemode,{'scaled','fixed'})
      error('LineMode must be ''scaled'' or ''fixed''.');
    end
   case 'linewidth'
    opts.linewidth = LocalToNum(value,auto.linewidth);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.linewidth)
	error('LineWidth must be a numeric scalar > 0');
      end
    end
   case 'defaultfixedlinewidth'
    opts.defaultlinewidth = LocalToNum(value,auto.defaultlinewidth);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.defaultlinewidth)
	error(['DefaultFixedLineWidth must be a numeric scalar >' ...
	       ' 0']);
      end
    end
   case 'linewidthmin'
    opts.linemin = LocalToNum(value,auto.linemin);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.linemin)
	error('LineWidthMin must be a numeric scalar > 0');
      end
    end
   case 'linewidthmax'
    opts.linemax = LocalToNum(value,auto.linemax);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.linemax)
	error('LineWidthMax must be a numeric scalar > 0');
      end
    end
   case 'linestylemap'
    opts.stylemap = LocalCheckAuto(value,auto.stylemap);
   case 'renderer'
    opts.renderer = LocalCheckAuto(lower(value),auto.renderer);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~strcmp(opts.renderer,{'painters','zbuffer','opengl'})
	error(['Renderer must be ''painters'', ''zbuffer'' or' ...
	       ' ''opengl''.']);
      end
    end
   case 'resolution'
    opts.resolution = LocalToNum(value,auto.resolution);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~(isnumeric(value) & (prod(size(value)) == 1) & (value >= 0));
	error('Resolution must be a numeric scalar >= 0');
      end
    end
   case 'applystyle' % means to apply the options and not export
    opts.applystyle = 1;
   case 'reference'
    if ischar(value)
      if strcmp(value,'auto')
	opts.refobj = auto.refobj;
      else
	opts.refobj = eval(value);
      end
    else
      opts.refobj = value;
    end
    if ~LocalIsHG(opts.refobj,'axes')
      error('Reference object must evaluate to an axes handle.');
    end
   case 'bounds'
    opts.bounds = LocalCheckAuto(lower(value),auto.bounds);
    explicitbounds = 1;
    if ~strcmp(opts.bounds,{'tight','loose'})
      error('Bounds must be ''tight'' or ''loose''.');
    end
   case 'boundscode'
    opts.boundscode = LocalCheckAuto(lower(value),auto.boundscode);
    if ~strcmp(opts.boundscode,{'internal','mcode'})
      error('BoundsCode must be ''internal'' or ''mcode''.');
    end
   case 'lockaxes'
    opts.lockaxes = LocalToNum(value,auto.lockaxes);
   case 'separatetext'
    opts.separatetext = LocalToNum(value,auto.separatetext);
   otherwise
    error(['Unrecognized option ' param '.']);
  end
end

% make sure figure is up-to-date
drawnow;

allLines  = findall(H, 'type', 'line');
allText   = findall(H, 'type', 'text');
allAxes   = findall(H, 'type', 'axes');
allImages = findall(H, 'type', 'image');
allLights = findall(H, 'type', 'light');
allPatch  = findall(H, 'type', 'patch');
allSurf   = findall(H, 'type', 'surface');
allRect   = findall(H, 'type', 'rectangle');
allFont   = [allText; allAxes];
allColor  = [allLines; allText; allAxes; allLights];
allMarker = [allLines; allPatch; allSurf];
allEdge   = [allPatch; allSurf];
allCData  = [allImages; allSurf];

% need to explicitly get ticks for some invisible figures
get(allAxes,'yticklabel');

old.objs = {};
old.prop = {};
old.values = {};

% Process format
if strncmp(opts.format,'eps',3) & ~strcmp(opts.preview,'none')
  args = {args{:}, ['-' opts.preview]};
end

hadError = 0;
oldwarn = warning;
try

  % lock axes limits, ticks and labels if requested
  if opts.lockaxes
    old = LocalManualAxesMode(old, allAxes, 'TickMode');
    old = LocalManualAxesMode(old, allAxes, 'TickLabelMode');
    old = LocalManualAxesMode(old, allAxes, 'LimMode');
  end

  % Process size parameters
  figurePaperUnits = get(H, 'PaperUnits');
  oldFigureUnits = get(H, 'Units');
  oldFigPos = get(H,'Position');
  set(H, 'Units', figurePaperUnits);
  figPos = get(H,'Position');
  refsize = figPos(3:4);
  if opts.refobj ~= -1
    oldUnits = get(opts.refobj, 'Units');
    set(opts.refobj, 'Units', figurePaperUnits);
    r = get(opts.refobj, 'Position');
    refsize = r(3:4);
    set(opts.refobj, 'Units', oldUnits);
  end
  aspectRatio = refsize(1)/refsize(2);
  if (opts.width == -1) & (opts.height == -1)
    opts.width = refsize(1);
    opts.height = refsize(2);
  elseif (opts.width == -1)
    opts.width = opts.height * aspectRatio;
  elseif (opts.height == -1)
    opts.height = opts.width / aspectRatio;
  end
  wscale = opts.width/refsize(1);
  hscale = opts.height/refsize(2);
  sizescale = min(wscale,hscale);
  old = LocalPushOldData(old,H,'PaperPositionMode', ...
			 get(H,'PaperPositionMode'));
  set(H, 'PaperPositionMode', 'auto');
  newPos = [figPos(1) figPos(2)+figPos(4)*(1-hscale) ...
	    wscale*figPos(3) hscale*figPos(4)];
  set(H, 'Position', newPos);
  set(H, 'Units', oldFigureUnits);

  old = LocalPushOldData(old,H,'Color', ...
			 get(H,'Color'));
  set(H,'Color','w');

  % process line-style map
  if ~isempty(opts.stylemap) & ~isempty(allLines)
    oldlstyle = LocalGetAsCell(allLines,'LineStyle');
    old = LocalPushOldData(old, allLines, {'LineStyle'}, ...
			   oldlstyle);
    newlstyle = oldlstyle;
    if ischar(opts.stylemap) & strcmpi(opts.stylemap,'bw')
      newlstyle = LocalMapColorToStyle(allLines);
    else
      try
	newlstyle = feval(opts.stylemap,allLines);
      catch
	warning(['Skipping stylemap. ' lasterr]);
      end
    end
    set(allLines,{'LineStyle'},newlstyle);
  end

  % Process rendering parameters
  switch (opts.color)
   case {'bw', 'gray'}
    if ~strcmp(opts.color,'bw') & strncmp(opts.format,'eps',3)
      opts.format = [opts.format 'c'];
    end
    args = {args{:}, ['-d' opts.format]};

    %compute and set gray colormap
    oldcmap = get(H,'Colormap');
    newgrays = 0.30*oldcmap(:,1) + 0.59*oldcmap(:,2) + 0.11*oldcmap(:,3);
    newcmap = [newgrays newgrays newgrays];
    old = LocalPushOldData(old, H, 'Colormap', oldcmap);
    set(H, 'Colormap', newcmap);

    %compute and set ColorSpec and CData properties
    old = LocalUpdateColors(allColor, 'color', old);
    old = LocalUpdateColors(allAxes, 'xcolor', old);
    old = LocalUpdateColors(allAxes, 'ycolor', old);
    old = LocalUpdateColors(allAxes, 'zcolor', old);
    old = LocalUpdateColors(allMarker, 'MarkerEdgeColor', old);
    old = LocalUpdateColors(allMarker, 'MarkerFaceColor', old);
    old = LocalUpdateColors(allEdge, 'EdgeColor', old);
    old = LocalUpdateColors(allEdge, 'FaceColor', old);
    old = LocalUpdateColors(allCData, 'CData', old);

    if strcmp(opts.color,'bw')
      lcolor = LocalGetAsCell(allLines,'color');
      n = length(lcolor);
      for k=1:n
	if (lcolor{k}(1) < 1-eps)
	  set(allLines(k),'color','k');
	end
      end
    end
   case {'rgb','cmyk'}
    if strncmp(opts.format,'eps',3)
      opts.format = [opts.format 'c'];
      args = {args{:}, ['-d' opts.format]};
      if strcmp(opts.color,'cmyk')
	args = {args{:}, '-cmyk'};
      end
    else
      args = {args{:}, ['-d' opts.format]};
    end
   otherwise
    error('Invalid Color parameter');
  end
  if (~isempty(opts.renderer))
    args = {args{:}, ['-' opts.renderer]};
  end
  if (~isempty(opts.resolution)) | ~strncmp(opts.format,'eps',3)
    if isempty(opts.resolution)
      opts.resolution = 0;
    end
    args = {args{:}, ['-r' int2str(opts.resolution)]};
  end

  % Process font parameters
  if ~isempty(opts.fontmode)
    oldfonts = LocalGetAsCell(allFont,'FontSize');
    oldfontunits = LocalGetAsCell(allFont,'FontUnits');
    set(allFont,'FontUnits','points');
    switch (opts.fontmode)
     case 'fixed'
      if (opts.fontsize == -1)
	set(allFont,'FontSize',opts.defaultfontsize);
      else
	set(allFont,'FontSize',opts.fontsize);
      end
     case 'scaled'
      if (opts.fontsize == -1)
	scale = sizescale;
      else
	scale = opts.fontsize;
      end
      newfonts = LocalScale(oldfonts,scale,opts.fontmin,opts.fontmax);
      set(allFont,{'FontSize'},newfonts);
     otherwise
      error('Invalid FontMode parameter');
    end
    old = LocalPushOldData(old, allFont, {'FontSize'}, oldfonts);
    old = LocalPushOldData(old, allFont, {'FontUnits'}, oldfontunits);
  end
  if strcmp(opts.fontencoding,'adobe') & strncmp(opts.format,'eps',3)
    args = {args{:}, '-adobecset'};
  end

  % Process line parameters
  if ~isempty(opts.linemode)
    oldlines = LocalGetAsCell(allMarker,'LineWidth');
    old = LocalPushOldData(old, allMarker, {'LineWidth'}, oldlines);
    switch (opts.linemode)
     case 'fixed'
      if (opts.linewidth == -1)
	set(allMarker,'LineWidth',opts.defaultlinewidth);
      else
	set(allMarker,'LineWidth',opts.linewidth);
      end
     case 'scaled'
      if (opts.linewidth == -1)
	scale = sizescale;
      else
	scale = opts.linewidth;
      end
      newlines = LocalScale(oldlines, scale, opts.linemin, opts.linemax);
      set(allMarker,{'LineWidth'},newlines);
    end
  end

  % adjust figure bounds to surround axes
  if strcmp(opts.bounds,'tight')
    if (~strncmp(opts.format,'eps',3) & LocalHas3DPlot(allAxes)) | ...
	  (strncmp(opts.format,'eps',3) & opts.separatetext)
      if (explicitbounds == 1)
	warning(['Cannot compute ''tight'' bounds. Using ''loose''' ...
		 ' bounds.']);
      end
      opts.bounds = 'loose';
    end
  end
  warning('off');
  if ~isempty(allAxes)
    usemcode = 1;
    if strncmp(opts.format,'eps',3)
      if strcmp(opts.boundscode,'internal')
	if strcmp(opts.bounds,'loose')
	  args = {args{:}, '-loose'};
	end
	old = LocalPushOldData(old,H,'Position', oldFigPos);
	usemcode = 0;
      else
	usemcode = 1;
      end
    end
    if strcmp(opts.bounds,'tight') & usemcode
      oldaunits = LocalGetAsCell(allAxes,'Units');
      oldapos = LocalGetAsCell(allAxes,'Position');
      oldtunits = LocalGetAsCell(allText,'units');
      oldtpos = LocalGetAsCell(allText,'Position');
      set(allAxes,'units','points');
      apos = LocalGetAsCell(allAxes,'Position');
      oldunits = get(H,'Units');
      set(H,'units','points');
      origfr = get(H,'position');
      fr = [];
      for k=1:length(allAxes)
	if ~strcmpi(get(allAxes(k),'Tag'),'legend')
	  axesR = apos{k};
	  r = LocalAxesTightBoundingBox(axesR, allAxes(k));
	  if ~isempty(r)
	    r(1:2) = r(1:2) + axesR(1:2);
	    fr = LocalUnionRect(fr,r);
	  end
	end
      end
      if isempty(fr)
	fr = [0 0 origfr(3:4)];
      end
      for k=1:length(allAxes)
	ax = allAxes(k);
	r = apos{k};
	r(1:2) = r(1:2) - fr(1:2);
	set(ax,'Position',r);
      end
      old = LocalPushOldData(old, allAxes, {'Position'}, oldapos);
      old = LocalPushOldData(old, allText, {'Position'}, oldtpos);
      old = LocalPushOldData(old, allText, {'Units'}, oldtunits);
      old = LocalPushOldData(old, allAxes, {'Units'}, oldaunits);
      old = LocalPushOldData(old, H, 'Position', oldFigPos);
      old = LocalPushOldData(old, H, 'Units', oldFigureUnits);
      r = [origfr(1) origfr(2)+origfr(4)-fr(4) fr(3:4)];
      set(H,'Position',r);
    end
  end

  % Process text in a separate file if needed
  if opts.separatetext & ~opts.applystyle
    % First hide all text and export
    oldtvis = LocalGetAsCell(allText,'visible');
    set(allText,'visible','off');
    oldax = LocalGetAsCell(allAxes,'XTickLabel',1);
    olday = LocalGetAsCell(allAxes,'YTickLabel',1);
    oldaz = LocalGetAsCell(allAxes,'ZTickLabel',1);
    null = cell(length(oldax),1);
    [null{:}] = deal([]);
    set(allAxes,{'XTickLabel'},null);
    set(allAxes,{'YTickLabel'},null);
    set(allAxes,{'ZTickLabel'},null);
    print(H, filename, args{:});
    set(allText,{'Visible'},oldtvis);
    set(allAxes,{'XTickLabel'},oldax);
    set(allAxes,{'YTickLabel'},olday);
    set(allAxes,{'ZTickLabel'},oldaz);
    % Now hide all non-text and export as eps in painters
    [path, name, ext] = fileparts(filename);
    tfile = fullfile(path,[name '_t.eps']);
    tfile2 = fullfile(path,[name '_t2.eps']);
    foundRenderer = 0;
    for k=1:length(args)
      if strncmp('-d',args{k},2)
	args{k} = '-deps';
      elseif strncmp('-zbuffer',args{k},8) | ...
	    strncmp('-opengl', args{k},6)
	args{k} = '-painters';
	foundRenderer = 1;
      end
    end
    if ~foundRenderer
      args = {args{:}, '-painters'};
    end
    allNonText = [allLines; allLights; allPatch; ...
		  allImages; allSurf; allRect];
    oldvis = LocalGetAsCell(allNonText,'visible');
    oldc = LocalGetAsCell(allAxes,'color');
    oldaxg = LocalGetAsCell(allAxes,'XGrid');
    oldayg = LocalGetAsCell(allAxes,'YGrid');
    oldazg = LocalGetAsCell(allAxes,'ZGrid');
    [null{:}] = deal('off');
    set(allAxes,{'XGrid'},null);
    set(allAxes,{'YGrid'},null);
    set(allAxes,{'ZGrid'},null);
    set(allNonText,'Visible','off');
    set(allAxes,'Color','none');
    print(H, tfile2, args{:});
    set(allNonText,{'Visible'},oldvis);
    set(allAxes,{'Color'},oldc);
    set(allAxes,{'XGrid'},oldaxg);
    set(allAxes,{'YGrid'},oldayg);
    set(allAxes,{'ZGrid'},oldazg);
    %hack up the postscript file
    fid1 = fopen(tfile,'w');
    fid2 = fopen(tfile2,'r');
    line = fgetl(fid2);
    while ischar(line)
      if strncmp(line,'%%Title',7)
	fprintf(fid1,'%s\n',['%%Title: ', tfile]);
      elseif (length(line) < 3)
	fprintf(fid1,'%s\n',line);
      elseif ~strcmp(line(end-2:end),' PR') & ...
	    ~strcmp(line(end-1:end),' L')
	fprintf(fid1,'%s\n',line);
      end
      line = fgetl(fid2);
    end
    fclose(fid1);
    fclose(fid2);
    delete(tfile2);

  elseif ~opts.applystyle
    drawnow;
    print(H, filename, args{:});
  end
  warning(oldwarn);

catch
  warning(oldwarn);
  hadError = 1;
end

% Restore figure settings
if opts.applystyle
  varargout{1} = old;
else
  for n=1:length(old.objs)
    if ~iscell(old.values{n}) & iscell(old.prop{n})
      old.values{n} = {old.values{n}};
    end
    set(old.objs{n}, old.prop{n}, old.values{n});
  end
end

if hadError
  error(deblank(lasterr));
end

%
%  Local Functions
%

function outData = LocalPushOldData(inData, objs, prop, values)
outData.objs = {objs, inData.objs{:}};
outData.prop = {prop, inData.prop{:}};
outData.values = {values, inData.values{:}};

function cellArray = LocalGetAsCell(fig,prop,allowemptycell);
cellArray = get(fig,prop);
if nargin < 3
  allowemptycell = 0;
end
if ~iscell(cellArray) & (allowemptycell | ~isempty(cellArray))
  cellArray = {cellArray};
end

function newArray = LocalScale(inArray, scale, minv, maxv)
n = length(inArray);
newArray = cell(n,1);
for k=1:n
  newArray{k} = min(maxv,max(minv,scale*inArray{k}(1)));
end

function gray = LocalMapToGray1(color)
gray = color;
if ischar(color)
  switch color(1)
   case 'y'
    color = [1 1 0];
   case 'm'
    color = [1 0 1];
   case 'c'
    color = [0 1 1];
   case 'r'
    color = [1 0 0];
   case 'g'
    color = [0 1 0];
   case 'b'
    color = [0 0 1];
   case 'w'
    color = [1 1 1];
   case 'k'
    color = [0 0 0];
  end
end
if ~ischar(color)
  gray = 0.30*color(1) + 0.59*color(2) + 0.11*color(3);
end

function newArray = LocalMapToGray(inArray);
n = length(inArray);
newArray = cell(n,1);
for k=1:n
  color = inArray{k};
  if ~isempty(color)
    color = LocalMapToGray1(color);
  end
  if isempty(color) | ischar(color)
    newArray{k} = color;
  else
    newArray{k} = [color color color];
  end
end

function newArray = LocalMapColorToStyle(inArray);
inArray = LocalGetAsCell(inArray,'Color');
n = length(inArray);
newArray = cell(n,1);
styles = {'-','--',':','-.'};
uniques = [];
nstyles = length(styles);
for k=1:n
  gray = LocalMapToGray1(inArray{k});
  if isempty(gray) | ischar(gray) | gray < .05
    newArray{k} = '-';
  else
    if ~isempty(uniques) & any(gray == uniques)
      ind = find(gray==uniques);
    else
      uniques = [uniques gray];
      ind = length(uniques);
    end
    newArray{k} = styles{mod(ind-1,nstyles)+1};
  end
end

function newArray = LocalMapCData(inArray);
n = length(inArray);
newArray = cell(n,1);
for k=1:n
  color = inArray{k};
  if (ndims(color) == 3) & isa(color,'double')
    gray = 0.30*color(:,:,1) + 0.59*color(:,:,2) + 0.11*color(:,:,3);
    color(:,:,1) = gray;
    color(:,:,2) = gray;
    color(:,:,3) = gray;
  end
  newArray{k} = color;
end

function outData = LocalUpdateColors(inArray, prop, inData)
value = LocalGetAsCell(inArray,prop);
outData.objs = {inData.objs{:}, inArray};
outData.prop = {inData.prop{:}, {prop}};
outData.values = {inData.values{:}, value};
if (~isempty(value))
  if strcmp(prop,'CData')
    value = LocalMapCData(value);
  else
    value = LocalMapToGray(value);
  end
  set(inArray,{prop},value);
end

function bool = LocalIsPositiveScalar(value)
bool = isnumeric(value) & ...
       prod(size(value)) == 1 & ...
       value > 0;

function value = LocalToNum(value,auto)
if ischar(value)
  if strcmp(value,'auto')
    value = auto;
  else
    value = str2num(value);
  end
end

%convert a struct to {field1,val1,field2,val2,...}
function c = LocalToCell(s)
f = fieldnames(s);
v = struct2cell(s);
opts = cell(2,length(f));
opts(1,:) = f;
opts(2,:) = v;
c = {opts{:}};

function c = LocalIsHG(obj,hgtype)
c = 0;
if (length(obj) == 1) & ishandle(obj)
  c = strcmp(get(obj,'type'),hgtype);
end

function c = LocalHas3DPlot(a)
zticks = LocalGetAsCell(a,'ZTickLabel');
c = 0;
for k=1:length(zticks)
  if ~isempty(zticks{k})
    c = 1;
    return;
  end
end

function r = LocalUnionRect(r1,r2)
if isempty(r1)
  r = r2;
elseif isempty(r2)
  r = r1;
elseif max(r2(3:4)) > 0
  left = min(r1(1),r2(1));
  bot = min(r1(2),r2(2));
  right = max(r1(1)+r1(3),r2(1)+r2(3));
  top = max(r1(2)+r1(4),r2(2)+r2(4));
  r = [left bot right-left top-bot];
else
  r = r1;
end

function r = LocalAxesTightBoundingBox(axesR, a)
if strcmp(get(a,'handlevisibility'),'on')
  r = get(a,'position');
else
  r = [];
end
atext = findall(a,'type','text','visible','on');
if ~isempty(atext)
  set(atext,'units','points');
  res=LocalGetAsCell(atext,'extent');
  for n=1:length(atext)
    r = LocalUnionRect(r,res{n});
  end
end
if strcmp(get(a,'visible'),'on')
  r = LocalUnionRect(r,[0 0 axesR(3:4)]);
  oldunits = get(a,'fontunits');
  set(a,'fontunits','points');
  label = text(0,0,'','parent',a,...
	       'units','points',...
	       'fontsize',get(a,'fontsize'),...
	       'fontname',get(a,'fontname'),...
	       'fontweight',get(a,'fontweight'),...
	       'fontangle',get(a,'fontangle'),...
	       'visible','off');
  fs = get(a,'fontsize');

  % handle y axis tick labels
  ry = [0 0 0 axesR(4)];
  ylabs = get(a,'yticklabels');
  yticks = get(a,'ytick');
  maxw = 0;
  if ~isempty(ylabs)
    for n=1:size(ylabs,1)
      if strcmp(get(a,'yscale'),'log')
	set(label,'string',['10^' ylabs(n,:)]);
      else
	set(label,'string',ylabs(n,:));
      end
      ext = get(label,'extent');
      maxw = max(maxw,ext(3));
    end
    if strcmp(get(a,'yaxislocation'),'left')
      ry(1) = -(maxw+5);
    else
      ry(1) = axesR(3);
    end
    ry(2) = -ext(4)/2;
    ry(3) = maxw+5;
    ry(4) = ry(4) + ext(4);
    if ~isempty(yticks)
      if ~strcmp(get(a,'yscale'),'log') & ...
            ((str2num(ylabs(1,:)) ~= yticks(1)) | ...
             (str2num(ylabs(end,:)) ~= yticks(end)))
        set(label,'string','10^1');
        ext = get(label,'extent');
        if strcmp(get(a,'xaxislocation'),'bottom')
          ry(4) = ry(4) + ext(4);
        end
      end
    end
    r = LocalUnionRect(r,ry);
  end

  % handle x axis tick labels
  rx = [0 0 axesR(3) 0];
  xlabs = get(a,'xticklabels');
  xticks = get(a,'xtick');
  if ~isempty(xlabs)
    if strcmp(get(a,'xscale'),'log')
      set(label,'string',['10^' xlabs(1,:)]);
    else
      set(label,'string',xlabs(1,:));
    end
    ext1 = get(label,'extent');
    rx(1) = -ext1(3)/2;
    if strcmp(get(a,'xscale'),'log')
      set(label,'string',['10^' xlabs(size(xlabs,1),:)]);
    else
      set(label,'string',xlabs(size(xlabs,1),:));
    end
    ext2 = get(label,'extent');
    rx(4) = max(ext1(4),ext2(4));
    rx(2) = -rx(4);
    rx(3) = axesR(3) + (ext2(3) + ext1(3))/2;
    if strcmp(get(a,'xaxislocation'),'top')
      rx(2) = rx(2) + axesR(4);
    end
    if ~isempty(xticks)
      if ~strcmp(get(a,'xscale'),'log') & ...
            ((str2num(xlabs(1,:)) ~= xticks(1)) | ...
             (str2num(xlabs(end,:)) ~= xticks(end)))
        set(label,'string','10^1');
        ext = get(label,'extent');
        if strcmp(get(a,'xaxislocation'),'bottom')
          rx(4) = rx(4) + ext(4)+5;
          rx(2) = rx(2) - ext(4)-5;
        else
          rx(4) = rx(4) + ext(4)+5 + axesR(4);
          rx(2) = rx(2) - ext(4)-5 - axesR(4);
        end
      end
    end
    r = LocalUnionRect(r,rx);
  end
  set(a,'fontunits',oldunits);
  delete(label);
end

function c = LocalManualAxesMode(old, allAxes, base)
xs = ['X' base];
ys = ['Y' base];
zs = ['Z' base];
xlog = LocalGetAsCell(allAxes,'xscale');
ylog = LocalGetAsCell(allAxes,'yscale');
zlog = LocalGetAsCell(allAxes,'zscale');
nonlogxflag = logical([]);
nonlogyflag = logical([]);
nonlogzflag = logical([]);
for k = 1:length(xlog)
  nonlogxflag(k) =  logical(~strcmp(xlog{k},'log'));
  nonlogyflag(k) =  logical(~strcmp(ylog{k},'log'));
  nonlogzflag(k) =  logical(~strcmp(zlog{k},'log'));
end
allNonLogXAxes = allAxes(nonlogxflag);
allNonLogYAxes = allAxes(nonlogyflag);
allNonLogZAxes = allAxes(nonlogzflag);
oldXMode = LocalGetAsCell(allNonLogXAxes,xs);
oldYMode = LocalGetAsCell(allNonLogYAxes,ys);
oldZMode = LocalGetAsCell(allNonLogZAxes,zs);
old = LocalPushOldData(old, allNonLogXAxes, {xs}, oldXMode);
old = LocalPushOldData(old, allNonLogYAxes, {ys}, oldYMode);
old = LocalPushOldData(old, allNonLogZAxes, {zs}, oldZMode);
set(allNonLogXAxes,xs,'manual');
set(allNonLogYAxes,ys,'manual');
set(allNonLogZAxes,zs,'manual');
c = old;

function val = LocalCheckAuto(val, auto)
if ischar(val) & strcmp(val,'auto')
  val = auto;
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
m = 100;
n = 1000;
k = 12;
A = spx.dict.simple.gaussian_dict(m, n);
gen = spx.data.synthetic.SparseSignalGenerator(n, k);
% create a sparse vector
x =  gen.biGaussian();
b = A*x;
A  = double(A);
result = spx.fast.omp(A, b, k, 1e-12);
cmpare = spx.commons.SparseSignalsComparison(x, result, k);
cmpare.summarize();

% we now run both algorithms 100 times
tic;
for i=1:1000
    result = spx.fast.omp(A, b, k, 1e-12);
end
elapsed_time1 = toc;
fprintf('Time taken %0.4f seconds\n', elapsed_time1);

tic;
for i=1:1000
    result = spx.pursuit.single.omp_chol(A, b, k, 1e-12);
end
elapsed_time2 = toc;
fprintf('Time taken %0.4f seconds\n', elapsed_time2);

improvement_factor = elapsed_time2 / elapsed_time1;
fprintf('improvement_factor: %0.2f \n', improvement_factor);

author: indigits

repo: sparse-plex


text:
classdef LMOCSO < ALGORITHM
% <multi/many> <real/integer> <large/none> <constrained/none>
% Large-scale multi-objective competitive swarm optimization algorithm

%------------------------------- Reference --------------------------------
% Y. Tian, X. Zheng, X. Zhang, and Y. Jin, Efficient large-scale multi-
% objective optimization based on a competitive swarm optimizer, IEEE
% Transactions on Cybernetics, 2020, 50(8): 3696-3708.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            [V,Problem.N] = UniformPoint(Problem.N,Problem.M);
            Population    = Problem.Initialization();
            Population    = EnvironmentalSelection(Population,V,(Problem.FE/Problem.maxFE)^2);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                Fitness = calFitness(Population.objs);
                if length(Population) >= 2
                    Rank = randperm(length(Population),floor(length(Population)/2)*2);
                else
                    Rank = [1,1];
                end
                Loser  = Rank(1:end/2);
                Winner = Rank(end/2+1:end);
                Change = Fitness(Loser) >= Fitness(Winner);
                Temp   = Winner(Change);
                Winner(Change) = Loser(Change);
                Loser(Change)  = Temp;
                Offspring      = Operator(Problem,Population(Loser),Population(Winner));
                Population     = EnvironmentalSelection([Population,Offspring],V,(Problem.FE/Problem.maxFE)^2);
            end
        end
    end
end

function Fitness = calFitness(PopObj)
% Calculate the fitness by shift-based density

    N      = size(PopObj,1);
    fmax   = max(PopObj,[],1);
    fmin   = min(PopObj,[],1);
    PopObj = (PopObj-repmat(fmin,N,1))./repmat(fmax-fmin,N,1);
    Dis    = inf(N);
    for i = 1 : N
        SPopObj = max(PopObj,repmat(PopObj(i,:),N,1));
        for j = [1:i-1,i+1:N]
            Dis(i,j) = norm(PopObj(i,:)-SPopObj(j,:));
        end
    end
    Fitness = min(Dis,[],2);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function ptch = showsurface( voxels )
%SHOWSURFACE: draw a surface based on some voxels
%
%   SHOWSURFACE(VOXELS) tries to render the supplied voxel structure as a
%   surface using MATLAB's ISOSURFACE command.
%
%   PTCH = SHOWSURFACE(VOXELS) also returns handles to the patches
%   created.

%   Copyright 2005-2009 The MathWorks, Inc.
%  $Revision: 1.0 $    $Date: 2006/06/30 00:00:00 $

% First grid the data
ux = unique(voxels.XData);
uy = unique(voxels.YData);
uz = unique(voxels.ZData);

% Expand the model by one step in each direction
ux = [ux(1)-voxels.Resolution; ux; ux(end)+voxels.Resolution];
uy = [uy(1)-voxels.Resolution; uy; uy(end)+voxels.Resolution];
uz = [uz(1)-voxels.Resolution; uz; uz(end)+voxels.Resolution];

% Convert to a grid
[X,Y,Z] = meshgrid( ux, uy, uz );

% Create an empty voxel grid, then fill only those elements in voxels
V = zeros( size( X ) );
N = numel( voxels.XData );
for ii=1:N
    ix = (ux == voxels.XData(ii));
    iy = (uy == voxels.YData(ii));
    iz = (uz == voxels.ZData(ii));
    V(iy,ix,iz) = voxels.Value(ii);
end

% Now draw it
ptch = patch( isosurface( X, Y, Z, V, 0.5 ) );
isonormals( X, Y, Z, V, ptch )
set( ptch, 'FaceColor', 'g', 'EdgeColor', 'none' );

set(gca,'DataAspectRatio',[1 1 1]);
xlabel('X');
ylabel('Y');
zlabel('Z');
view(-140,22)
lighting( 'gouraud' )
camlight( 'right' )
axis( 'tight' )

author: Sable

repo: mcbench-benchmarks


text:
function y = vl_nnpdist(x, x0, p, varargin)
%VL_NNPDIST CNN p-distance from target.
%   VL_NNPDIST(X, X0, P) computes the P distance raised of each feature
%   vector in X to the corresponding feature vector in X0:
%
%     Y(i,j,1) = (SUM_d (X(i,j,d) - X0(i,j,d))^P)^(1/P)
%
%   X0 should have the same size as X; the outoput Y has the same
%   height and width as X, but depth equal to 1. Optionally, X0 can
%   be a 1 x 1 x D x N array, in which case the same target feature
%   vector in X0 is compared to all feature vectors in X.
%
%   Setting the `noRoot` option to `true` does not take the 1/P power
%   in the formula, computing instead
%
%     Y(i,j,1) = SUM_d (X(i,j,d) - X0(i,j,d))^P
%
%   For example, `vl_nnpdist(x, x0, 2, 'noRoot', true)` computes the
%   squared L2 distance.
%
%   DZDX = VL_NNPDISTP(X, X0, P, DZDY) computes the derivative of the
%   block projected onto DZDY. DZDX and DZDY have the same dimensions
%   as X and Y, respectively.
%
%   VL_NNPDIST(___, 'OPT', VAL, ...) accepts the following options:
%
%   `NoRoot`:: `false`
%      If set to true, compute the P-distance to the P-th power.
%
%   `Epsilon`:: 1e-6
%      When computing derivatives, quantities that are divided in are
%      lower boudned by this value. For example, the L2 distance is
%      not smooth at the origin; this option prevents the
%      derivative from diverging.
%
%   `Aggregate`:: false
%      Instead of returning one scalar for each spatial location in
%      the inputs, sum all of them into a single scalar.
%
%   `InstanceWeights``:: `[]`
%      Optionally weight individual instances. This parameter can be
%      eigther a scalar or a weight mask, one for each pixel in the
%      input tensor.

% Copyright (C) 2015  Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% -------------------------------------------------------------------------
%                                                             Parse options
% -------------------------------------------------------------------------

opts.noRoot = false ;
opts.epsilon = 1e-6 ;
opts.aggregate = false ;
opts.instanceWeights = [] ;
backMode = numel(varargin) > 0 && ~ischar(varargin{1}) ;
if backMode
  dzdy = varargin{1} ;
  opts = vl_argparse(opts, varargin(2:end)) ;
else
  dzdy = [] ;
  opts = vl_argparse(opts, varargin) ;
end

% -------------------------------------------------------------------------
%                                                             Parse options
% -------------------------------------------------------------------------

d = bsxfun(@minus, x, x0) ;

if ~isempty(dzdy) && ~isempty(opts.instanceWeights)
  dzdy = bsxfun(@times, opts.instanceWeights, dzdy) ;
end

if ~opts.noRoot
  if isempty(dzdy)
    if p == 1
      y = sum(abs(d),3) ;
    elseif p == 2
      y = sqrt(sum(d.*d,3)) ;
    else
      y = sum(abs(d).^p,3).^(1/p) ;
    end
  else
    if p == 1
      y = bsxfun(@times, dzdy, sign(d)) ;
    elseif p == 2
      y = max(sum(d.*d,3), opts.epsilon).^(-0.5) ;
      y = bsxfun(@times, bsxfun(@times, dzdy, y),  d) ;
    elseif p < 1
      y = sum(abs(d).^p,3).^((1-p)/p) ;
      y = bsxfun(@times, bsxfun(@times, dzdy, y), max(abs(d), opts.epsilon).^(p-1) .* sign(d)) ;
    else
      y = max(sum(abs(d).^p,3), opts.epsilon).^((1-p)/p) ;
      y = bsxfun(@times, bsxfun(@times, dzdy, y), abs(d).^(p-1) .* sign(d)) ;
    end
  end
else
  if isempty(dzdy)
    if p == 1
      y = sum(abs(d),3) ;
    elseif p == 2
      y = sum(d.*d,3) ;
    else
      y = sum(abs(d).^p,3) ;
    end
  else
    if p == 1
      y = bsxfun(@times, dzdy, sign(d)) ;
    elseif p == 2
      y = bsxfun(@times, 2 * dzdy, d) ;
    elseif p < 1
      y = bsxfun(@times, p * dzdy, max(abs(d), opts.epsilon).^(p-1) .* sign(d)) ;
    else
      y = bsxfun(@times, p * dzdy, abs(d).^(p-1) .* sign(d)) ;
    end
  end
end

if isempty(dzdy)
  if ~isempty(opts.instanceWeights)
    y = bsxfun(@times, opts.instanceWeights, y) ;
  end
  if opts.aggregate
    y = sum(sum(y)) ;
  end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% SCRIPT TEST FOR THE KUKA LBR ROBOT KINEMATICS

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
robot = load_robot('KUKA', 'LBR_IIWA_R820_COP');
close all

% test both solutions: based on the transpose an on the Moore-Penrose
q = -pi/2*[1,1,1,1,1,1,1]';
drawrobot3d(robot, q)
targetposition = [0.5, -0.5, 0.5]';
vmag = 0.3;
delta_time = 50/1000;
qs = []
for i=1:20
    T = directkinematic(robot, q);
    J = manipulator_jacobian(robot, q);
    J
    error = targetposition-T(1:3, 4);
    vref = vmag*error/norm(error);
    wref = [0 0 0]';
    vwref = [vref; wref];
    iJ = pinv(J)
    qd = iJ*vwref;
    q = q + qd*delta_time;
    qs = [qs q];
    drawrobot3d(robot, q)
    pause(0.01)
    norm(error)
    if norm(error) < 0.01
        break
    end
end

norm(error)
figure, plot(qs', 'LineWidth',3.0)
legend('q1', 'q2','q3','q4','q5','q6','q7')

repo: ARTE


text:
function obj=stat_bin2d(obj,varargin)
% stat_bin2d() Makes 2D bins of X and Y data and displays count
%
% Parameters as 'name',value pairs:
% - 'nbins': Array in the form of [nxbins nybins] to set the
% number of bins in each dimension
% - 'edges': Cell in the form of {[x__edges] [y_edges]} to set
% custom bin edges for each dimension
% - 'geom': Set how results are displayed. 'image' uses a
% heatmap (default), 'contour' uses a contour plot. 'point'
% uses circles of varying size.

p=inputParser;
my_addParameter(p,'nbins',[30 30]);
my_addParameter(p,'edges',{});
my_addParameter(p,'geom','image'); %contour

parse(p,varargin{:});

obj.geom=vertcat(obj.geom,{@(dobj,dd)my_bin2d(dobj,dd,p.Results)});
obj.results.stat_bin2d={};
end

function hndl=my_bin2d(obj,draw_data,params)

x=comb(draw_data.x);
y=comb(draw_data.y);

if isempty(params.edges)
    [N,C] = hist3([shiftdim(x),shiftdim(y)],params.nbins);
else
    [N,C] = hist3([shiftdim(x),shiftdim(y)],'Edges',params.edges);

    obj.plot_lim.minx(obj.current_row,obj.current_column)=params.edges{1}(1);
    obj.plot_lim.maxx(obj.current_row,obj.current_column)=params.edges{1}(end);
    obj.plot_lim.miny(obj.current_row,obj.current_column)=params.edges{2}(1);
    obj.plot_lim.maxy(obj.current_row,obj.current_column)=params.edges{2}(end);

    %Put values on the upper edges as if they were in the last
    %bin
    N(:,end-1)=N(:,end-1)+N(:,end);
    N(end-1,:)=N(end-1,:)+N(end,:);

    %Remove upper edge
    N(:,end)=[];
    N(end,:)=[];

end

obj.results.stat_bin2d{obj.result_ind,1}.edges=C;
obj.results.stat_bin2d{obj.result_ind,1}.counts=N;

switch params.geom
    case 'contour'

        [~,hndl]=contour(C{1},C{2},N',5,'Color',draw_data.color);

    case 'image'

        if ~obj.continuous_color_options.active
            obj.continuous_color_options.active = true;
        end

        Nr=reshape(N',1,numel(N));
        sel=Nr>0;
        %sel=true(size(Nr));

        if isempty(params.edges)
            %Get polygon half widths
            wx=(C{1}(2)-C{1}(1))/2;
            wy=(C{2}(2)-C{2}(1))/2;

            %Generate polygon edges
            [X,Y] = meshgrid(C{1},C{2});

            X=reshape(X,1,numel(X));
            Y=reshape(Y,1,numel(Y));

            patchesx=[X(sel)-wx ; X(sel)-wx ; X(sel)+wx ; X(sel)+wx ];
            patchesy=[Y(sel)-wy ; Y(sel)+wy ; Y(sel)+wy ; Y(sel)-wy ];

        else
            [Xs, Ys]=meshgrid(params.edges{1}(1:end-1),params.edges{2}(1:end-1));
            [Xe, Ye]=meshgrid(params.edges{1}(2:end),params.edges{2}(2:end));

            Xs=reshape(Xs,1,numel(Xs));
            Ys=reshape(Ys,1,numel(Ys));
            Xe=reshape(Xe,1,numel(Xe));
            Ye=reshape(Ye,1,numel(Ye));

            patchesx=[Xs(sel) ; Xs(sel) ; Xe(sel) ; Xe(sel)];
            patchesy=[Ys(sel) ; Ye(sel) ; Ye(sel) ; Ys(sel)];

            %If we have varied-size patches (we use rounding to
            %get away with numerical issues of unique
            if length(unique(round(diff(params.edges{1})*1e10)))>1 || length(unique(round(diff(params.edges{2})*1e10)))>1
                %Correct values by the area of each patch ?
                Nr=Nr./((Xe-Xs).*(Ye-Ys));
                obj.aes_names.color='Count/area';
            else
                obj.aes_names.color='Count';
            end

        end

        %patchesz=[Nr(sel) ; Nr(sel) ; Nr(sel) ; Nr(sel) ];

        %p=patch(patchesx,patchesy,patchesz,Nr(sel));
        hndl=patch(patchesx,patchesy,Nr(sel));
        set(hndl,'edgeColor','none')

        %Store color values
        obj.plot_lim.maxc(obj.current_row,obj.current_column)=max(Nr);
        obj.plot_lim.minc(obj.current_row,obj.current_column)=min(Nr);

    case 'point'
        [X,Y] = meshgrid(C{1},C{2});
        X=reshape(X,1,numel(X));
        Y=reshape(Y,1,numel(Y));
        Nr=reshape(N',1,numel(N));
        sel=Nr>0;
        %hndl=point_patch(X(sel),Y(sel),Nr(sel)*(C{1}(2)-C{1}(1))/30,draw_data.color,20,ratio);
        hndl=scatter(X(sel),Y(sel),Nr(sel)*10,draw_data.marker,'MarkerEdgeColor',draw_data.color,'MarkerFaceColor','none');
end

obj.results.stat_bin2d{obj.result_ind,1}.handle=hndl;

end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text:
% clear all;
% close all;
% clc;

display('Nth (2) order Reconstruction');

N = 2;                  % Nth order nonuniform sampling
TQ = 1;Fs = 1/TQ;                 % Nyquist Period

T = [1.5*TQ 3*TQ];      % Decimation Periods
K = 0.5*lcm(2*T(1), 2*T(2))/TQ; capT = K*TQ; M = capT./T;
ML = (100:100:2500); % number of slices
w_c = 0.85;
NS = 100;  % Number of Sinusoids

LF = lcm(M(1),M(2))*2*K+1;      % length of LF should be Multiple of LCM{M(p)}*2*K
n = -(LF-1)/2:1:(LF-1)/2;

std = 1e-1;
taus = [0 1+std*randn]*TQ;
tausI = sort([taus(1) taus(2) T(1)+taus(1)]);
display(tausI);

Frq = rand(1,NS)*w_c/2;
Amp = rand(1,NS)/(sqrt(NS)*2);
Phi = rand(1,NS)*2*pi;

timeP = zeros(size(ML));
timeE = zeros(size(ML));
timeI = zeros(size(ML));
timeV = zeros(size(ML));
timePr = zeros(size(ML));
timeJ = zeros(size(ML));

MC_runs = 1;
for rrr = 1:length(ML)
display(rrr);
input = zeros(1,ML(rrr)*K);
for k = 1:NS
  input = input + Amp(k)*sin(2*pi*Frq(k)*(0:ML(rrr)*K-1)*TQ+Phi(k));
end;
for tt = 1:MC_runs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A realization of Digital Filter Banks for Reconstruction of Uniformly
% sampled signals from nonuniform samples
% Authors: Itami, Watanabe, Nishihara
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tauI = zeros(K,ML(rrr));
for p = 1:K
    tauI(p,:) = tausI(p)+(0:ML(rrr)-1)*capT;
end;

x11 = zeros(K,ML(rrr));
for k = 1:NS
    x11 = x11 + Amp(k)*sin(2*pi*Frq(k)*tauI+Phi(k));
end;

tic
a = zeros(1,K);
for p = 1:K
    a(p) = 1;
    for q = 1:K
        if q ~= p
                a(p) = a(p)/sin(pi*(tausI(p)-tausI(q))/capT);
        end;
    end;
end;
% c = sin(pi*(tausI(0+1))/capT);
% s = cos(pi*(tausI(0+1))/capT);
% b(1,1) = 0.5*(c+1i*s);
% c = sin(pi*(tausI(1+1))/capT);
% s = cos(pi*(tausI(1+1))/capT);
% b(1,2) = 0.5*(c+1i*s);
% b(3,:) = conj(b(1,:));
% b(2,:) = 0;

b = zeros(2*K-1,K);
c = -0.5*cos(pi*(tausI(1+1)+tausI(2+1))/capT);
s = -0.5*sin(pi*(tausI(1+1)+tausI(2+1))/capT);
b(1,1) = 0.5*(c+1i*s);
c = -0.5*cos(pi*(tausI(0+1)+tausI(2+1))/capT);
s = -0.5*sin(pi*(tausI(0+1)+tausI(2+1))/capT);
b(1,2) = 0.5*(c+1i*s);
c = -0.5*cos(pi*(tausI(0+1)+tausI(1+1))/capT);
s = -0.5*sin(pi*(tausI(0+1)+tausI(1+1))/capT);
b(1,3) = 0.5*(c+1i*s);
b(5,:) = conj(b(1,:));
b(2,:) = 0;
b(4,:) = conj(b(2,:));
c = 0.5*cos(pi*(tausI(2+1)-tausI(1+1))/capT);
b(3,1) = c;
c = 0.5*cos(pi*(tausI(2+1)-tausI(0+1))/capT);
b(3,2) = c;
c = 0.5*cos(pi*(tausI(1+1)-tausI(0+1))/capT);
b(3,3) = c;

y1=upsample(x11.',K).';

k = -(K-1):1:(K-1);
m = (0:1:(2*K-1))';
F = exp(1i*(pi/K).*kron(m,k));

y = zeros(K,size(y1,2));
for r = 1:K
    tempI = a(r)*y1(r,:);
    tempI2 = F*b(:,r);
    y2 = tempI2*tempI;
%     h = sinc((n/K)-tausI(r)/capT).*knab(LF,4,-tausI(r)/TQ).';
    h = sinc((n/K)-tausI(r)/capT).*kaiser_mine1(LF,3,-tausI(r)/TQ);
    for i=1:2*K
            h1 = upsample(downsample(h,2*K,i-1),2*K);
            temp = filter([zeros(1,i-1),1],1,h1);
%             y = y + filter(temp,1,y2(i,:));
            y2(i,:) = filter(temp,1,y2(i,:));
%             y2(i,:) = filter(h1,1,y2(i,:));
%             y2(i,:) = filter([zeros(1,i-1),1],1,y2(i,:));%,zeros(1,2*K-i)
    end;
    y(r,:) = sum(y2,1);
end;
y = (sum(y,1));
timeI(rrr) = timeI(rrr)+toc;
delayI = (length(n)-1)/2;
x=input(1:end-delayI);
y=y(1+delayI:end);
y = y(160:end);
x = x(160:end);
serI = 20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of Nonuniformly Sampled Band-Limited Signals
% Using a Differentiator-Multiplier Cascade
% Authors: Stefan Tertinek and Christian Vogel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r = tausI-TQ*(0:K-1);
rr = r(mod((0:ML(rrr)*K-1),K)+1);
x1 = reshape(x11,1,K*size(x11,2));
% x1 = x11;% LF = (0:11); Fs = 1;
% Differentiator Design
% figure();
% NFFT = 2^nextpow2(length(Hd)); % Next power of 2 from length of Hd
% HD = fftshift(fft(Hd,NFFT))/length(Hd);
% f = Fs*linspace(-1,1,NFFT);
% % Plot double-sided amplitude spectrum.
% plot(f,2*abs(HD(1:NFFT)))
% title('Double-Sided Amplitude Spectrum of Hd(n)')
% xlabel('Frequency (Hz)')
% ylabel('|HD(f)|')
tic
Hd = firpm(LF-1,[0 w_c],[0 w_c*pi],'differentiator');
delayV = (LF-1)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y1 = filter(Hd,1,x1);
x1 = filter([zeros(1,delayV),1],1,x1);
r2 = filter([zeros(1,delayV),1],1,rr);
e = y1.*r2;
y1 = x1-e;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y2 = filter(Hd,1,y1);
temp = filter([zeros(1,delayV),1],1,y2);
x1 = filter([zeros(1,2*delayV),1],1,x1);
r2 = filter([zeros(1,2*delayV),1],1,r2);
e1 = temp.*r2;
y2 = filter(Hd,1,y2);
e2 = 0.5*y2.*r2.^2;
y2 = x1-e1-e2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y3 = filter(Hd,1,y2);
temp = filter([zeros(1,2*delayV),1],1,y3);
x1 = filter([zeros(1,3*delayV),1],1,x1);
r2 = filter([zeros(1,3*delayV),1],1,r2);
e1 = temp.*r2;
y3 = filter(Hd,1,y3);
temp = filter([zeros(1,delayV),1],1,y3);
e2 = 0.5*temp.*r2.^2;
y3 = filter(Hd,1,y3);
e3 = (y3.*r2.^3)/6;
y3 = x1-e1-e2-e3;
timeV(rrr) = timeV(rrr)+toc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y = real(y3(1+6*delayV:end));
x = input(1:end-6*delayV);
serV = 20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of Band-Limited Periodic Nonuniformly Sampled Signals
% Through Multirate Filter Banks
% Ryan S Prendergast, Bernard C Levy, Paul J Hurst
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tauPr = zeros(K,ML(rrr));
for p = 1:K
    tauPr(p,:) = -tausI(p)+(0:ML(rrr)-1)*capT;
end;
x1 = zeros(K,ML(rrr));
for k = 1:NS
    x1 = x1 + Amp(k)*sin(2*pi*Frq(k)*tauPr+Phi(k));
end;
tic
H = zeros(K,LF);
for i=1:K
%     H(i,:) = sinc(n-tausI(i)).*conv(sinc(n-tausI(i)),kaiser(LF,34).','same');
%     H(i,:) = sinc(n-tausI(i)).*knab(LF,34,-tausI(i)).';
    H(i,:) = sinc(n-tausI(i)).*kaiser_mine1(LF,3,-tausI(i));
end;
H = H(:,1:end-1);
EP = reshape(H.',K,length(H(1,:))/K,K);
capE = [];
for k=1:K
    temp = [];
    for i=1:K
        temp=[temp,toeplitz([EP(i,1,k),zeros(1,(length(H(1,:))/K)-1)],[EP(i,:,k),zeros(1,(length(H(1,:))/K)-1)])];
    end
    capE = [capE;temp];
end
d = ceil(size(capE,2)/(2*K));
P = kron(eye(K),[zeros(1,d),1,zeros(1,(size(capE,2)/K)-d-1)]);
R = P/capE;
% size(capE)
% size(zeros(LF-1,2*LF-2-K))
% size(R)
% size(zeros(K,LF-1))
% size(P)
% size(zeros(K,2*LF-2-K))
R = upsample(R.',K).';
rt = size(R,2)/K;
FR = zeros(K,rt);
for j=1:K
    for i=1:K
        temp = filter([zeros(1,K-i),1],1,R(i,(j-1)*rt+1:j*rt));
        FR(j,:) = FR(j,:)+temp;
    end
end
yb = upsample(x1.',K).';
for i = 1:K
    yb(i,:) = filter(FR(i,:),1,yb(i,:));
end;
y = sum(yb,1);
timePr(rrr) = timePr(rrr)+toc;
delayPr = (size(FR,2))/2+K-1;
x=input(1:end-delayPr);
y=y(1+delayPr:end);
serPr = 20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of Periodically Nonuniformly Sampled Bandlimited Signals
% Using Time-Varying FIR Filters
% Authors: H. Johansson and Per Lowenborg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r = r.';
x1 = reshape(x11,1,size(x11,2)*K);
tic
w_o = w_c*pi*TQ;
hJ = zeros(K,LF);
C = zeros(1,LF);
Nt = (LF-1)/2;
for i = 1:K
    C = -2*sin(w_o*(n-r(1+(mod(i-1-n,K)))'))./(pi*(n-r(1+(mod(i-1-n,K)))'));
    C(isnan(C)==1)=-2*w_o/pi;
    C = C.';
    S = zeros(LF,LF);
    for k = 1:LF
        S(k,:) = sin(w_o*(-Nt+k-1-r(1+(mod(i-1-(-Nt+k-1),K)))-(n-r(1+(mod(i-1-n,K)))')))./(pi*(-Nt+k-1-r(1+(mod(i-1-(-Nt+k-1),K)))-(n-r(1+(mod(i-1-n,K)))')));
    end;
    S(isnan(S)==1)=w_o/pi;
    hJ(i,:) = -0.5*S\C;
end;

y1 = zeros(K,length(x1));
for j=1:K
    y1(j,:) = filter(hJ(j,:),1,x1);
    y1(j,:) = upsample(downsample(y1(j,:),K,j-1),K)/K;
    y1(j,:) = filter([zeros(1,j-1),1],1,y1(j,:));
end;
y = K*0.25*sum(y1,1);
timeJ(rrr) = timeJ(rrr)+toc;
delayJ = (size(hJ,2)-1)/2;
y = real(y(1+delayJ:end));
x = input(1:end-delayJ);
serJ = 20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filterbank Reconstruction of Bandlimited Signals from Nonuniform and
% Generalized Samples
% Authors: Y C Eldar and A V Oppenheim
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

y = zeros(N,ML(rrr)*K);
for p = 1:N
    tau = taus(p)+(0:ML(rrr)*M(p)-1)*T(p);
    x1 = zeros(1,ML(rrr)*M(p));
    for k = 1:NS
        x1 = x1 + Amp(k)*sin(2*pi*Frq(k)*tau+Phi(k));
    end;
    tic
    y1 = upsample(x1,K);

    LFE = M(p)*lcm(M(1),M(2))*2*K+1;      % length of LF should be Multiple of LCM{M(p)}*2*K
    nE = -(LFE-1)/2:1:(LFE-1)/2;
    h = sinc((nE/K)-(taus(p)/T(p))).*kaiser_mine1(LFE,3,-K*(taus(p)/T(p)));

%%%%%%%%%%%%%%%%%%%%%%%%%%% Implementation 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%
%     bb = zeros(2*(K-M(p))+1,M(p));
%     aaa = ones(1,M(p));
%     bbvl = zeros(M(p),LFE);
%     for l = 1:M(p)
%         for q = 1:N
%             if q ~= p
%                 aaa(l) = aaa(l)*sin(pi*(taus(p)-taus(q)+(l-1)*T(p))/T(q));
%                 c = sin(pi*M(q)*(-taus(q)+(l-1)*T(p))/capT);
%                 s = cos(pi*M(q)*(taus(q)-(l-1)*T(p))/capT);
%                 bb(2*M(q)+1,l) = 0.5*(c-1i*s);
%                 bb(1,l) = conj(bb(2*M(q)+1,l));
%             end;
%         end;
%         bbv = zeros(2*(K-M(p))+1,LFE);
%         for v = -(K-M(p)):K-M(p)
%             bbv(K-M(p)+1+v,:) = bb(K-M(p)+1+v,l)*exp(1i*(pi/(K*M(p)))*v*nE);
%         end;
%         bbvl(l,:) = sum(bbv,1)/aaa(l);
%     end;
%     bbb = zeros(M(p),LFE);
%     bbn = zeros(M(p),LFE);
%     for m = 1:M(p)
%         for l = 1:M(p)
%             bbb(l,:) = bbvl(l,:)*exp(1i*(2*pi/M(p))*(m-1)*(l-1));
%         end
%         bbb = sum(bbb,1);
%         bbn(m,:) = bbb.*exp(1i*(2*pi/M(p))*(m-1)*nE);
%     end;
%     bbn = sum(bbn,1);
%     bbn = bbn.*h/M(p);

%%%%%%%%%%%%%%%%%%%%%%%%%%% Implementation 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%
%     bb = zeros(2*(K-M(p))+1,M(p));
%     aaa = ones(1,M(p));
%     for l = 1:M(p)
%         for q = 1:N
%             if q ~= p
%                 aaa(l) = aaa(l)/sin(pi*(taus(p)-taus(q)+(l-1)*T(p))/T(q));
%                 c = sin(pi*M(q)*(-taus(q)+(l-1)*T(p))/capT);
%                 s = cos(pi*M(q)*(taus(q)-(l-1)*T(p))/capT);
%                 bb(2*M(q)+1,l) = 0.5*(c-1i*s);
%                 bb(1,l) = conj(bb(2*M(q)+1,l));
%             end;
%         end;
%     end;
%     AE = diag(aaa);
%     BE = bb.';
%
%     dim = K-M(p); w = -dim:1:dim;
%     FE = exp(1i*(pi/(K*M(p))).*kron(nE,w'));
%
%     E1E = exp(1i*(2*pi/M(p)).*kron((0:M(p)-1),(0:M(p)-1)'))/M(p);
%
%     E2E = exp(1i*(2*pi/M(p)).*kron(nE,(0:M(p)-1)'));
%
%     temp = E1E*AE*BE*FE;
%     bbn = h.*sum(E2E.*temp,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%% Implementation 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%
    aaa = ones(1,M(p));
    bb = ones(M(p),LFE);
    for l = 1:M(p)
        for q = 1:N
            if q ~= p
                aaa(l) = aaa(l)*sin(pi*(taus(p)-taus(q)+(l-1)*T(p))/T(q));
                bb(l,:) = bb(l,:).*sin(pi*((nE*TQ/M(p))-taus(q)+(l-1)*T(p))/T(q));
            end;
        end;
        bb(l,:) = bb(l,:)/aaa(l);
    end;

    bbb = zeros(M(p),LFE);
    bbn = zeros(M(p),LFE);
    for m = 1:M(p)
        for l = 1:M(p)
            bbb(l,:) = bb(l,:)*exp(1i*(2*pi/M(p))*(m-1)*(l-1));
        end
        bbb = sum(bbb,1);
        bbn(m,:) = bbb.*exp(1i*(2*pi/M(p))*(m-1)*nE);
    end;
    bbn = sum(bbn,1);
    bbn = bbn.*h/M(p);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    y1 = conv(y1,bbn);
    delayE = (length(h)-1)/2;
    y(p,:) = y1(1+delayE:M(p):end-delayE);
    timeE(rrr) = timeE(rrr)+toc;
end;
y = sum(real(y),1);
x = input;
serE = 20*log10(norm(x,2)/norm(y-x,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reconstruction of N-th order nonuniformly sampled bandlimited signals
% using digital filter banks
% Authors: S. K. Sindhi, K. M. M. Prabhu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xp = zeros(N,ML(rrr)*K);
for p = 1:N
    tau = taus(p)+(0:ML(rrr)*M(p)-1)*T(p);
    x1 = zeros(1,ML(rrr)*M(p));
    for k = 1:NS
        x1 = x1 + Amp(k)*sin(2*pi*Frq(k)*tau+Phi(k));
    end;

    tic
    m = (0:1:M(p)-1)'; lemda = 0:1:M(p)-1;
    W = exp(1i*(2*pi/M(p)).*kron(m,lemda)); % m*lemda

    bb = zeros(2*(K-M(p))+1,M(p));
    aaa = ones(1,M(p));
    for l = 1:M(p)
        for q = 1:N
            if q ~= p
                    aaa(l) = aaa(l)/sin(pi*M(q)*(taus(p)-taus(q)+(l-1)*T(p))/capT);
                    c = sin(pi*M(q)*(-taus(q)+(l-1)*T(p))/capT);
                    s = cos(pi*M(q)*(taus(q)-(l-1)*T(p))/capT);
                    bb(2*M(q)+1,l) = 0.5*(c-1i*s);
                    bb(1,l) = conj(bb(2*M(q)+1,l));
            end;
        end;
    end;
    A = diag(aaa); % display(A);
    B = bb; % display(B);

    y1 = upsample(x1,K);
    y1 = reshape(y1,M(p),length(y1)/M(p));

    if M(p)>1
        y1(2:end,:) = flipud(y1(2:end,:));
        for i = 1:M(p)-1
            y1(i+1,:) = filter([0,1],1,y1(i+1,:));
        end;
    end;

    dim = K-M(p); w = -dim:1:dim;

    xlemda = zeros(M(p),size(y1,2));
    for lemda = 0:M(p)-1

        rP = (lemda/M(p))+(0:1:(2*K-1))';
        Fshift = exp(1i*(pi/K).*kron(rP,w));   % r*w

%         h = sinc((n*TQ/T(p))+(lemda/K)-(taus(p)/T(p))).*knab(LF,4,(lemda/M(p))-(taus(p)/TQ)).';
        h = sinc((n*TQ/T(p))+(lemda/K)-(taus(p)/T(p))).*kaiser_mine1(LF,3,(lemda/M(p))-(taus(p)/TQ));
        mtemp = A*W;
        mtemp = B*mtemp;
        mtemp = Fshift*mtemp;
%         mtemp = mtemp*W(:,lemda+1);
%         y2 = mtemp*y1(lemda+1,:);
%         for i=1:2*K
%             h1 = upsample(downsample(h,2*K,i-1),2*K);
% %             y2(i,:) = filter(h1,1,y2(i,:));
% %             y2(i,:) = filter([zeros(1,i-1),1],1,y2(i,:));%,zeros(1,2*K-i)
%             temp = filter([zeros(1,i-1),1],1,h1);
%             y2(i,:) = filter(temp,1,y2(i,:));
%         end;
%         xlemda(lemda+1,:) = sum(y2,1);

        Htemp = mtemp*W(:,lemda+1);
        h1 = zeros(2*K, length(h)+2*K-1);
        for i = 2:2*K
            h1(i,:) = [filter([zeros(1,i-1),1],1,upsample(downsample(h,2*K,i-1),2*K)) zeros(1,2*K)]*Htemp(i);
        end;
        h1(1,:) = upsample(downsample(h,2*K),2*K)*Htemp(1);
        h1 = sum(h1,1);
        xlemda(lemda+1,:) = filter(h1,1,y1(lemda+1,:));
    end;
    xp(p,:) = sum(xlemda,1)/M(p);
    clear bb;
    timeP(rrr) = timeP(rrr)+toc;
end;
y = real(sum(xp,1));
delayP = (length(h)-1)/2;
y = y(1+delayP:end);
x = input(1:end-delayP);
serP = 20*log10(norm(x,2)/norm(y-x,2));
end;
end;
timeP = timeP/MC_runs;
timeE = timeE/MC_runs;
timeI = timeI/MC_runs;
timeV = timeV/MC_runs;
timePr = timePr/MC_runs;
timeJ = timeJ/MC_runs;

figure();hold on;
plot(ML,timeJ,'kp-','LineWidth',2);
plot(ML,timePr,'ko-','LineWidth',2);
plot(ML,timeV,'ks-','LineWidth',2);
plot(ML,timeP,'kd-','LineWidth',2);
plot(ML,timeI,'k>-','LineWidth',2);
plot(ML,timeE,'k+-','LineWidth',2);
legend('Johansson','Prendergast','Tertinek','Proposed','Itami','Eldar');
xlabel('Input signal length','fontsize',14,'fontweight','b');
ylabel('Time in seconds','fontsize',14,'fontweight','b');
grid on;box on;
set(gca,'fontsize',14,'fontweight','b')

display(serI)
display(serV)
display(serPr)
display(serJ)
display(serP)
display(serE)
% figure();
% subplot(2,1,1);
% plot(([x' y']));
% title('input / output signals');
% xlabel('sample');
% ylabel('signal value');
% grid on;
% subplot(2,1,2);
% plot((x'-y'));
% xlabel('time (sample)');
% ylabel('error value');
% grid on;

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function [ X ] = uniformCircle(c,R,N)
%  X  = uniformCircle(c,R,N)
% inputs : c = [c1,c2] center of the circle, R radius, N number of points
% output : X of size Nx2, points uniformly distributed on the cricle
theta = 2*pi*rand(N,1);
X = [c(1)+R*cos(theta) c(2)+R*sin(theta)];

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [c] = power(a,b,varargin)
% C = A.^B
%   [C]=POWER(A,B) Compute A.^B for a TT-tensor A and natural B.
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

if isa(a,'tt_tensor') && numel(b) == 1
    c=a;
    b=b-1;
    while(b>0)
      c=c.*a;
      b=b-1;
    end

else
    error('Incorrect usage for TT-power');
end

author: oseledets

repo: TT-Toolbox


text:
classdef TestTickMeter
    %TestTickMeter

    methods (Static)
        function test_1
            tm = cv.TickMeter();
            tm.reset();
            tm.start();
            pause(0.1);
            tm.stop();
            validateattributes(tm.TimeTicks, {'int64'}, {'scalar'});
            validateattributes(tm.TimeMicro, {'double'}, {'scalar'});
            validateattributes(tm.TimeMilli, {'double'}, {'scalar'});
            validateattributes(tm.TimeSec, {'double'}, {'scalar'});
            validateattributes(tm.Counter, {'int64'}, {'scalar'});
            %assert(abs(tm.TimeSec - 0.1) < 0.02);
        end

        function test_2
            t1 = cv.TickMeter.getTickCount();
            pause(0.1);
            t2 = cv.TickMeter.getTickCount();
            f = cv.TickMeter.getTickFrequency();
            tsec = double(t2 - t1) / f;
            validateattributes(t1, {'int64'}, {'scalar'});
            validateattributes(t2, {'int64'}, {'scalar'});
            validateattributes(f, {'double'}, {'scalar'});
            %assert(abs(tsec - 0.1) < 0.02);
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


function models=train_svms_box(imnames, feat_dir, region_meta_info, cachedir, run_name, varargin)
%random seed
rng(3);

%input parser to allow additional inputs
p = inputParser;
addOptional(p, 'to_train', [1:20], @isnumeric);
addOptional(p, 'svm_C', 10^-3, @isscalar);
addOptional(p,'bias_mult',10,@isscalar);
addOptional(p, 'pos_loss_weight', 2, @isscalar);
addOptional(p, 'do_latent', 0, @isscalar);
addOptional(p, 'pos_ov_thresh', 0.5, @isscalar);
addOptional(p, 'neg_ov_thresh', 0.3, @isscalar);
addOptional(p, 'hard_neg_thresh', -1.0001, @isscalar);
addOptional(p, 'evict_thresh', -1.2, @isscalar);
addOptional(p, 'retrain_limit', 2000, @isscalar);
parse(p, varargin{:});
opts=p.Results;

to_train=opts.to_train;

%if cachedir does not exist, create it
if(~exist(cachedir, 'file')) mkdir(cachedir); end
cachedir=fullfile(cachedir, run_name);
if(~exist(cachedir, 'file')) mkdir(cachedir); end

%collect feature stats
fprintf('Collecting feature stats...\n');
stats=collect_feature_stats(imnames, feat_dir, cachedir, region_meta_info);

%initialize caches
for i=1:numel(to_train)
    caches(i) = init_cache;
    models(i).categid=to_train(i);
    models(i).w=[];
    models(i).b=[];
    models(i).meannrm=stats.meannrm;
end

%load GT positives
pos_file=fullfile(cachedir, 'pos.mat');
if(exist(pos_file, 'file'))
    tmp=load(pos_file);
    X_pos=tmp.X_pos;
else
    X_pos=get_positives(imnames, feat_dir, to_train, region_meta_info, false, models, opts.pos_ov_thresh);
    save(pos_file, 'X_pos');
end
for i=1:numel(to_train)
    caches(i).X_pos=X_pos{i};
    fprintf('Category %d has %d positives\n', to_train(i), size(X_pos{i},1));
end

%train
first_time=true;
force_update=false;
max_outer_iter=opts.do_latent+1;
for iter=1:max_outer_iter
    if(iter~=max_outer_iter)
        %subsample negatives
        negidx=randperm(numel(imnames));
        negidx=negidx(1:min(1000, numel(negidx)));
    else
        negidx=1:numel(imnames);
    end
    if(iter>1)
        %latent update
        X_pos=get_positives(imnames, feat_dir, to_train, region_meta_info, true, models, opts.pos_ov_thresh);
        for i=1:numel(to_train)
            caches(i).X_pos=X_pos{i};
            fprintf('Latent update: Category %d has %d positives\n', to_train(i), size(X_pos{i},1));

        end
        force_update=true;
    end
    %for every negative
    for i=1:numel(negidx)
        fprintf('Hard negatives iter: %d, image number: %d/%d\n',iter, i, numel(negidx));
        [X_neg, curr_keys]=get_hard_negatives(negidx(i), imnames, feat_dir, to_train,...
                             region_meta_info, models, {caches.keys}, opts.neg_ov_thresh, opts.hard_neg_thresh, first_time);
       very_last_iter=(iter==max_outer_iter) && (i==numel(negidx));
       for j=1:numel(to_train)
            caches(j).X_neg=cat(1, caches(j).X_neg,X_neg{j});
            caches(j).keys=cat(1, caches(j).keys, curr_keys{j});
            caches(j).num_added=caches(j).num_added+size(X_neg{j},1);
            if(caches(j).num_added>opts.retrain_limit || first_time ||very_last_iter || force_update)
                %time to retrain
                fprintf('Retraining category %d with %d positives and %d negatives\n', to_train(j), size(caches(j).X_pos,1), size(caches(j).X_neg,1));
                models(j)=update_model(models(j), caches(j), opts.svm_C, opts.bias_mult, opts.pos_loss_weight);
                [reg, posloss, negloss, neg_scores]=compute_obj_val(models(j),...
                                                   caches(j), opts.svm_C, opts.pos_loss_weight);
                caches(j).pos_loss(end+1)=posloss;
                caches(j).neg_loss(end+1)=negloss;
                caches(j).reg_loss(end+1)=reg;
                caches(j).tot_loss(end+1)=posloss+negloss+reg;
                for t=1:numel(caches(j).tot_loss)
                    fprintf('%d: %f = %f + %f + %f\n',t,caches(j).tot_loss(t),...
                    caches(j).reg_loss(t),caches(j).pos_loss(t), caches(j).neg_loss(t));
                end

                %evict
                keep=neg_scores>opts.evict_thresh;
                caches(j).X_neg=caches(j).X_neg(keep,:);
                caches(j).keys=caches(j).keys(keep,:);
                fprintf('Kept: %d negatives \n', sum(keep));
                caches(j).num_added=0;
            end

        end
        first_time=false;
        force_update=false;
    end
end

%--------------------------------
%   Feature stats
%--------------------------------
function stats = collect_feature_stats(imnames, feat_dir, cachedir, region_meta_info)
stat_file=fullfile(cachedir, 'stats.mat');
if(exist(stat_file, 'file'))
    tmp=load(stat_file);
    stats=tmp.stats; clear tmp;
else

    num_images=100;
    regs_per_img=100;
    ri=randperm(numel(imnames));
    ri=ri(1:num_images);
    totalnrm=0;
    totalcnt=0;
    for i=1:num_images
        fprintf('%d/%d\n', i, num_images);

        %load random image
        tmp=load(fullfile(feat_dir, [imnames{ri(i)} '.mat']));

        %pick random regions
        num=min(regs_per_img, region_meta_info.num_regs(ri(i)));
        ri2=randperm(region_meta_info.num_regs(ri(i)));
        ri2=ri2(1:num);

        %get norm
        nrm=sqrt(sum(tmp.feats(ri2,:).^2,2));
        totalnrm=totalnrm+sum(nrm);
        totalcnt=totalcnt+num;
    end
    stats.meannrm=totalnrm./totalcnt;
    save(stat_file, 'stats');
end
fprintf('\n\nMean norm=%f\n\n',stats.meannrm);

%------------------------------
%   Init cache
%------------------------------
function cache=init_cache
cache.X_pos = single([]);
cache.X_neg = single([]);
cache.keys = [];
cache.num_added = 0;
cache.retrain_limit = 2000;
cache.evict_thresh = -1.2;
cache.hard_thresh = -1.0001;
cache.pos_loss = [];
cache.neg_loss = [];
cache.reg_loss = [];
cache.tot_loss = [];

%------------------------------
%   Get positives from image
%------------------------------

function X_pos=get_positives(imnames, featdir, to_train, region_meta_info, islatent, models, pos_ov_thresh)
for k=1:numel(to_train)
    X_pos{k}=single([]);
end

for i=1:numel(imnames)
    fprintf('Positives: %d/%d\n',i, numel(imnames));
    tmp=load(fullfile(featdir, [imnames{i} '.mat']));
    feats=xform_feat(tmp.feats, models(1).meannrm);

    %for every category
    for j=1:numel(to_train)
        sel=[];
        if(islatent)
            %latent update: consider all features greater than thresh
            %find gt for this category and get corresponding overlaps
            gtidx=find(region_meta_info.gt{i}==to_train(j));
            nongt=find(region_meta_info.gt{i}==0);

            ovs=region_meta_info.box_overlaps{i}(gtidx,:);
            overlapping=find(max(ovs,[],1)>=pos_ov_thresh);
            ovs=ovs(:,overlapping);

            scores=feats*models(j).w+models(j).b;

            %for every gt, find the highest scoring gt that overlaps by >70
            for k=1:numel(gtidx)
                ovs_this=ovs(k,:);
                idx=find(ovs_this>pos_ov_thresh);
                if(isempty(idx)) continue; end
                [m, argmax]=max(scores(nongt(overlapping(idx))));
                sel=[sel nongt(overlapping(idx(argmax)))];
            end
        else
            sel=find(region_meta_info.gt{i}==to_train(j));
        end
        X_pos{j}=cat(1, X_pos{j},feats(sel,:));
    end
end

%------------------------------
%   Hard negatives
%------------------------------

function [X_neg, curr_keys]=get_hard_negatives(i, imnames, featdir, to_train, region_meta_info, models, keys, neg_ov_thresh, score_thresh, first_time)
tmp=load(fullfile(featdir, [imnames{i} '.mat']));
feats=xform_feat(tmp.feats, models(1).meannrm);

if(~first_time)
    %concatenate all models
    w=cat(2,models(:).w);
    b=cat(2,models(:).b);

    %compute all scores
    scores=bsxfun(@plus,feats*w, b);
else
    scores=(score_thresh+1)*ones(size(feats,1), numel(models));
end

%for every category
for j=1:numel(to_train)
    %get the groundtruths that don't belong to this category
    wronggtidx=find(region_meta_info.gt{i}~=to_train(j) & region_meta_info.gt{i}~=0);
    %get the non-groundtruths that overlap by less than threshold

    nongt=find(region_meta_info.gt{i}==0);
    gtidx=find(region_meta_info.gt{i}==to_train(j));
    if(~isempty(gtidx))
        idx=find(max(region_meta_info.box_overlaps{i}(gtidx,:),[],1)<neg_ov_thresh);
        nongtidx=nongt(idx);
    else
        nongtidx=nongt;
    end
    allidx=[wronggtidx(:); nongtidx(:)];
    allidx=allidx(scores(allidx,j)>=score_thresh);

    %get the keys for each of these
    tmp_keys=[i*ones(numel(allidx),1) allidx];

    %check that duplicate keys are not added
    keep=~ismember(tmp_keys, keys{j}, 'rows');

    %add these feats
    X_neg{j}=feats(allidx(keep),:);
    curr_keys{j}=tmp_keys(keep,:);
end

%------------------------------
%   Update model
%------------------------------
function model=update_model(model, cache, svm_c, bias_mult, pos_weight)
opts = sprintf('-w1 %.10f -c %.10f -s 3 -B %.10f', ...
                      pos_weight, svm_c, bias_mult);

X=cat(2, cache.X_pos', cache.X_neg');
y=[ones(size(cache.X_pos,1),1); zeros(size(cache.X_neg,1),1)];
fprintf('calling liblinear with:%s\n',opts);
llm = liblinear_train(y, sparse(double(X)), opts, 'col');
model.w = single(llm.w(1:end-1)');
model.b = single(llm.w(end)*bias_mult);

%---------------------------
% Compute objective value
%---------------------------
function [reg, posloss, negloss, neg_scores]=compute_obj_val(model,cache, svm_c, pos_weight)
pos_scores=cache.X_pos*model.w+model.b;
neg_scores=cache.X_neg*model.w+model.b;
posloss=sum(max(1-pos_scores,0))*svm_c*pos_weight;
negloss=sum(max(1+neg_scores,0))*svm_c;
reg=0.5*model.w'*model.w;

%---------------------------
% xform_feat
%---------------------------
function feats=xform_feat(feats, nrm)
feats=feats*20./nrm;


text:
% FPMat Interface
% Example of Surface, Contour plot data

% Read Surface, Contour plot data.
% The data has been generated in FlexPDE
% by running Example01.PDE

fpreadsc('TempC.flx');
fpreadsc('TempS.flx');
fpreadsc('VS.flx');


text:
function output = callqpas(interfacedata)

% Retrieve needed data
options = interfacedata.options;
model = yalmip2quadprog(interfacedata);

if size(model.Aeq,1)==0
    model.Aeq = [];
    model.beq = [];
end

if options.savedebug
    save debugfile model
end

if options.showprogress;showprogress(['Calling ' interfacedata.solver.tag],options.showprogress);end
solvertime = tic;
[x,flag,lm] = qpas(model.Q, model.c, model.A, model.b, model.Aeq, model.beq, model.lb, model.ub, options.verbose);
solvertime = toc(solvertime);

% Internal format for duals
if ~isempty(lm)
    D_struc = [lm.equality;lm.inequality];
else
    D_struc = [];
end

if isempty(flag)
    problem = 9;
else
    % Check, currently not exhaustive...
    switch flag
        case 0
            problem = 0;
        case 2
            problem = 1;
        otherwise
            problem = 1;
    end
end

% Save all data sent to solver?
if options.savesolverinput
    solverinput = model;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.x = x;
    solveroutput.flag = flag;
    solveroutput.lm = lm;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x(:),D_struc,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
clear all; close all; clc;

%%
L=10; x3=-L:0.1:L; n=length(x3)-1; % define domain
x2=x3(1:n); k=(2*pi/(2*L))*[0:n/2-1 -n/2:-1]; % k-vector
ye=exp(-(x2.^2)); ye2=exp((x2.^2)/2); % define Gaussians
for j=0:9      % loop through 10 modes
  yd=real(ifft(((i*k).^j).*fft(ye))); % 2nd derivative
  mode=((-1)^(j))*(((2^j)*factorial(j)*sqrt(pi))^(-0.5))*ye2.*yd;
  y(:,j+1)=(mode).';  % store modes as columns
end

x=x2(n/2+1-40:n/2+1+40);  % keep only -4<x<4
yharm=y(n/2+1-40:n/2+1+40,:);

%%
n=length(x);
f=(exp(-(x-0.5).^2)+3*exp(-2*(x+1.5).^2))';

for j=1:10  % full reconstruction
  a(j,1)=trapz(x,f.*yharm(:,j));
end
f2=yharm*a;
Efull(1)=log(norm(f2-f)+1);  % reconstruction error

for j=1:10  % matrix M reconstruction
   for jj=1:j
       Area=trapz(x,yharm(:,j).*yharm(:,jj));
       M(j,jj)=Area;
       M(jj,j)=Area;
   end
end
Cfull=cond(M)   % get condition number

%% Test Random trials with P% of measurements
for thresh=1:5;
  for jloop=1:1000  % 1000 random trials
    n2=randsample(n,8*thresh);  % random sampling
    P=zeros(n,1); P(n2)=1;
    for j=1:10
        for jj=1:j   % compute M matrix
            Area=trapz(x,P.*(yharm(:,j).*yharm(:,jj)));
            M2(j,jj)=Area; M2(jj,j)=Area;
        end
    end
    for j=1:10  % reconstruction using gappy
        ftild(j,1)=trapz(x,P.*(f.*yharm(:,j)));
    end
    atild=M2\ftild;   % compute error
    f2=yharm*atild;   % compute reconstruction
    Err(jloop)=norm(f2-f); % L2 error
    con(jloop)=cond(M2);   % condition number
  end
  % mean and variance
  E(thresh)=mean(log(Err+1)); V(thresh)=(var(log(Err+1)));
  Ec(thresh)=mean(log(con)); Vc(thresh)=(var(log(con)));
end
E=[E Efull]; V=[V 0];
Ec=[Ec log(Cfull)]; Vc=[Vc 0];
subplot(2,1,1), bar(E), hold on, errorbar(E,V,'r.')
subplot(2,1,2), bar(Ec), hold on, errorbar(Ec,Vc,'r.')

%%  For 20% measurements, sort great from bad

Pmaster=[];
for jloop=1:200  % 200 random trials
n2=randsample(n,20);
P=zeros(n,1); P(n2)=1;

Pmaster=[Pmaster P];

for j=1:10
   for jj=1:j
       Area=trapz(x,P.*(yharm(:,j).*yharm(:,jj)));
       M2(j,jj)=Area; M2(jj,j)=Area;
   end
end
for j=1:10  % reconstruction using gappy
  ftild(j,1)=trapz(x,P.*(f.*yharm(:,j)));
end
atild=M2\ftild;   % compute error
f2=yharm*atild;
Er(jloop)=norm(f2-f);
co(jloop)=cond(M2);

end

%%
figure(5)
subplot(3,1,1), bar(log(co))
subplot(3,1,2), hist(log(Er+1),50)
subplot(3,1,3), hist(log(co),50)

figure(7)
pcolor(-Pmaster.'), colormap(hot)

[Cloc,jloc]=sort(co);
figure(8)
%bar(log(Cloc));

jbest=jloc(1:11); jworst=jloc(end-10:end);
Sbest(:,1:11)=Pmaster(:,jbest);
Sworst(:,1:11)=Pmaster(:,jworst);
csor=[100*Cloc(1:10) zeros(1,10) Cloc(end-9:end)];

subplot(3,1,1), pcolor(-Sbest.'); colormap(hot)
subplot(3,1,2), pcolor(-Sworst.'); colormap(hot)
subplot(3,1,3), bar(csor)

%%
figure(7), axis off
figure(8), subplot(3,1,1), axis off
figure(8), subplot(3,1,2), axis off
figure(8), subplot(3,1,3), set(gca,'Xticklabel',{},'Yticklabel',{},'Xlim',[0 31])
figure(5), subplot(3,1,1), set(gca,'Xticklabel',{},'Yticklabel',{},'Xlim',[0 200])
figure(5), subplot(3,1,2), set(gca,'Xticklabel',{},'Yticklabel',{})
figure(5), subplot(3,1,3), set(gca,'Xticklabel',{},'Yticklabel',{})
figure(3), subplot(2,1,1), set(gca,'Xticklabel',{},'Yticklabel',{})
figure(3), subplot(2,1,2), set(gca,'Xticklabel',{},'Yticklabel',{})

%%

clear q
clear q2
figure(9)
m=30;
q=rand(m,m); n=m^2;
p=randsample(n,200);
q2=zeros(m,m);
q2(p)=q(p);
p2=randsample(n,40);
q3=zeros(m,m);
q3(p2)=q(p2);
pcolor(-q), colormap(hot), caxis([-1 0]), axis off
figure(10)
pcolor(-q2), colormap(hot), caxis([-1 0]), axis off
figure(11)
pcolor(-q3), colormap(hot), caxis([-1 0]), axis off

author: dynamicslab


text:
%% Rotation Demonstration: Discontinous Data
%
% A demonstration of how the timestep affects the error of the splitting
% method for a rotation. Similar results can be expected in general,
% although nonlinear sharpening may give smaller errors.
%
% The equation we solve is
%
% $$u_t - y u_x + x u_y = 0, \qquad u(x,y,0)=u_0(x,y)$$
%
% where the initial function u0 is discontinuous. To solve the equation, we
% will use Godunov splitting and front tracking to approximate the 1-D
% solution operators.

%% Initial setup
T     = 4*pi;
xmin  = -1.5; xmax=1.5;
N     = 128;
h     = (xmax-xmin)/N;
x     = xmin+h*(0:N);
y     = 0.5*(x(1:end-1)+x(2:end));
[X,Y] = meshgrid(y,y);
u0    = adiscontnuous_function(X,Y);
surfl(X,Y,u0); shading interp, view(-15,60), colormap(gray), axis tight

%% Error versus CFL number
% First we study the pointwise and the L1 error for CFL numbers 2.^(0:5)
colormap(jet)
for i=1:6
	nu=2^(i-1);
	err= abs( rotrack(u0,y,y,nu,T) - u0);
   subplot(2,3,i), pcolor(y,y,err), axis equal image, shading interp;
   set(colorbar('West'),'XColor',[1 1 1],'YColor',[1 1 1]);
   caxis([0 2]), title(['CFL = ', num2str(nu)]);
   set(gca,'XTick',[],'YTick',[])
   xlabel(['L1 error: ', num2str(h*h*sum(abs(err(:))))]);
	drawnow;
end;
%%
% The error is determined by two error mechanisms that work in opposite
% directions. The splitting error increases with increasing splitting
% steps, whereas the smoothing error caused by the projection operator
% increases with decreasing splitting steps. The total splitting error is
% therefore a convex function of the time step, and in the figure above,
% the minimum error is observed for CFL number 16.

%% Testing convergence
% We choose the best CFL number vu=16, and check for convergence rate.
% To this end, we use only one lap (it takes too long otherwise).
T    = 2*pi;
Lerr = zeros(1,7);
wbar = waitbar(0,'Convergence study');
for i=1:7,
	N    = 2^(i+3);
	h    = (xmax-xmin)/N; x=xmin+h*(1:N);
	[X,Y]=meshgrid(x,x);
	u0   = adiscontnuous_function(X,Y);
	err  = abs(rotrack(u0,x,x,nu,T) - u0);
	Lerr(i)=sum(err(:)) / sum(abs(u0(:)));
	waitbar(i/7,wbar);
end;
close(wbar);
clf, semilogy(1:7,Lerr,'o--');
set(gca,'XTick',1:7,'XTickLabel',2.^((1:7)+3))
xlabel('Grid size'); ylabel('Log(error)'); axis tight
title('Errors for \nu = 16');

repo: Aero-matlab


text:
function [s, S_l] = hmgLin2seg(l)

%   HMGLIN2SEG HMG line to segment conversion
%   HMGLIN2SEG  returns a 3d segment with the two suppor points of the HMG
%   line L.
%
%   [s, S_l] = IDPLIN2SEG(L) returns the Jacobian wrt L.

%   Copyright 2009 Teresa Vidal.

if nargout == 1

    [e1, e2] = hmgLin2hmgPnts(l);
    p1       = hmg2euc(e1);
    p2       = hmg2euc(e2);

    s        = [p1;p2];

else

    [e1, e2, E1_l, E2_l] = hmgLin2hmgPnts(l);
    [p1, P1_e1]          = hmg2euc(e1);
    [p2, P2_e2]          = hmg2euc(e2);

    s   = [p1;p2];
    S_l = [P1_e1*E1_l ; P2_e2*E2_l];

end

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb