text:
function [cfg] = ft_spike_plot_raster(cfg, spike, timelock)

% FT_SPIKE_PLOT_RASTER makes a raster plot of spike-trains and allows for a
% spike-density or a PSTH plot on top.
%
% Use as
%   ft_spike_plot_raster(cfg, spike)
% or
%   ft_spike_plot_raster(cfg, spike, timelock)
%
% The input SPIKE data structure should be organized as the spike or the
% raw datatype The optional input TIMELOCK should be organized as the
% timelock datatype, e.g. the output from FT_SPIKE_PSTH or FT_SPIKEDENSITY,
% having the average firing rate / spike count per time-point / time-bin.
% However, timelock could also be the output from FT_TIMELOCKANALYSIS.
%
% Configuration options
%   cfg.spikechannel     =  see FT_CHANNELSELECTION for details
%   cfg.latency          =  [begin end] in seconds, 'maxperiod' (default), 'minperiod',
%                           'prestim' (all t<=0), or 'poststim' (all t>=0).
%   cfg.linewidth        =  number indicating the width of the lines (default = 1);
%   cfg.cmapneurons      =  'auto' (default), or nUnits-by-3 matrix.
%                           Controls coloring of spikes and psth/density
%                           data if multiple cells are present.
%   cfg.spikelength      =  number >0 and <=1 indicating the length of the spike. If
%                           cfg.spikelength = 1, then no space will be left between
%                           subsequent rows representing trials (row-unit is 1).
%   cfg.trialborders     =  'yes' or 'no'. If 'yes', borders of trials are
%                           plotted
%   cfg.plotselection   =  'yes' or 'no' (default). If yes plot Y axis only for selection in cfg.trials
%   cfg.topplotsize      =  number ranging from 0 to 1, indicating the proportion of the
%                           rasterplot that the top plot will take (e.g., with 0.7 the top
%                           plot will be 70% of the rasterplot in size). Default = 0.5.
%   cfg.topplotfunc      =  'bar' (default) or 'line'.
%   cfg.errorbars        = 'no', 'std', 'sem' (default), 'conf95%','var'
%
%   cfg.interactive      = 'yes' (default) or 'no'. If 'yes', zooming and panning operate via callbacks.
%   cfg.trials           =  numeric or logical selection of trials (default = 'all').

% Copyright (C) 2010-2013, Martin Vinck
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% 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$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble provenance spike timelock

% check if input spike structure is indeed a spike structure
spike = ft_checkdata(spike,'datatype', 'spike', 'feedback', 'yes'); % converts raw as well

% get the default options
cfg.spikechannel = ft_getopt(cfg,'spikechannel', 'all');
cfg.trials     = ft_getopt(cfg,'trials', 'all');
cfg.latency      = ft_getopt(cfg,'latency','maxperiod');
cfg.linewidth    = ft_getopt(cfg,'linewidth', 1);
cfg.cmapneurons  = ft_getopt(cfg,'cmapneurons', 'auto');
cfg.spikelength  = ft_getopt(cfg,'spikelength', 0.9);
cfg.topplotsize  = ft_getopt(cfg,'topplotsize', 0.5);
cfg.topplotfunc  = ft_getopt(cfg,'topplotfunc', 'bar');
cfg.errorbars    = ft_getopt(cfg,'errorbars', 'sem');
cfg.trialborders = ft_getopt(cfg,'trialborders','yes');
cfg.plotselection = ft_getopt(cfg,'plotselection','no');
cfg.interactive  = ft_getopt(cfg,'interactive','yes');

% ensure that the options are valid
cfg = ft_checkopt(cfg,'spikechannel',{'cell', 'char', 'double'});
cfg = ft_checkopt(cfg,'latency', {'char', 'ascendingdoublebivector'});
cfg = ft_checkopt(cfg,'trials', {'char', 'doublevector', 'logical'});
cfg = ft_checkopt(cfg,'linewidth', 'doublescalar');
cfg = ft_checkopt(cfg,'cmapneurons', {'char', 'doublematrix', 'doublevector'});
cfg = ft_checkopt(cfg,'spikelength', 'doublescalar');
cfg = ft_checkopt(cfg,'topplotsize', 'doublescalar');
cfg = ft_checkopt(cfg,'topplotfunc', 'char', {'bar', 'line'});
cfg = ft_checkopt(cfg,'errorbars', 'char', {'sem', 'std', 'conf95%', 'no', 'var'});
cfg = ft_checkopt(cfg,'trialborders', 'char', {'yes', 'no'});
cfg = ft_checkopt(cfg,'plotselection', 'char', {'yes', 'no'});
cfg = ft_checkopt(cfg,'interactive', 'char', {'yes', 'no'});

cfg = ft_checkconfig(cfg, 'allowed', {'spikechannel', 'latency', 'trials', 'linewidth', 'cmapneurons', 'spikelength', 'topplotsize', 'topplotfunc', 'errorbars', 'trialborders', 'plotselection', 'interactive'});

% check if a third input is present, and check if it's a timelock structure
if nargin==3
  doTopData = true;
  timelock  = ft_checkdata(timelock, 'datatype', 'timelock', 'hastrials', 'no', 'feedback', 'yes');
else
  doTopData = false;
end

% get the spikechannels
cfg.spikechannel = ft_channelselection(cfg.spikechannel, spike.label);
spikesel         = match_str(spike.label, cfg.spikechannel);
nUnits           = length(spikesel); % number of spike channels
if nUnits==0, error('No spikechannel selected by means of cfg.spikechannel'); end

% get the number of trials and set the cfg.trials field
nTrialsOrig = size(spike.trialtime,1);
nTrialsShown = nTrialsOrig;
if  strcmp(cfg.trials,'all')
  cfg.trials = 1:nTrialsOrig;
elseif islogical(cfg.trials) || all(cfg.trials==0 | cfg.trials==1)
  cfg.trials = find(cfg.trials);
end
cfg.trials = sort(cfg.trials(:));

if max(cfg.trials)>nTrialsOrig
  error('maximum trial number in cfg.trials should not exceed length of spike.trial')
end
if isempty(cfg.trials), errors('No trials were selected in cfg.trials'); end

% determine the duration of each trial
begTrialLatency = spike.trialtime(cfg.trials,1);
endTrialLatency = spike.trialtime(cfg.trials,2);

% select the latencies
if strcmp(cfg.latency,'minperiod')
  cfg.latency = [max(begTrialLatency) min(endTrialLatency)];
elseif strcmp(cfg.latency,'maxperiod') || strcmp(cfg.latency,'all')
  cfg.latency = [min(begTrialLatency) max(endTrialLatency)];
elseif strcmp(cfg.latency,'prestim')
  cfg.latency = [min(begTrialLatency) 0];
elseif strcmp(cfg.latency,'poststim')
  cfg.latency = [0 max(endTrialLatency)];
end

% check whether the time window fits with the data
if (cfg.latency(1) < min(begTrialLatency)), cfg.latency(1) = min(begTrialLatency);
  warning('Correcting begin latency of averaging window');
end
if (cfg.latency(2) > max(endTrialLatency)), cfg.latency(2) = max(endTrialLatency);
  warning('Correcting begin latency of averaging window');
end

% delete trials that are not in our window
overlaps      = endTrialLatency>=cfg.latency(1) & begTrialLatency<=cfg.latency(2);
trialSel      = overlaps(:);
cfg.trials    = cfg.trials(trialSel); %update the trial selection
if isempty(cfg.trials),warning('No trials were selected');end

% create the data that should be plotted
[unitX,unitY] = deal(cell(1,nUnits));

% find the spiketimes per neuron and make one vector of them with y value per element
for iUnit = 1:nUnits
  unitIndx       = spikesel(iUnit);
  latencySel     = spike.time{unitIndx}>=cfg.latency(1) & spike.time{unitIndx} <= cfg.latency(2);
  isInTrials     = ismember(spike.trial{unitIndx},cfg.trials);
  unitX{iUnit}   = spike.time{unitIndx}(isInTrials(:) & latencySel(:));
  unitY{iUnit}   = spike.trial{unitIndx}(isInTrials(:) & latencySel(:));
  if strcmp(cfg.plotselection,'yes')
   tempY{iUnit} = zeros(size(unitY{iUnit}));
   u = unique(unitY{iUnit});
   for i = 1:length(u)
     idx = find(unitY{iUnit} == u(i));
     tempY{iUnit}(idx) = i;
   end
   nTrialsShown = length(u);
   unitY{iUnit} = tempY{iUnit};
  end
end

% some error checks on spike length
if (cfg.spikelength<=0 || cfg.spikelength>1)
  error('cfg.spikelength should be a single number >0 and <=1. 1 row (1 trial) = 1');
end

% some error checks on the size of the top figure
if (cfg.topplotsize<=0 || cfg.topplotsize>1)
  error('cfg.topplotsize should be a single number >0 and <=1. 0.7 = 70%');
end

% plot the trial borders if desired
if strcmp(cfg.trialborders,'yes')
  begTrialLatency = begTrialLatency(overlaps(:));
  endTrialLatency = endTrialLatency(overlaps(:));

  % make yellow fills after trial end
  X = endTrialLatency;
  Y = cfg.trials;
  Y = [Y(:);Y(end);Y(1);(Y(1))];
  X = [X(:);cfg.latency(2);cfg.latency(2);endTrialLatency(1)];
  f = fill(X,Y,'y');
  set(f,'EdgeColor', [1 1 1]);

  % make yellow fills before trial beginning
  X = begTrialLatency;
  Y = cfg.trials;
  Y = [Y(:);Y(end);Y(1);(Y(1))];
  X = [X(:);cfg.latency(1);cfg.latency(1);begTrialLatency(1)];
  hold on, fill(X,Y,'y'), hold on
  set(f,'EdgeColor', [1 1 1]);
end

% start plotting the rasterplot
for iUnit = 1:nUnits

  % create the x and y value to be plotted
  x = [unitX{iUnit}(:)'; unitX{iUnit}(:)']; % create x duplicates, spike times
  y = unitY{iUnit}(:); % these are the trial values
  y = [y' - cfg.spikelength/2; y' + cfg.spikelength/2];

  % process the color information for the different units
  if strcmp(cfg.cmapneurons, 'auto'), cfg.cmapneurons = colormap_cgbprb(nUnits); end
  if all(size(cfg.cmapneurons) ./ [nUnits 3])
    color = cfg.cmapneurons(iUnit,:);
  else
    error('cfg.cmapneurons should be nUnits-by-3 matrix or "auto"');
  end

  % create axes for the rasterplot, all go to the same position, so do this for unit 1
  if iUnit==1
    ax(1) = newplot;
    posOrig = get(ax(1), 'Position'); % original position for a standard plot
    pos     = posOrig;
    if doTopData % if topdata, we leave space on top
      posRaster   = [pos(1:3) pos(4)*(1-cfg.topplotsize)]; % and decrease the size of width and height
    else
      posRaster   = pos;
    end
    set(ax(1), 'ActivePositionProperty', 'position', 'Position', posRaster)
  end

  % make the raster plot and hold on for the next plots
  rasterHdl = line(x, y,'LineWidth', cfg.linewidth,'Color', color);
  cfg.hdl.raster = rasterHdl;
  set(ax(1),'NextPlot', 'add')
  set(ax(1),'Box', 'off')
end

% create the labels for the first axes
xlabel('time (sec)')
ylabel('Trial Number')
axis ij

% plot the top data
if doTopData

  % match on the basis of labels, and specify this in the documentary
  unitIndx = find(ismember(timelock.label,spike.label(spikesel)));
  if sum(unitIndx)<nUnits, error('timelock.label should contain all label of selected units'); end

  % select timepoints  and get the data to be plotted
  binSel = timelock.time>=cfg.latency(1) & timelock.time<=cfg.latency(2);
  dataY  = timelock.avg(unitIndx,binSel);
  dataX  = timelock.time(binSel);

  % create the position for the topplot and axes with this position
  posTopPlot = [0 pos(4)*(1-cfg.topplotsize)  0 0] + pos.*[1 1 1 cfg.topplotsize];
  ax(2)  = axes('Units', get(ax(1), 'Units'), 'Position', posTopPlot,...
    'Parent', get(ax(1), 'Parent'));

  % make the bar or the line plot
  if strcmp(cfg.topplotfunc,'bar')

    % plot thje density / psth as a bar plot
    avgHdl  = feval(cfg.topplotfunc,dataX,dataY', 'stack');

    % color the different bars
    for iUnit = 1:nUnits
      set(avgHdl(iUnit),'FaceColor',cfg.cmapneurons(iUnit,:),'EdgeColor', cfg.cmapneurons(iUnit,:),'LineWidth', 3);
    end

    % set the bar width to 1
    set(avgHdl,'BarWidth', 1);

    if ~strcmp(cfg.errorbars, 'no') && nUnits>1
      warning('error bars can only be plotted with bar plot if one unit is plotted');
    end

    % add standard error bars to it
    if ~strcmp(cfg.errorbars, 'no') && nUnits==1
      x = [timelock.time(binSel);timelock.time(binSel)]; % create x doublets
      if strcmp(cfg.errorbars, 'sem')
        err = sqrt(timelock.var(unitIndx,binSel)./timelock.dof(unitIndx,binSel));
      elseif strcmp(cfg.errorbars, 'std')
        err = sqrt(timelock.var(unitIndx,binSel));
      elseif strcmp(cfg.errorbars, 'var')
        err = timelock.var(unitIndx,binSel);
      elseif strcmp(cfg.errorbars, 'conf95%')
        % use a try statement just in case the statistics toolbox doesn't work.
        try
          tCrit = tinv(0.975,timelock.dof(unitIndx,binSel));
          err = tCrit.*sqrt(timelock.var(unitIndx,binSel)./timelock.dof(unitIndx,binSel)); % assuming normal distribution, SHOULD BE REPLACED BY STUDENTS-T!
        catch % this is in case statistics toolbox is lacking
          err = 0;
          warning('error with computing 95% conf limits, probably issue with statistics toolbox');
        end
      end

      y = dataY + err;

      % plot the errorbars as line plot on top of the bar plot
      hold on
      y = [zeros(1,length(y));y]; % let lines run from 0 to error values
      hd = plot(x,y,'k'); % plot the error bars
    end

  else

    % plot the density / psth
    avgHdl  = feval(cfg.topplotfunc,dataX,dataY');

    % set the color for every unit
    for iUnit = 1:nUnits
      set(avgHdl(iUnit),'Color',cfg.cmapneurons(iUnit,:));
    end

    % ensure that the data limits are increasing
    mnmax = [nanmin(dataY(:))-eps nanmax(dataY(:))+eps];

    % plot the errorbars
    if ~strcmp(cfg.errorbars,'no')

      % check if the right error information is there
      if ~isfield(timelock, 'var') || ~isfield(timelock, 'dof')
        error('timelock should contain field .var and .dof for errorbars');
      end

      % select the degrees of freedom
      df = timelock.dof(unitIndx,binSel);

      % plot the different error bars
      if strcmp(cfg.errorbars, 'sem')
        err = sqrt(timelock.var(unitIndx,binSel)./df);
      elseif strcmp(cfg.errorbars, 'std')
        err = sqrt(timelock.var(unitIndx,binSel));
      elseif strcmp(cfg.errorbars, 'var')
        err = timelock.var(unitIndx,binSel);
      elseif strcmp(cfg.errorbars, 'conf95%')
        tCrit = tinv(0.975,df);
        err   = tCrit.*sqrt(timelock.var(unitIndx,binSel)./df);
      end

      % ensure division by zero does not count, should not happen however
      err(~isfinite(err)) = NaN;

      % make a polygon of the error bars that allows easy filling in adobe
      for iUnit = 1:nUnits
        upb  = dataY(iUnit,:) + err(iUnit,:);
        lowb = dataY(iUnit,:) - err(iUnit,:);
        sl   = ~isnan(upb);
        [X,Y] = polygon(dataX(sl),upb(sl)+0.0001,lowb(sl)-0.0001);
        hold on
        hd = plot(X,Y,'--');
        set(hd,'Color', cfg.cmapneurons(iUnit,:));
        hold on
      end
      mnmax = [nanmin(dataY(:) - err(:)) nanmax(dataY(:) + err(:))];
    end

    % set the y limits explicitly
    ylim = mnmax;
    d    = ylim(2)-ylim(1);
    yl = ylim(1)-0.05*d;
    if ylim(1)>0
      if yl<0, yl = 0; end
    end
    ylim(1) = yl;

    yl = ylim(2)+0.05*d;
    if ylim(2)<0
      if yl>0, yl = 0; end
    end
    ylim(2) = yl;
  if ylim(2) == ylim(1) %if the plot is empty
    ylim(2)=ylim(1)+1; %be nice to set
  end
    set(gca,'YLim', ylim)
  end

  % store the handle
  cfg.hdl.psth = avgHdl;

  % modify the axes
  set(ax(2),'YAxisLocation', 'right') % swap y axis location
  set(ax(2),'XTickLabel', {}) % remove ticks and labels for x axis
  set(ax(2),'XTick', [])
  set(ax(2), 'Box', 'off') % turn the box off

  % change the axis settings
  try
    ylabel(timelock.cfg.outputunit)
  catch
    warning('unit of the y-axis is unknown');
    ylabel('Signal intensity')
  end
end

% set the limits for the axis
set(ax,'XLim', [cfg.latency])
if nTrialsShown==0; nTrialsShown = 1; end %
set(ax(1), 'YLim', [0.5 nTrialsShown+0.5]); % number of trials
set(ax,'TickDir','out') % put the tickmarks outside

% now link the axes, constrain zooming and keep ticks intact
limX       = [cfg.latency];
limY       = get(ax,'YLim');
if ~iscell(limY), limY = {limY}; end

% constrain the zooming and zoom psth together with the jpsth, remove
% ticklabels jpsth
if strcmp(cfg.interactive,'yes')
  set(zoom,'ActionPostCallback',{@mypostcallback,ax,limX,limY});
  set(pan,'ActionPostCallback',{@mypostcallback,ax,limX,limY});
end

% pass positions and axis handles so downstream
% functions can respect the positions set here.
cfg.pos.posRaster = posRaster;
cfg.hdl.axRaster = ax(1);
if doTopData
 cfg.pos.posTopPlot = posTopPlot;
 cfg.hdl.axTopPlot = ax(2);
end

% do the general cleanup and bookkeeping at the end of the function

ft_postamble previous spike
ft_postamble provenance

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [] = mypostcallback(fig,evd,ax,lim,ylim)
currentAxes = evd.Axes;
indx = find(currentAxes==ax);

% keep the y axis within boundaries
origLim = ylim{indx};
ylim = get(ax(indx), 'YLim');
if origLim(1)>ylim(1), ylim(1) = origLim(1); end
if origLim(2)<ylim(2), ylim(2) = origLim(2); end
set(ax(indx), 'YLim', ylim)

% keep the x axis within boundaries and reset for both
xlim = get(ax(indx), 'XLim');
if lim(1)>xlim(1), xlim(1) = lim(1); end
if lim(2)<xlim(2), xlim(2) = lim(2); end
set(ax,'XLim', xlim)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [colorspace] = colormap_cgbprb(N)
% COLORMAP_SEPARATION returns a color map with well separated colors that does not include
% white, yellow or orange since these are not well visible in white plots.
%
% Inputs: N is the number of colors desired
% Outputs: COLORSPACE is an N-by-3 colorspace.
% Chosen colors are green, cyan, blue, purple, red and black
% Script was adapted from varycolors.m from FEX
% Copyright (c) 2009, Martin Vinck

if nargin<1
  error('specify the number of colors that you want')
end

% create a well visible color space
%        green     cyan  blue    purple    red     black
colors = [[0 1 0]; [0 1 1] ; [0 0 1] ; [1 0 1]; [1 0 0]; [0 0 0]];
order  = [1 5 3 2 4 6]; % our preference order when plotting different lines
rm{2} = [2:5];
rm{3} = [2 4 5];
rm{4} = [2 4];
rm{5} = [4];

if N==1
  colorspace = colors(end,:);
elseif N>1&&N<6
  colors(rm{N},:) = [];
  colorspace = colors;
else
  n      = floor(N/5)*ones(1,5);
  rest   = mod(N,5);
  order  = [1 5 3 2 4];
  % if we have some rest, add them starting oppositly
  n(order(1:rest)) = n(order(1:rest)) + 1;

  colorspace = [];
  for iColor = 1:5
    corStart  = colors(iColor,:);
    corEnd    = colors(iColor+1,:);
    dim       = corStart~=corEnd;
    subColors = corStart(ones(n(iColor)+1,1),:);
    if iColor>1
      subColors(:,dim) = linspace(corStart(dim),corEnd(dim),n(iColor)+1);
      subColors(1,:) = [];
    else
      subColors(1:end-1,dim) =   linspace(corStart(dim),corEnd(dim),n(iColor));
      subColors(end,:) = [];
    end
    colorspace = [colorspace;subColors];
  end
end

function [X,Y] = polygon(x,sm1,sm2,multiplier)

x = x(:)';
if nargin < 4
    multiplier = 1;
end
X = [x x(end:-1:1) x(1)];
up = sm1(:)';
down = sm2(:)';
Y = [down up(end:-1:1) down(1)];

repo: fieldtrip


text:
function p = cs_nd (A)
%CS_ND generalized nested dissection ordering.
%   p = cs_nd(A) computes the nested dissection ordering of a matrix.  Small
%   submatrices (order 500 or less) are ordered via cs_amd.  A must be sparse
%   and symmetric (use p = cs_nd(A|A') if it is not symmetric).
%
%   Example:
%       A = delsq (numgrid ('L', 300)) ;    % matrix used in 'bench'
%       p = cs_nd (A) ;
%       cspy (A (p,p)) ;
%
%   See also CS_AMD, CS_SEP, CS_ESEP, CS_NSEP, AMD.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

n = size (A,1) ;
if (n == 1)
    p = 1 ;
elseif (n < 500)
    p = cs_amd (A) ;                % use cs_amd on small graphs
else
    [s a b] = cs_nsep (A) ;         % find a node separator
    a = a (cs_nd (A (a,a))) ;       % order A(a,a) recursively
    b = b (cs_nd (A (b,b))) ;       % order A(b,b) recursively
    p = [a b s] ;                   % concatenate to obtain the final ordering
end

author: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [esTheta, esPhi] = mie(radius, frequency, theta, phi, nMax)

% Compute the complex-value scattered electric far field of a perfectly
% conducting sphere using the mie series. Follows the treatment in
% Chapter 3 of
%
% Ruck, et. al. "Radar Cross Section Handbook", Plenum Press, 1970.
%
% The incident electric field is in the -z direction (theta = 0) and is
% theta-polarized. The time-harmonic convention exp(jwt) is assumed, and
% the Green's function is of the form exp(-jkr)/r.
%
% Inputs:
%   radius: Radius of the sphere (meters)
%   frequency: Operating frequency (Hz)
%   theta: Scattered field theta angle (radians)
%   phi: Scattered field phi angle (radians)
%   nMax: Maximum mode for computing Bessel functions
% Outputs:
%   esTheta: Theta-polarized electric field at the given scattering angles
%   esPhi: Phi-polarized electric field at the given scattering angles
%
%   Output electric field values are normalized such that the square of the
%   magnitude is the radar cross section (RCS) in square meters.
%
%   Author: Walton C. Gibson, email: kalla@tripoint.org

% speed of light
c = 299792458.0;

% radian frequency
w = 2.0*pi*frequency;

% wavenumber
k = w/c;

% conversion factor between cartesian and spherical Bessel/Hankel function
s = sqrt(0.5*pi/(k*radius));

% mode numbers
mode = 1:nMax;

% compute spherical bessel, hankel functions
[J(mode)] = besselj(mode + 1/2, k*radius); J = J*s;
[H(mode)] = besselh(mode + 1/2, 2, k*radius); H = H*s;
[J2(mode)] = besselj(mode + 1/2 - 1, k*radius); J2 = J2*s;
[H2(mode)] = besselh(mode + 1/2 - 1, 2, k*radius); H2 = H2*s;

% derivatives of spherical bessel and hankel functions
% Recurrence relationship, Abramowitz and Stegun Page 361
kaJ1P(mode) = (k*radius*J2 - mode .* J );
kaH1P(mode) = (k*radius*H2 - mode .* H );

% Ruck, et. al. (3.2-1)
An = -((i).^mode) .* ( J ./ H ) .* (2*mode + 1) ./ (mode.*(mode + 1));

% Ruck, et. al. (3.2-2), using derivatives of bessel functions
Bn = ((i).^(mode+1)) .* (kaJ1P ./ kaH1P) .* (2*mode + 1) ./ (mode.*(mode + 1));

[esTheta esPhi] = mieScatteredField(An, Bn, theta, phi, frequency);

return

author: Sable

repo: mcbench-benchmarks


text:
function res = hlp_varargin2struct(args, varargin)
% struct = hlp_varargin2struct(varargin, defaults)
% converts a list of name-value pairs into a struct with values assigned to names
%
% In:
%   varargin: cell array of alternating names and values
%             instead of any name-value pair, a struct can be given, which is substituted in-place into the argument list
%   defaults: optional name-value list of defaults; if a default value is '<mandatory>', an error
%             is raised when no value is specified in varargin
%             instead of any name-value pair, a struct can be specified, which is substituted in-place into the argument list
%
% Notes:
%   some parameters may have multiple alternative names, which shall be remapped to the
%   standard name within opts; alternative names are given together with the defaults,
%   by specifying a cell array of names instead of the name, as in the following example:
%   ... ,{'standard_name','alt_name_x','alt_name_y'}, default_value, ...
%
% Out:
%   struct: structure with values assigned to names
%
%                               Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                               2010-04-05

res = struct();
% it is allowed to specify just a struct in place of the varargin clause... (for further convenience)
if isstruct(args)
    args = {args}; end
% splice all substructs into the name-value lists
args = flatten_structs(args);
varargin = flatten_structs(varargin);

% go through the defaults (varargin) and apply them (also, build a remapping table for additional names along the way)
default_name_for = struct();
for i=1:2:length(varargin)
    if iscell(varargin{i})
        % alternative names were specified, remember them
        for j=2:length(varargin{i})
            default_name_for.(varargin{i}{j}) = varargin{i}{1}; end
        varargin{i} = varargin{i}{1};
    end
    try
        res.(varargin{i}) = varargin{i+1};
    catch
        error(['invalid field name specified in defaults: ' exp_fullform(varargin{i})]);
    end
end

% go through the arguments and override
for i=1:2:length(args)
    if isfield(default_name_for,args{i})
        % rewrite alternative names
        args{i} = default_name_for.(args{i}); end
    try
        res.(args{i}) = args{i+1};
    catch
        error(['invalid field name specified in args: ' exp_fullform(args{i})]);
    end
end

% check for missing but mandatory args
for fn = fieldnames(res)'
    if strcmp(res.(fn{1}), '<mandatory>')
        error(['the parameter ' fn{1} ' was unspecified but is mandatory']); end
end

% substitute any structs in place of a name-value pair into the name-value list
function outargs = flatten_structs(inargs)
outargs = {};
k=1;
while k <= length(inargs)
    if isstruct(inargs{k})
        for fn=fieldnames(inargs{k})'
            outargs([end+1 end+2]) = {fn{1} inargs{k}.(fn{1})}; end
        k = k+1;
    else
        outargs([end+1 end+2]) = inargs([k k+1]);
        k = k+2;
    end
end
if mod(length(outargs),2)
    error('a non-matching number of names and values was specified.'); end

repo: REST


text:
function [elon, elat] = eqecl (tjd, icoord, ra, dec)

% this function converts right ascension and declination
% to ecliptic longitude and latitude

% input

%  tjd    = tt julian date of equator, equinox, and ecliptic
%           used for coordinates

%  icoord = coordinate system selection

%               set icoord = 0 for mean equator and equinox of date
%               set icoord = 1 for true equator and equinox of date

%               (ecliptic is always the mean plane)

%  ra     = right ascension in hours, referred to specified
%           equator and equinox of date

%  dec    = declination in degrees, referred to specified
%           equator and equinox of date

% output

%  elon = ecliptic longitude in degrees, referred to specified
%         ecliptic and equinox of date

%  elat = ecliptic latitude in degrees, referred to specified
%         ecliptic and equinox of date

% note:  to convert icrs ra and dec to ecliptic coordinates (mean
% ecliptic and equinox of j2000.0), set tjd = 0.d0 and icoord = 0.
% except for the input to this case, all coordinates are dynamical.

% ported from NOVAS3.0

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

radcon = pi / 180.0d0;

% form position vector in equatorial system from input coordinates

r = ra * 15.0d0 * radcon;

d = dec * radcon;

pos1(1) = cos(d) * cos(r);

pos1(2) = cos(d) * sin(r);

pos1(3) = sin(d);

% convert the vector from equatorial to ecliptic system

pos2 = eqec (tjd, icoord, pos1);

% decompose ecliptic vector into ecliptic longitude and latitude

xyproj = sqrt(pos2(1)^2 + pos2(2)^2);

e = 0.0d0;

if (xyproj > 0.0d0)
    e = atan2 (pos2(2), pos2(1));
end

elon = e / radcon;

if (elon < 0.0d0)
    elon = elon + 360.0d0;
end

e = atan2 (pos2(3), xyproj);

elat = e / radcon;

author: goGPS-Project


text:
function fom = figuremerit(ideal,detected)

    Ni = numel(find(ideal == 1));
    Nd = numel(find(detected == 1));

    [h w c] = size(ideal);

    alpha = 1/9;
    sum = 0.0;
    for r = 1:h
        for c = 1:w
            if (detected(r,c) == 1)
                dist = distance(ideal,detected,r,c)^2;
                sum = sum + (1/(1+alpha*dist));
            end
        end
    end

    fom = (1/max(Ni,Nd))*sum;

function d = distance(ideal,detected,rind,cind)

    [h w c] = size(ideal);

    diff  = 1;
    found = 0;
    d = 10.0^100;

    if (ideal(rind,cind) == detected(rind,cind))
        d     = 0.0;
        found = 1;
    end

    left   = cind - diff;
    right  = cind + diff;
    top    = rind - diff;
    bottom = rind + diff;
    while (~found)

        for c = left:right
            if ((top >= 1) && (c >= 1) && (c <= w))
                if (ideal(top,c) == detected(rind,cind))
                    tmp = sqrt((top-rind)^2 + (c-cind)^2);
                    found = 1;

                    if (tmp < d)
                        d = tmp;
                    end
                end
            end
        end

        for c = left:right
            if ((bottom <= h) && (c >= 1) && (c <= w))
                if (ideal(bottom,c) == detected(rind,cind))
                    tmp = sqrt((bottom-rind)^2 + (c-cind)^2);
                    found = 1;

                    if (tmp < d)
                        d = tmp;
                    end
                end
            end
        end

        for r = top+1:bottom-1
            if ((left >= 1) && (r >= 1) && (r <= h))
                if (ideal(r,left) == detected(rind,cind))
                    tmp = sqrt((r-rind)^2 + (left-cind)^2);
                    found = 1;

                    if (tmp < d)
                        d = tmp;
                    end
                end
            end
        end

        for r = top+1:bottom-1
            if ((right <= w) && (r >= 1) && (r <= h))
                if (ideal(r,right) == detected(rind,cind))
                    tmp = sqrt((r-rind)^2 + (right-cind)^2);
                    found = 1;

                    if (tmp < d)
                        d = tmp;
                    end
                end
            end
        end

        diff = diff + 1;
        left   = cind - diff;
        right  = cind + diff;
        top    = rind - diff;
        bottom = rind + diff;
        if ((left < 1) && (right > w) && (top < 1) && (bottom > h))
            found = 1;
        end
    end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function pass = test_addBreaksAtRoots(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

%%
% Test that pointValues are exactly zero at the new roots.

% Scalar:
f = chebfun(@(x) sin(x)-.5, pref);
g = addBreaksAtRoots(f);
pass(1) = g.pointValues(2) == 0;

% Array-valued:
f = chebfun(@(x) [sin(x), sin(x)-.5], pref);
g = addBreaksAtRoots(f);
pass(2) = g.pointValues(2,1) == 0 && g.pointValues(3,2) == 0;

%% piecewise smooth chebfun: smoothfun + singfun & splitting off.

% define the domain:
dom = [-2 7];
domCheck = [dom(1)+0.1 dom(2)-0.1];

pow1 = -0.5;
pow2 = -1.2;
op = @(x) cos(40*x).*((x-dom(1)).^pow1).*((x-dom(2)).^pow2);
f = chebfun(op, dom, 'exps', [pow1 pow2], 'splitting', 'off');
g = addBreaksAtRoots(f);

% check values:

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

vals_g = feval(g, x);
vals_check = feval(op, x);
err = vals_g - vals_check;

r_exact = (((-25:88)+1/2)*pi/40).';

pass(3) = ( norm(err, inf) < 1e4*eps*norm(vals_check, inf) ) && ...
    ( norm( [dom(1); r_exact; dom(2)] - g.domain.', inf) < ...
    1e4*eps*norm(r_exact, inf) );

%% Tests for functions defined on unbounded domain:

% Functions on [-inf inf]:

% Set the domain:
dom = [-Inf Inf];
domCheck = [-1e2 1e2];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

op = @(x) (1-exp(-x.^2))./x;
f = chebfun(op, dom);
g = addBreaksAtRoots(f);
rExact = 0;

vals_g = feval(g, x);
vals_check = feval(op, x);
err = vals_g - vals_check;
pass(4) = ( norm(err, inf) < 1e2*eps*vscale(f) ) && ...
    ( norm( rExact - g.domain(2:end-1).', inf) < 1e2*eps*vscale(f) );

% Blow-up function:
op = @(x) x.^2.*(1-exp(-x.^2))-2;
f = chebfun(op, dom, 'exps', [2 2]);
g = addBreaksAtRoots(f);
rExact = [-1.4962104914103104707 ; 1.4962104914103104707];

vals_g = feval(g, x);
vals_check = feval(op, x);
err = vals_g - vals_check;
pass(5) = ( norm(err, inf) < 1e7*eps*vscale(f) ) && ...
    ( norm( rExact - g.domain(2:end-1).', inf) < 1e5*eps*vscale(f) );

% Functions on [a inf]:
dom = [0 Inf];
domCheck = [0 100];

% Generate a few random points to use as test values:
x = diff(domCheck) * rand(100, 1) + domCheck(1);

op = @(x) 0.15+sin(10*x)./exp(x);
f = chebfun(op, dom);
g = addBreaksAtRoots(f);
rExact = [0.33529141416564289113;
          0.60061694515161002799;
          0.98375750309184861332;
          1.2042605667187311146;
          1.6619482204330474390;
          1.7760894757659030239];

vals_g = feval(g, x);
vals_check = feval(op, x);
err1 = norm(vals_g - vals_check, inf);
tol1 = 1e2*eps*vscale(f);
err2 = norm( rExact - g.domain(2:end-1).', inf);
tol2 = 1e2*eps*vscale(f);
pass(6) = ( err1 < tol1 ) && ( err2 < tol2 );

% TODO: Add more tests.

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function p = getPosition( o, u )
%getPosition  Get position of graphics object
%
%  p = getPosition(o,u) sets the position of a graphics object o with units
%  u.
%
%  In contrast to getting the Position property directly, this function
%  honors the ActivePositionProperty of axes.

%  Copyright 2009-2016 The MathWorks, Inc.
%  $Revision: 1435 $ $Date: 2016-11-17 17:50:34 +0000 (Thu, 17 Nov 2016) $

% Get position
if ~isempty( findprop( o, 'ActivePositionProperty' ) )
    switch o.ActivePositionProperty
        case 'position'
            q = o.Position;
        case 'outerposition'
            q = o.OuterPosition;
        otherwise
            error( 'uix:InvalidState', ...
                'Unknown value ''%s'' for property ''ActivePositionProperty'' of %s.', ...
                o.ActivePositionProperty, class( o ) )
    end
else
    q = o.Position;
end

% Get units
if ~isempty( findprop( o, 'Units' ) )
    v = o.Units;
else
    v = 'pixels';
end

% Convert
if strcmp( u, v ) % trivial
    p = q;
else % conversion required
    p = hgconvertunits( ancestor( o, 'figure' ), q, v, u, o.Parent );
end

end % getPosition


text:
function [U,s,V] = csvd(A,tst)
%CSVD Compact singular value decomposition.
%
% s = csvd(A)
% [U,s,V] = csvd(A)
% [U,s,V] = csvd(A,'full')
%
% Computes the compact form of the SVD of A:
%    A = U*diag(s)*V',
% where
%    U  is  m-by-min(m,n)
%    s  is  min(m,n)-by-1
%    V  is  n-by-min(m,n).
%
% If a second argument is present, the full U and V are returned.

% Per Christian Hansen, IMM, 06/22/93.

if (nargin==1)
  if (nargout > 1)
    [m,n] = size(A);
    if (m >= n)
      [U,s,V] = svd(full(A),0); s = diag(s);
    else
      [V,s,U] = svd(full(A)',0); s = diag(s);
    end
  else
    U = svd(full(A));
  end
else
  if (nargout > 1)
    [U,s,V] = svd(full(A)); s = diag(s);
  else
    U = svd(full(A));
  end
end

author: indigits

repo: sparse-plex


text:
classdef Zero < CSMaster
    % wrapper class for zero-padded reconstruction
    %
    % (c) Thomas Kuestner
    % ---------------------------------------------------------------------

    properties

    end

    methods
        function obj = Zero(measPara,espresso)
            obj = obj@CSMaster(1, 1, [], [], measPara);
            obj.espresso = espresso;
        end
    end

end

author: thomaskuestner


text:
classdef KRVEA < ALGORITHM
% <multi/many> <real/integer> <expensive>
% Surrogate-assisted RVEA
% alpha ---  2 --- The parameter controlling the rate of change of penalty
% wmax  --- 20 --- Number of generations before updating Kriging models
% mu    ---  5 --- Number of re-evaluated solutions at each generation

%------------------------------- Reference --------------------------------
% T. Chugh, Y. Jin, K. Miettinen, J. Hakanen, and K. Sindhya, A surrogate-
% assisted reference vector guided evolutionary algorithm for
% computationally expensive many-objective optimization, IEEE Transactions
% on Evolutionary Computation, 2018, 22(1): 129-142.
%------------------------------- 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".
%--------------------------------------------------------------------------

% This function is written by Cheng He

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [alpha,wmax,mu] = Algorithm.ParameterSet(2,20,5);

            %% Generate the reference points and population
            [V0,Problem.N] = UniformPoint(Problem.N,Problem.M);
            V     = V0;
            NI    = 11*Problem.D-1;
            P     = UniformPoint(NI,Problem.D,'Latin');
            A2    = Problem.Evaluation(repmat(Problem.upper-Problem.lower,NI,1).*P+repmat(Problem.lower,NI,1));
            A1    = A2;
            THETA = 5.*ones(Problem.M,Problem.D);
            Model = cell(1,Problem.M);

            %% Optimization
            while Algorithm.NotTerminated(A2)
                % Refresh the model and generate promising solutions
                A1Dec = A1.decs;
                A1Obj = A1.objs;
                for i = 1 : Problem.M
                    % The parameter 'regpoly1' refers to one-order polynomial
                    % function, and 'regpoly0' refers to constant function. The
                    % former function has better fitting performance but lower
                    % efficiency than the latter one
                    dmodel     = dacefit(A1Dec,A1Obj(:,i),'regpoly1','corrgauss',THETA(i,:),1e-5.*ones(1,Problem.D),100.*ones(1,Problem.D));
                    Model{i}   = dmodel;
                    THETA(i,:) = dmodel.theta;
                end
                PopDec = A1Dec;
                w      = 1;
                while w <= wmax
                    drawnow('limitrate');
                    OffDec = OperatorGA(Problem,PopDec);
                    PopDec = [PopDec;OffDec];
                    [N,~]  = size(PopDec);
                    PopObj = zeros(N,Problem.M);
                    MSE    = zeros(N,Problem.M);
                    for i = 1: N
                        for j = 1 : Problem.M
                            [PopObj(i,j),~,MSE(i,j)] = predictor(PopDec(i,:),Model{j});
                        end
                    end
                    index  = KEnvironmentalSelection(PopObj,V,(w/wmax)^alpha);
                    PopDec = PopDec(index,:);
                    PopObj = PopObj(index,:);
                    % Adapt referece vectors
                    if ~mod(w,ceil(wmax*0.1))
                        V(1:Problem.N,:) = V0.*repmat(max(PopObj,[],1)-min(PopObj,[],1),size(V0,1),1);
                    end
                    w = w + 1;
                end

                % Select mu solutions for re-evaluation
                [NumVf,~] = NoActive(A1Obj,V0);
                PopNew    = KrigingSelect(PopDec,PopObj,MSE(index,:),V,V0,NumVf,0.05*Problem.N,mu,(w/wmax)^alpha);
                New       = Problem.Evaluation(PopNew);
                A2        = [A2,New];
                A1        = UpdataArchive(A1,New,V,mu,NI);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function phi = CharacteristicFunctionLib(model,u,lnS,T,r,d,varargin)
%---------------------------------------------------------
% Characteristic Function Library of the following models:
%---------------------------------------------------------
% Black Scholes
% Merton Jump Diffusion
% Heston Stochastic Volatility Model
% Bates Stochastic Volatility / Jump Diffusion Model
% Variance Gamma
% Normal Inverse Gaussian
% Meixner
% Generalized Hyperbolic
% CGMY
% Variance Gamma with Gamma Ornstein Uhlenbeck clock
% Variance Gamma with CIR clock
%---------------------------------------------------------

optAlfaCalculation = true;

ME1 = MException('VerifyInput:InvalidNrOfArguments',...
    'Invalid number of Input arguments');
ME2 = MException('VerifyInput:InvalidModel',...
    'Undefined Model');

if strcmp(model,'BlackScholes')
    if nargin == 7
        funobj = @Black_characteristicFn;
    else
        throw(ME1)
    end
elseif strcmp(model,'Heston')
    if nargin == 11
        funobj = @Heston_characteristicFn;
    else
        throw(ME1)
    end
elseif strcmp(model,'Merton')
    if nargin == 10
        funobj = @Merton_characteristicFn;
    else
        throw(ME1)
    end
elseif strcmp(model,'Bates')
    if nargin == 14
        funobj = @Bates_characteristicFn;
    else
        throw(ME1)
    end
elseif strcmp(model,'VarianceGamma')
    if nargin == 9
        funobj = @VG_characteristicFn;
    else
        throw(ME1)
    end
elseif strcmp(model,'NIG')
    if nargin == 10
        funobj = @NIG_characteristicFn;
    else
        throw(ME1)
    end
elseif strcmp(model,'CGMY')
    if nargin == 10
        funobj = @CGMY_characteristicFn;
    else
        throw(ME1)
    end
elseif strcmp(model,'Meixner')
    if nargin == 9
        funobj = @Meixner_characteristicFn;
    else
        throw(ME1)
    end
elseif strcmp(model,'GH')
    if nargin == 10
        funobj = @GH_characteristicFn;
    else
        throw(ME1)
    end
elseif strcmp(model,'IntegratedCIR')
    if nargin == 13
        funobj = @VarianceGammaCIR_characteristicFn;
    else
        throw(ME1)
    end
elseif strcmp(model,'VarianceGammaOU')
    if nargin == 12
        funobj = @VarianceGammaOU_characteristicFn;
    else
        throw(ME1)
    end
else
    throw(ME2)
end

fval = feval(funobj,u,lnS,T,r,d,varargin{:});

if optAlfaCalculation == true
    phi = fval;
else
    phi = exp(fval);
end

end

%% Explicit Implementation of the characteristic Functions E[exp(iu*lnS_T)]
%-----------------------------------------------------------------------

% Black Scholes
function phi = Black_characteristicFn(u,lnS,T,r,d,sigma)
    %phi = exp(i*u*(lnS+(r-d-0.5*sigma*sigma)*T) - 0.5*sigma*sigma*u.*u*T);
    phi = 1i*u*(lnS+(r-d-0.5*sigma*sigma)*T) - 0.5*sigma*sigma*u.*u*T;
end

% Merton Jump Diffusion
function phi = Merton_characteristicFn(u,lnS,T,r,d,sigma,a,b,lambda)
    phi = Black_characteristicFn(u,lnS,T,r,d,sigma) + LogNormalJump_characteristicFn(u,a,b,lambda,T);
end

% Heston
function phi = Heston_characteristicFn(u,lnS,T,r,d,V0,theta,kappa,omega,rho)

alfa = -0.5*(u.*u + u*1i);
beta = kappa - rho*omega*u*1i;

omega2 = omega * omega;
gamma = 0.5 * omega2;

D = sqrt(beta .* beta - 4.0 * alfa .* gamma);

bD = beta - D;
eDt = exp(- D * T);

G = bD ./ (beta + D);
B = (bD ./ omega2) .* ((1.0 - eDt) ./ (1.0 - G .* eDt));
psi = (G .* eDt - 1.0) ./(G - 1.0);
A = ((kappa * theta) / (omega2)) * (bD * T - 2.0 * log(psi));

phi = A + B*V0 + 1i*u*(lnS+(r-d)*T);

end

% Bates
function phi = Bates_characteristicFn(u,lnS,T,r,d,V0,theta,kappa,omega,rho,a,b,lambda)
    phiHes = Heston_characteristicFn(u,lnS,T,r,d,V0,theta,kappa,omega,rho);
    phi = phiHes + LogNormalJump_characteristicFn(u,a,b,lambda,T);
end

% LogNormalJump for Merton and Bates
function phiJump = LogNormalJump_characteristicFn(u,a,b,lambda,T)
    %phiJump = exp(lambda*T*(-a*u*i + (exp(u*i*log(1.0+a)+0.5*b*b*u*i.*(u*i-1.0))-1.0)));
    phiJump = lambda*T*(-a*u*1i + (exp(u*1i*log(1.0+a)+0.5*b*b*u*1i.*(u*1i-1.0))-1.0));
end

% Variance Gamma
function phi = VG_characteristicFn(u,lnS,T,r,d,sigma,nu,theta)
    omega = (1/nu)*( log(1-theta*nu-sigma*sigma*nu/2) );
    tmp = 1 - 1i * theta * nu * u + 0.5 * sigma * sigma * u .* u * nu;
    %tmp = tmp.^(-T / nu);
    %phi = exp( i * u * (lnS + (r + omega - d) * T )) .* tmp;
    phi = 1i * u * (lnS + (r + omega - d) * T ) - T*log(tmp)/nu;

end

% Normal Inverse Gaussian
function phi = NIG_characteristicFn(u,lnS,T,r,d,alfa,beta,mu,delta)
    m = delta*(sqrt(alfa*alfa-(beta+1)^2)-sqrt(alfa*alfa-beta*beta));
    %tmp = exp(i*u*mu*T-delta*T*(sqrt(alfa*alfa-(beta+i*u).^2)-sqrt(alfa*alfa-beta*beta)));
    %phi = exp( i*u*(lnS + (r-d+m)*T)).*tmp;
    tmp = 1i*u*mu*T-delta*T*(sqrt(alfa*alfa-(beta+1i*u).^2)-sqrt(alfa*alfa-beta*beta));
    phi = 1i*u*(lnS + (r-d+m)*T) + tmp;
end

% Meixner
function phi = Meixner_characteristicFn(u,lnS,T,r,d,alfa,beta,delta)
    m = -2*delta*(log(cos(0.5*beta)) - log(cos((alfa+beta)/2)));
    tmp = (cos(0.5*beta)./cosh(0.5*(alfa*u-1i*beta))).^(2*T*delta);
    %phi = exp( i*u*(lnS + (r-d+m)*T)).*tmp;
    phi = 1i*u*(lnS + (r-d+m)*T) + log(tmp);
end

% Generalized Hyperbolic
function phi = GH_characteristicFn(u,lnS,T,r,d,alfa,beta,delta,nu)
    arg1 = alfa*alfa-beta*beta;
    arg2 = arg1-2*1i*u*beta + u.*u;
    argm = arg1-2*beta-1;
    m = -log((arg1./argm).^(0.5*nu).*besselk(nu,delta*sqrt(argm))./besselk(nu,delta*sqrt(arg1)));
    tmp = (arg1./arg2).^(0.5*nu).*besselk(nu,delta*sqrt(arg2))./besselk(nu,delta*sqrt(arg1));
    %phi = exp( i*u*(lnS + (r-d+m)*T)).*tmp.^T;
    phi = 1i*u*(lnS + (r-d+m)*T) + log(tmp).*T;
end

% CGMY
function phi = CGMY_characteristicFn(u,lnS,T,r,d,C,G,M,Y)
    m = -C*gamma(-Y)*((M-1)^Y-M^Y+(G+1)^Y-G^Y);
    %tmp = exp(C*T*gamma(-Y)*((M-i*u).^Y-M^Y+(G+i*u).^Y-G^Y));
    %phi = exp( i*u*(lnS + (r-d+m)*T)).*tmp;
    tmp = C*T*gamma(-Y)*((M-1i*u).^Y-M^Y+(G+1i*u).^Y-G^Y);
    phi = 1i*u*(lnS + (r-d+m)*T) + tmp;
end

%function phi = IntegratedCIR_characteristicFn(u,lnS,T,r,d,sigma,nu,theta,kappa,eta,lambda,y0)

%v1 = i*log(1 - i * theta * nu * u + 0.5 * sigma * sigma * u .* u * nu)/nu;
%v2 = i*log(1 - i * theta * nu *(-i) + 0.5 * sigma * sigma * (-i) .* (-i) * nu)/nu;

function phi = VarianceGammaCIR_characteristicFn(u,lnS,T,r,d,C,G,M,kappa,eta,lambda,y0)

v1 = -1i*C*(log(G*M)-log(G*M+(M-G)*1i*u + u.*u));
v2 = -1i*C*(log(G*M)-log(G*M+(M-G)*1i*(-1i) + (-1i).*(-1i)));

gamma1 = sqrt(kappa^2 - 2*lambda^2*1i*v1);
gamma2 = sqrt(kappa^2 - 2*lambda^2*1i*v2);
phi1 = kappa^2*eta*T/lambda^2 + 2*y0*1i*v1 ./ (kappa + gamma1.*coth(0.5*gamma1*T))...
    - 2*kappa*eta/lambda^2*log(cosh(0.5*gamma1*T) + kappa*sinh(0.5*gamma1*T)./gamma1);
phi2 = kappa^2*eta*T/lambda^2 + 2*y0*1i*v2 / (kappa + gamma2*coth(0.5*gamma2*T))...
    - 2*kappa*eta/lambda^2*log(cosh(0.5*gamma2*T) + kappa*sinh(0.5*gamma2*T)/gamma2);

%gam = sqrt(kappa^2-2*lambda^2);
m = 0;%- (kappa^2*eta*T/lambda^2 + 2*y0/(kappa+gam*coth(gam*T/2))) + 2*kappa*eta/lambda^2*log(cosh(gam*T/2)+kappa*sinh(gam*T/2)/gam);
phi = 1i*u*(lnS + (r-d+m)*T) + phi1 - 1i*u*phi2;%0.5*(log(phi1.^2) - i*u*log(phi2.^2));

end

function phi = VarianceGammaOU_characteristicFn(u,lnS,T,r,d,C,G,M,lambda,a,b)
    psiX1 = (-1i)*log((G*M./(G*M+(M-G)*1i*u+u.*u)).^C);
    psiX2 = (-1i)*log((G*M/(G*M+(M-G)-1)).^C);
    phiCIR1 = 1i*psiX1*lambda^(-1)*(1-exp(-lambda*T)) + lambda*a./(1i*psiX1-lambda*b).*(b*log(b./(b-1i*psiX1*lambda^(-1)*(1-exp(-lambda*T))))-1i*psiX1*T);
    phiCIR2 = 1i*psiX2*lambda^(-1)*(1-exp(-lambda*T)) + lambda*a./(1i*psiX2-lambda*b).*(b*log(b./(b-1i*psiX2*lambda^(-1)*(1-exp(-lambda*T))))-1i*psiX2*T);
    phi = 1i*u*(lnS + (r-d)*T)+phiCIR1 - 1i*u.*phiCIR2;
end


text:
function [coherogram,phase,t,f] = MTCoherogram(lfp1,lfp2,varargin)

%MTCoherogram - Compute LFP coherogram by multi-taper estimation.
%
%  USAGE
%
%    [coherogram,phase,t,f] = MTCoherogram(lfp1,lfp2,<options>)
%
%    lfp1,lfp2      wide-band LFPs (one channel each).
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'frequency'   sampling rate (in Hz) (default = from timestamps if
%                   available, otherwise 1250Hz)
%     'range'       frequency range (in Hz) (default = all)
%     'window'      duration (in s) of the time window (default = 5)
%     'overlap'     overlap between successive windows (default = window/2)
%     'step'        step between successive windows (default = window/2)
%     'tapers'      relative resolution and order of the tapers [NW K]
%                   (default = [3 5])
%     'pad'         FFT padding (see help for cohgramc) (default = 0)
%     'show'        plot results (default = 'off')
%     'cutoffs'     cutoff values for color plot (default = [0 1])
%    =========================================================================
%
%  NOTES
%
%    The LFP can be provided either as a time stamped matrix (list of time-voltage
%    pairs), or as a voltage vector - in which case the frequency must be specified.
%
%    The time displacement between successive short time coherences can be supplied
%    either as a 'step' (explicit time difference) or as an 'overlap' (between
%    successive time windows).
%
%  OUTPUT
%
%    coherogram     coherogram magnitude
%    phase          coherogram phase
%    t              time bins
%    f              frequency bins
%
%  DEPENDENCIES
%
%    This function requires the chronux toolbox.
%
%  SEE
%
%    See also MTCoherence, MTSpectrum, MTSpectrogram, PlotColorMap.

% Copyright (C) 2010-2014 by Michaël Zugaro
%
% 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.
warning('This function has been deprecated, and will be removed in the future. Try using bz_MTCoherogram.m')
% Make sure chronux is installed and functional
CheckChronux('cohgramc');

% Defaults
f = 1250;
frequency = [];
window = 5;
range = [];
overlap = [];
step = [];
show = 'off';
tapers = [3 5];
pad = 0;
cutoffs = [0 1];

% Check number of parameters
if nargin < 2 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help MTCoherogram'' for details).');
end

% Check parameter sizes
if size(lfp1,2) ~= 1 && size(lfp1,2) ~= 2,
	error('Parameter ''lfp1'' is not a vector or a Nx2 matrix (type ''help MTCoherogram'' for details).');
end
if size(lfp2,2) ~= 1 && size(lfp2,2) ~= 2,
	error('Parameter ''lfp2'' is not a vector or a Nx2 matrix (type ''help MTCoherogram'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+2) ' is not a property (type ''help MTCoherogram'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'frequency',
			frequency = varargin{i+1};
			if ~isdscalar(frequency,'>0'),
				error('Incorrect value for property ''frequency'' (type ''help MTCoherogram'' for details).');
			end
		case 'range',
			range = varargin{i+1};
			if ~isdvector(range,'#2','<','>=0'),
				error('Incorrect value for property ''range'' (type ''help MTCoherogram'' for details).');
			end
		case 'window',
			window = varargin{i+1};
			if ~isdscalar(window,'>0'),
				error('Incorrect value for property ''window'' (type ''help MTCoherogram'' for details).');
			end
		case 'overlap',
			overlap = varargin{i+1};
			if ~isdscalar(overlap,'>0'),
				error('Incorrect value for property ''overlap'' (type ''help MTCoherogram'' for details).');
			end
		case 'step',
			step = varargin{i+1};
			if ~isdscalar(step,'>0'),
				error('Incorrect value for property ''step'' (type ''help MTCoherogram'' for details).');
			end
		case 'tapers',
			tapers = varargin{i+1};
			if ~isivector(tapers,'#2','>0'),
				error('Incorrect value for property ''tapers'' (type ''help MTCoherogram'' for details).');
			end
		case 'pad',
			pad = varargin{i+1};
			if ~isiscalar(pad,'>-1'),
				error('Incorrect value for property ''pad'' (type ''help MTCoherogram'' for details).');
			end
		case 'show',
			show = varargin{i+1};
			if ~isstring_FMAT(show,'on','off'),
				error('Incorrect value for property ''show'' (type ''help MTCoherogram'' for details).');
			end
		case 'cutoffs',
			cutoffs = varargin{i+1};
			if ~isdvector(cutoffs,'#2','>=0','<'),
				error('Incorrect value for property ''cutoffs'' (type ''help MTCoherogram'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help MTCoherogram'' for details).']);
	end
end

% Determine LFP frequency
if isempty(frequency),
	if size(lfp1,2) == 2,
		frequency = 1/median(diff(lfp1(:,1)));
	else
		frequency = f;
	end
end

% Determine step/overlap
if isempty(step),
	if isempty(overlap),
		overlap = window/2;
	end
else
	if isempty(overlap),
		overlap = window-step;
	elseif overlap ~= window-step,
		error('Incompatible ''step'' and ''overlap'' parameters (type ''help MTCoherogram'' for details).');
	end
end

% Compute and plot coherogram
parameters.Fs = frequency;
if ~isempty(range), parameters.fpass = range; end
parameters.tapers = tapers;
parameters.pad = pad;
[coherogram,phase,~,~,~,t,f] = cohgramc(lfp1(:,2),lfp2(:,2),[window window-overlap],parameters);
t = t'+lfp1(1,1);
f = f';
coherogram = coherogram';
% coherogram = permute(coherogram,[2 1 3]);  % Previous code by Gabrielle Girardeau, keep it around just in case
phase = phase';
if strcmp(lower(show),'on'),
  figure;hold on;
  subplot(2,1,1);
  PlotColorMap(coherogram,'x',t,'y',f,'cutoffs',cutoffs,'newfig','off');
  xlabel('Time (s)');
  ylabel('Frequency (Hz)');
  title('Coherogram Amplitude');
  subplot(2,1,2);
  PlotColorMap(phase,'x',t,'y',f,'cutoffs',[-pi pi],'newfig','off');
  xlabel('Time (s)');
  ylabel('Frequency (Hz)');
  title('Coherogram Phase');
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function plotci2()
    %  plot a circle containing ni events
    %  around each grid point
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    st = 2;
    [X,Y] = meshgrid(gx,gy);
    [m,n]= size(r);
    set(gca,'NextPlot','add')
    x = -pi-0.1:0.1:pi;
    for i = 1:st:m
        for k = 1:st:n
            if r(i,k) <= ZG.tresh_km
                plot(X(i,k)+r(i,k)*sin(x)/(cosd(ya0)*111),Y(i,k)+r(i,k)*cos(x)/(cosd(ya0)*111) ,'k')
                plot(X(i,k),Y(i,k),'+k')
            end
        end
    end

end

author: CelsoReyes


text:
classdef ShFunc_Volume < ShapeFunctional

    properties (Access = public)
        geometricVolume
    end

    methods (Access = public)

        function obj = ShFunc_Volume(cParams)
            cParams.filterParams.quadratureOrder = 'LINEAR';
            obj.init(cParams);
            obj.geometricVolume = sum(obj.dvolu(:));
        end

        function computeFunctionAndGradient(obj)
            obj.nVariables = obj.designVariable.nVariables;
            obj.updateHomogenizedMaterialProperties();
            obj.computeFunction();
            obj.computeGradient();
        end

        function computeFunction(obj)
            density = obj.homogenizedVariablesComputer.rho;
            obj.computeFunctionFromDensity(density);
        end

        function computeFunctionFromDensity(obj,dens)
            densV = dens;
            volume = (sum(obj.dvolu,2)'*mean(densV,2));
            volume = volume/(obj.geometricVolume);
            obj.value = volume;
        end

        function v = getVariablesToPlot(obj)
            v{1} = obj.value*obj.value0;
        end

        function t = getTitlesToPlot(obj)
            t{1} = 'Volume non scaled';
        end

    end

    methods (Access = private)

        function computeGradient(obj)
            drho = obj.homogenizedVariablesComputer.drho;
            g = drho;
            gf = zeros(size(obj.Msmooth,1),obj.nVariables);
            for ivar = 1:obj.nVariables
                gs = g{ivar}/obj.geometricVolume;
                gf(:,ivar) = obj.filter.getP1fromP0((gs));
            end
            g = obj.Msmooth*gf;
            %g = gf;
            obj.gradient = g(:);
            %obj.gradient = ones(size(g(:)))/obj.geometricVolume;
        end

        function updateHomogenizedMaterialProperties(obj)
            nx = length(obj.designVariable.value)/obj.designVariable.nVariables;
            x  = obj.designVariable.value;
            xf = cell(obj.designVariable.nVariables,1);
            for ivar = 1:obj.nVariables
                i0 = nx*(ivar-1) + 1;
                iF = nx*ivar;
                xs = x(i0:iF);
                xf{ivar} = obj.filter.getP0fromP1(xs);
            end
            obj.homogenizedVariablesComputer.computeDensity(xf);
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef nme_bus_ld_acp_node_test < mp.nme_bus_nld_acp_node_test

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        function name = name(obj)
            name = 'bus_ld';
        end

        function np = np(obj)
            np = 1;     %% this is a 1 port element
        end

        function obj = build_params(obj, nm, dm)
            build_params@mp.nme_bus_nld_acp_node_test(obj, nm, dm); %% call parent

            dme = obj.data_model_element(dm);

            %% bus shunts
            Ysh = dme.gs + 1j * dme.bs;             %% shunt admittances
            obj.Y = sparse(1:obj.nk, 1:obj.nk, Ysh, obj.nk, obj.nk);

            %% constant power loads
            obj.s = dme.pd + 1j * dme.qd;           %% complex power demand
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [t_stat,B,SE,F_stat] = gen_bf_ttest(X,Y,c,U)
% Run a t-test

if nargin<4,
    U=[];
end;
if isempty(U),
    U=eye(size(Y,2));
end;

X0  = X - X*c*pinv(c);  %% make sure X0 is orthogonal to X
Xred   = full(X*c); %% reduced design matrix
X0  = spm_svd(X0); %% X0 is null space i.e. everything that is happening in other columns of X

%   ==========================================================================
% remove null space of contrast
%--------------------------------------------------------------------------
Y     = Y - X0*(X0'*Y); %% eg remove DC level or drift terms from all of Y
Xred     = Xred - X0*(X0'*Xred);

P     = pinv(Xred);

[n,b] = size(Xred);
[n,m] = size(Y); %% n is number of epochs, m is number of features
b     = rank(Xred);
h     = min(b,m); %% either number of features or rank of X

Ym=mean(Y,2);
B  = pinv(Xred)*Ym;
RSS   = sum((Ym - Xred*B).^2);
MRSS  = RSS / (n-b);
SE    = sqrt(MRSS*(pinv(Xred'*Xred)));
t_stat=B./SE;
F_stat=(B./SE).^2;

end

author: spm


text:
function z = generateSimulatedBouncingPath(iterationNum, roundTimes)
% This function generates the Z-coordinate path for the center of heart.

% by Xin Zhao
% Feb 14, 2008

% define the z path of the center of ball for rountTimes of bouncing along z-axis
if nargin < 2
    roundTimes = 3;
end

if nargin < 1
    iterationNum = 40;
end

radius = 1;
deltaH = 1.5; % the highest point comapred with radius
deltaL = -0.3;% the lowest point comapred with radius

rangeY = [-2, 1];

% gravity accelerator constant
% Note: this is smaller to the real 9.8 value, but this makes the bouncing slower
g = 5;

% define the accelerator for the z < 0
% k is similar to spring coefficient, but it's more complicated than that.
% k just needs to satisfy one condition here, i.e. it will stop at location defined by
% deltaL*radius
k = g * deltaH/(-deltaL);

% start point
startZ = deltaH*radius;

% speed when z =0
V0 = sqrt(2*startZ*g);

% calcualte the total time for rountTimes of bouncing
% time for z > 0
plusZtime = sqrt(2*deltaH*radius/g);
% time for z < 0
minusZtime = sqrt(-2*deltaL*radius/k);

% half period
qPeriod = plusZtime + minusZtime;

% define totalTime
qPeriodTime = linspace(0, qPeriod, iterationNum);

% z path when z > 0
zPlus = startZ - 0.5 * g * qPeriodTime(qPeriodTime<plusZtime).^2;

% z path when z < 0
zMinus = -V0 * (qPeriodTime(qPeriodTime >=plusZtime)-plusZtime) + 0.5 * k* (qPeriodTime(qPeriodTime >=plusZtime)-plusZtime).^2;

z = [zPlus, zMinus, fliplr(zMinus), fliplr(zPlus)];
% the last zPlus is for the stop action
z = [repmat(z, 1, roundTimes), zPlus];
y = linspace(rangeY(1), rangeY(2), numel(z));


text:
%%
%% Run this script in Matlab command window
%%

   function Installmex
   curdir = pwd;
   fprintf(' current directory is:  %s\n',curdir);
   fsp = filesep;
%%
%%
%%
   src = [curdir,fsp,'Solver',fsp,'Mexfun'];
   eval(['cd ','Solver',fsp,'Mexfun']);
   fprintf ('\n Now compiling the mexFunctions in:\n');
   fprintf (' %s\n',src);

   fname{1} = 'mexProd2';
   fname{2} = 'mexProd2nz';
   fname{3} = 'mexinprod';
   fname{4} = 'mexmat';
   fname{5} = 'mexsmat';
   fname{6} = 'mexsvec';
   fname{7} = 'mexschur';
   fname{8} = 'mexqops';
   fname{9} = 'mexexpand';
   fname{10} = 'mexskron';
   fname{11} = 'mexnnz';
   fname{12} = 'mexschurfun';
   fname{13} = 'mexMatvec';
   fname{14} = 'mextriang';
   fname{15} = 'mextriangsp';

   if exist('OCTAVE_VERSION'),
      mexcmd = 'mex ';
   else
      mexcmd = 'mex  -O -largeArrayDims  ';
   end
   for k = 1:length(fname)
       cmd = [mexcmd,fname{k},'.c'];
       disp( cmd );
       eval( cmd );
   end
   cd ..
   cd ..


text:
function y = bin2polar(x)
y = 2*x-1;


text:
function demo_svg_water
% Create a new figure and remember the figure handle
fig = figure;
% Set default font and font size
set(fig, 'DefaultAxesFontName', 'Arial')
set(fig, 'DefaultAxesFontSize', 16)
% This is the data [year water_consumption]
data = [
1980	230
1982	265
1984	265
1986	266
1988	262
1990	259
1992	277
1994	247
1996	238
1998	250
2000	250
2002	235
2004	233
2006	228
];
% Let's plot the data
hold on
s = bar(data(:, 1), data(:, 2));
% We do not want to apply the filter to the edge of the bars. Therefore, we
% plot them in addition on top
h = bar(data(:, 1), data(:, 2));
set(h, 'FaceColor', 'none');
set(h, 'EdgeColor', 'black');
% Add the text inside each bar with correct rotation and alignment
t = text(data(:, 1), data(:, 2) - 10, num2str(data(:, 2)));
set(t, 'Rotation', -90);
set(t, 'FontSize', 16);
set(t, 'FontName', 'Arial');
set(t, 'VerticalAlignment', 'middle');
set(t, 'FontWeight', 'bold');
% Adding labels and tick marks
axis([1979 2008 0 400])
set(gca, 'XTick', data(:, 1));
set(gca, 'XTickLabel', num2str(data(:, 1)));
title('Water Consumption per Person and Year')
ylabel('Water consumption [l]')
grid on
box on
set(gca, 'Position', [0.13 0.17 0.80 0.72]);
% Now we add the filters
% The bounding box with extension axes makes sure that we cover the full
% axis region with the background images. Due to the shadow we have to
% define an overlap region of 12px. Otherwise, distortions at the border of
% the axis reagion may be visible.
svgBoundingBox(s, 'axes', 12, 'off');
% Now we add the image. To keep the aspect ratio we use 'xMidYMid slice'.
% This setting centers the picture for both x and y. The kewyord 'slice'
% makes sure that the picture is scaled so that the full axis region is
% covered.
svgImage(s, 'water_stones.jpg', 'xMidYMid slice', 'pic');
% Create a composite between bars and picture by union.
svgComposite(s, 'pic', 'SourceGraphic', 'atop', 'obj');
% To create the shadow we use the bars and use a blur filter
svgGaussianBlur(s, 'SourceAlpha', 5, 'blur');
% The blurred bars are now shifted by [10 10].
svgOffset(s, 'blur', [10 10], 'shade');
% Combine the shadow and 'picture bars' by applying the later on top of the
% shadow.
svgComposite(s, 'obj', 'shade', 'over', 'front');
% For the background we use the same picture but make it brighter
svgComposite(s, 'pic', 'pic', 'arithmetic', 'background', [0 0.2 0 0]);
% As last step, we combine the foreground and background
svgComposite(s, 'front', 'background', 'over', 'final');
% Beautify the tick mark labels by rotating them by 90 deg
setting.svg.XTickLabelAngle = -90;
set(gca, 'UserData', setting);
% Finally, we save the result into a SVG file
plot2svg('demo_svg_water.svg')

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function s = CubicTimeScaling(Tf, t)
% *** CHAPTER 9: TRAJECTORY GENERATION ***
% Takes Tf: Total time of the motion in seconds from rest to rest,
%       t: The current time t satisfying 0 < t < Tf.
% Returns s: The path parameter s(t) corresponding to a third-order
%            polynomial motion that begins and ends at zero velocity.
% Example Input:
%
% clear; clc;
% Tf = 2;
% t = 0.6;
% s = CubicTimeScaling(Tf,t)
%
% Output:
% s =
%    0.2160

s = 3 * (t / Tf) ^ 2 - 2 * (t / Tf) ^ 3;
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function t_auction_minopf(quiet)
%T_AUCTION_MINOPF  Tests for code in auction.m, using MINOPF solver.

%   MATPOWER
%   Copyright (c) 2004-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

n_tests = 183;

t_begin(n_tests, quiet);

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;

if ~have_feature('smartmarket')
    t_skip(n_tests, 'smartmarket code not available');
elseif ~have_feature('minopf')
    t_skip(n_tests, 't_auction_minopf requires MINOPF');
else
    mpopt = mpoption('opf.ac.solver', 'MINOPF', 'out.lim.all', 1, 'out.branch', 0, 'out.sys_sum', 0, 'out.all', 0, 'verbose', 0);
    q = [
        12 24 24;
        12 24 24;
        12 24 24;
        12 24 24;
        12 24 24;
        12 24 24;
        10 10 10;
        10 10 10;
        10 10 10;
    ];

    %%-----  one offer block marginal @ $50  -----
    p = [
        20 50 60;
        20 40 70;
        20 42 80;
        20 44 90;
        20 46 75;
        20 48 60;
        100 70 60;
        100 50 20;
        100 60 50;
    ];

    t = 'one marginal offer @ $50, auction_type = 5';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1150, 100, [], [], mpopt);
    cq5 = cq;
    cp5 = cp;
    i2e = bus(:, BUS_I);
    e2i = sparse(max(i2e), 1);
    e2i(i2e) = (1:size(bus, 1))';
    G = find( ~isload(gen) );   %% real generators
    L = find(  isload(gen) );   %% dispatchable loads
    Gbus = e2i(gen(G,GEN_BUS));
    Lbus = e2i(gen(L,GEN_BUS));
    Qfudge =  zeros(size(p));
    Qfudge(L,:) = diag(gen(L,QG) ./ gen(L,PG) .* bus(Lbus, LAM_Q)) * ones(size(p(L,:)));

    t_is( cq(G(1),2:3), [23.32 0], 2, t );
    t_is( cp(G(1),:), 50, 4, t );
    t_is( cq(L(2),1:2), [10 0], 2, t );
    t_is( cp(L(2),:), 54.0312, 4, t );
    t_is( cp(G,1), bus(Gbus, LAM_P), 8, [t ' : gen prices'] );
    t_is( cp(L,1), bus(Lbus, LAM_P) + Qfudge(L,1), 8, [t ' : load prices'] );

    lao_X = p(G(1),2)/bus(Gbus(1), LAM_P);
    fro_X = p(G(6),3)/bus(Gbus(6), LAM_P);
    lab_X = p(L(3),2)/(bus(Lbus(3), LAM_P) + Qfudge(L(3),1));
    frb_X = p(L(2),2)/(bus(Lbus(2), LAM_P) + Qfudge(L(2),1));

    t_is( lao_X, 1, 4, 'lao_X');
    t_is( fro_X, 1.1324, 4, 'fro_X');
    t_is( lab_X, 1.0787, 4, 'lab_X');
    t_is( frb_X, 0.9254, 4, 'frb_X');

    t = 'one marginal offer @ $50, auction_type = 1';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1110, 100, [], [], mpopt);
    cp1 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5, 8, [t ' : prices'] );

    t = 'one marginal offer @ $50, auction_type = 2';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1120, 100, [], [], mpopt);
    cp2 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G,:), cp5(G,:)*fro_X, 8, [t ' : gen prices'] );
    t_is( cp(L(1:2),:), cp5(L(1:2),:)*fro_X, 8, [t ' : load 1,2 prices'] );
    t_is( cp(L(3),:), 60, 5, [t ' : load 3 price'] );   %% clipped by accepted bid

    t = 'one marginal offer @ $50, auction_type = 3';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1130, 100, [], [], mpopt);
    cp3 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5*lab_X, 8, [t ' : prices'] );

    t = 'one marginal offer @ $50, auction_type = 4';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1140, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G(1),:), p(G(1),2), 8, [t ' : gen 1 price'] );
    t_is( cp(G(2:6),:), cp5(G(2:6),:)*frb_X, 8, [t ' : gen 2-6 prices'] );
    t_is( cp(L,:), cp5(L,:)*frb_X, 8, [t ' : load prices'] );

    t = 'one marginal offer @ $50, auction_type = 6';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1160, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp3, 8, [t ' : prices'] );
    p2 = p;
    p2(L,:) = [ 100 100 100;
                100   0   0;
                100 100   0 ];
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1160, 100, [], [], mpopt);
    t_is( cq, cq5, 5, [t ' : quantities'] );
    t_is( cp(G,:), cp5(G,:)*fro_X, 4, [t ' : gen prices'] );
    t_is( cp(L,:), cp5(L,:)*fro_X, 4, [t ' : load prices'] ); %% load 3 not clipped as in FRO

    t = 'one marginal offer @ $50, auction_type = 7';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1170, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5 * (lao_X+lab_X)/2, 8, [t ' : prices'] );
    t_is( cp, (cp1 + cp3) / 2, 8, [t ' : prices'] );

    t = 'one marginal offer @ $50, auction_type = 8';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1180, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G,:), cp1(G,:), 8, [t ' : gen prices'] );
    t_is( cp(L,:), cp3(L,:), 8, [t ' : load prices'] );

    t = 'one marginal offer @ $50, auction_type = 0';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1100, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, p, 8, [t ' : prices'] );

    %%-----  one bid block marginal @ $55  -----
    p(L(2),2) = 55;
    t = 'one marginal bid @ $55, auction_type = 5';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1150, 100, [], [], mpopt);
    cq5 = cq;
    cp5 = cp;
    Qfudge =  zeros(size(p));
    Qfudge(L,:) = diag(gen(L,QG) ./ gen(L,PG) .* bus(Lbus, LAM_Q)) * ones(size(p(L,:)));

    t_is( cq(G(1),2:3), [24 0], 2, t );
    t_is( cp(G(1),:), 50.016, 3, t );
    t_is( cq(L(2),1:2), [10 0.63], 2, t );
    t_is( cp(L(2),:), 55, 4, t );
    t_is( cp(G,1), bus(Gbus, LAM_P), 8, [t ' : gen prices'] );
    t_is( cp(L,1), bus(Lbus, LAM_P) + Qfudge(L,1), 8, [t ' : load prices'] );

    lao_X = p(G(1),2)/bus(Gbus(1), LAM_P);
    fro_X = p(G(6),3)/bus(Gbus(6), LAM_P);
    lab_X = p(L(2),2)/(bus(Lbus(2), LAM_P) + Qfudge(L(2),1));
    frb_X = p(L(3),3)/(bus(Lbus(3), LAM_P) + Qfudge(L(3),1));

    t_is( lao_X, 0.9997, 4, 'lao_X');
    t_is( fro_X, 1.1111, 4, 'fro_X');
    t_is( lab_X, 1, 4, 'lab_X');
    t_is( frb_X, 0.8960, 4, 'frb_X');

    t = 'one marginal bid @ $55, auction_type = 1';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1110, 100, [], [], mpopt);
    cp1 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5*lao_X, 8, [t ' : prices'] );

    t = 'one marginal bid @ $55, auction_type = 2';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1120, 100, [], [], mpopt);
    cp2 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G,:), cp5(G,:)*fro_X, 8, [t ' : gen prices'] );
    t_is( cp(L(1),:), cp5(L(1),:)*fro_X, 8, [t ' : load 1 price'] );
    t_is( cp(L(2),:), 55, 5, [t ' : load 2 price'] );
    t_is( cp(L(3),:), 60, 5, [t ' : load 3 price'] );

    t = 'one marginal bid @ $55, auction_type = 3';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1130, 100, [], [], mpopt);
    cp3 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5, 8, [t ' : prices'] );

    t = 'one marginal bid @ $55, auction_type = 4';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1140, 100, [], [], mpopt);
    cp4 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G(1),:), 50, 5, [t ' : gen 1 price'] );
    t_is( cp(G(2:6),:), cp5(G(2:6),:)*frb_X, 8, [t ' : gen 2-6 prices'] );
    t_is( cp(L,:), cp5(L,:)*frb_X, 8, [t ' : load prices'] );

    t = 'one marginal bid @ $55, auction_type = 6';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1160, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp1, 8, [t ' : prices'] );

    p2 = p;
    p2(G,:) = [ 0 0 100;
                0 0 100;
                0 0 100;
                0 0 100;
                0 0 100;
                0 0 100 ];
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1160, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G,:), cp5(G,:)*frb_X, 4, [t ' : gen prices'] );  %% gen 1, not clipped this time
    t_is( cp(L,:), cp4(L,:), 4, [t ' : load prices'] );

    t = 'one marginal bid @ $55, auction_type = 7';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1170, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5 * (lao_X+lab_X)/2, 8, [t ' : prices'] );
    t_is( cp, (cp1 + cp3) / 2, 8, [t ' : prices'] );

    t = 'one marginal bid @ $55, auction_type = 8';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1180, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G,:), cp1(G,:), 8, [t ' : gen prices'] );
    t_is( cp(L,:), cp3(L,:), 8, [t ' : load prices'] );

    t = 'one marginal bid @ $55, auction_type = 0';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1100, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, p, 8, [t ' : prices'] );

    %%-----  one bid block marginal @ $54.50 and one offer block marginal @ $50  -----
    p(L(2),2) = 54.5;
    t = 'marginal offer @ $50, bid @ $54.50, auction_type = 5';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1150, 100, [], [], mpopt);
    cq5 = cq;
    cp5 = cp;
    Qfudge =  zeros(size(p));
    Qfudge(L,:) = diag(gen(L,QG) ./ gen(L,PG) .* bus(Lbus, LAM_Q)) * ones(size(p(L,:)));

    t_is( cq(G(1),2:3), [23.74 0], 2, t );
    t_is( cp(G(1),:), 50, 4, t );
    t_is( cq(L(2),1:2), [10 0.39], 2, t );
    t_is( cp(L(2),:), 54.5, 4, t );
    t_is( cp(G,1), bus(Gbus, LAM_P), 8, [t ' : gen prices'] );
    t_is( cp(L,1), bus(Lbus, LAM_P) + Qfudge(L,1), 8, [t ' : load prices'] );

    lao_X = p(G(1),2)/bus(Gbus(1), LAM_P);
    fro_X = p(G(6),3)/bus(Gbus(6), LAM_P);
    lab_X = p(L(2),2)/(bus(Lbus(2), LAM_P) + Qfudge(L(2),1));
    frb_X = p(L(3),3)/(bus(Lbus(3), LAM_P) + Qfudge(L(3),1));

    t_is( lao_X, 1, 4, 'lao_X');
    t_is( fro_X, 1.1221, 4, 'fro_X');
    t_is( lab_X, 1, 4, 'lab_X');
    t_is( frb_X, 0.8976, 4, 'frb_X');

    t = 'marginal offer @ $50, bid @ $54.50, auction_type = 1';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1110, 100, [], [], mpopt);
    cp1 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5, 4, [t ' : prices'] );

    t = 'marginal offer @ $50, bid @ $54.50, auction_type = 2';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1120, 100, [], [], mpopt);
    cp2 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G,:), cp5(G,:)*fro_X, 5, [t ' : gen prices'] );
    t_is( cp(L(1),:), cp5(L(1),:)*fro_X, 5, [t ' : load 1 price'] );
    t_is( cp(L(2),:), 54.5, 5, [t ' : load 2 price'] );
    t_is( cp(L(3),:), 60, 5, [t ' : load 3 price'] );

    t = 'marginal offer @ $50, bid @ $54.50, auction_type = 3';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1130, 100, [], [], mpopt);
    cp3 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5, 8, [t ' : prices'] );

    t = 'marginal offer @ $50, bid @ $54.50, auction_type = 4';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1140, 100, [], [], mpopt);
    cp4 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G(1),:), 50, 5, [t ' : gen 1 price'] );
    t_is( cp(G(2:5),:), cp5(G(2:5),:)*frb_X, 8, [t ' : gen 2-5 prices'] );
    t_is( cp(G(6),:), 48, 5, [t ' : gen 6 price'] );
    t_is( cp(L,:), cp5(L,:)*frb_X, 8, [t ' : load prices'] );

    t = 'marginal offer @ $50, bid @ $54.50, auction_type = 6';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1160, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5, 4, [t ' : prices'] );

    t = 'marginal offer @ $50, bid @ $54.50, auction_type = 7';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1170, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5, 4, [t ' : prices'] );

    t = 'marginal offer @ $50, bid @ $54.50, auction_type = 8';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1180, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5, 4, [t ' : prices'] );

    t = 'marginal offer @ $50, bid @ $54.50, auction_type = 0';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1100, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, p, 8, [t ' : prices'] );

    %%-----  gen 1 at Pmin, load 3 block 2 marginal @ $60  -----
    t = 'gen 1 @ Pmin, marginal bid @ $60, auction_type = 5';
    p(L(2),2) = 50;     %% undo previous change
    p2 = p;
    p2(G(1),2:3) = [65 65];
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1150, 100, [], [], mpopt);
    Qfudge =  zeros(size(p));
    Qfudge(L,:) = diag(gen(L,QG) ./ gen(L,PG) .* bus(Lbus, LAM_Q)) * ones(size(p(L,:)));

    t_is( cp(G(1),:), 65, 4, [t ' : gen 1 price'] );
    t_is( cp(G(2),:), 54.2976, 4, [t ' : gen 2 price'] );
    cq5 = cq;
    cp5 = cp;
    cp_lam = cp5;
    cp_lam(1,:) = bus(Gbus(1), LAM_P);  %% unclipped

    lao_X = p2(G(6),2)/bus(Gbus(6), LAM_P);
    fro_X = p2(G(6),3)/bus(Gbus(6), LAM_P);
    lab_X = p2(L(3),2)/(bus(Lbus(3), LAM_P) + Qfudge(L(3),1));
    frb_X = p2(L(2),2)/(bus(Lbus(2), LAM_P) + Qfudge(L(2),1));

    t_is( lao_X, 0.8389, 4, 'lao_X');
    t_is( fro_X, 1.0487, 4, 'fro_X');
    t_is( lab_X, 1, 4, 'lab_X');
    t_is( frb_X, 0.8569, 4, 'frb_X');

    t = 'gen 1 @ Pmin, marginal bid @ $60, auction_type = 1';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1110, 100, [], [], mpopt);
    cp1 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G(1),:), 65, 8, [t ' : gen 1 price'] );
    t_is( cp(G(2:6),:), cp_lam(G(2:6),:)*lao_X, 8, [t ' : gen 2-6 prices'] );
    t_is( cp(L,:), cp_lam(L,:)*lao_X, 8, [t ' : load prices'] );

    t = 'gen 1 @ Pmin, marginal bid @ $60, auction_type = 2';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1120, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G(1),:), 65, 8, [t ' : gen 1 price'] );
    t_is( cp(G(2:6),:), cp_lam(G(2:6),:)*fro_X, 8, [t ' : gen 2-6 prices'] );
    t_is( cp(L(1:2),:), cp_lam(L(1:2),:)*fro_X, 8, [t ' : load 1-2 prices'] );
    t_is( cp(L(3),:), 60, 8, [t ' : load 3 price'] );

    t = 'gen 1 @ Pmin, marginal bid @ $60, auction_type = 3';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1130, 100, [], [], mpopt);
    cp3 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G(1),:), 65, 8, [t ' : gen 1 price'] );
    t_is( cp(G(2:6),:), cp_lam(G(2:6),:), 8, [t ' : gen 2-6 prices'] );
    t_is( cp(L,:), cp_lam(L,:), 8, [t ' : load prices'] );

    t = 'gen 1 @ Pmin, marginal bid @ $60, auction_type = 4';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1140, 100, [], [], mpopt);
    cp4 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G(1),:), 65, 5, [t ' : gen 1 price'] );
    t_is( cp(G(2:6),:), cp5(G(2:6),:)*frb_X, 8, [t ' : gen 2-6 prices'] );
    t_is( cp(L,:), cp5(L,:)*frb_X, 8, [t ' : load prices'] );

    t = 'gen 1 @ Pmin, marginal bid @ $60, auction_type = 6';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1160, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp4, 8, [t ' : prices'] );

    t = 'gen 1 @ Pmin, marginal bid @ $60, auction_type = 7';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1170, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G(1),:), 65, 4, [t ' : gen 1 price'] );
    t_is( cp(G(2:6),:), cp_lam(G(2:6),:) * (lao_X+lab_X)/2, 8, [t ' : gen 2-6 prices'] );
    t_is( cp(L,:), cp_lam(L,:) * (lao_X+lab_X)/2, 8, [t ' : load prices'] );
    t_is( cp, (cp1 + cp3) / 2, 8, [t ' : prices'] );

    t = 'gen 1 @ Pmin, marginal bid @ $60, auction_type = 8';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1180, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G,:), cp1(G,:), 8, [t ' : prices'] );
    t_is( cp(L,:), cp3(L,:), 8, [t ' : prices'] );

    t = 'gen 1 @ Pmin, marginal bid @ $60, auction_type = 0';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1100, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, p2, 8, [t ' : prices'] );

    %%-----  gen 1 at Pmin, gen 6 block 3 marginal @ $60  -----
    t = 'gen 1 @ Pmin, marginal offer @ $60, auction_type = 5';
    p2(L,:) = [ 100 100 100;
                100   0   0;
                100 100   0 ];
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1150, 100, [], [], mpopt);
    Qfudge =  zeros(size(p));
    Qfudge(L,:) = diag(gen(L,QG) ./ gen(L,PG) .* bus(Lbus, LAM_Q)) * ones(size(p(L,:)));

    t_is( cp(G(1),:), 65, 4, [t ' : gen 1 price'] );
    t_is( cp(G(2),:), 57.1616, 4, [t ' : gen 2 price'] );
    cq5 = cq;
    cp5 = cp;
    cp_lam = cp5;
    cp_lam(1,:) = bus(Gbus(1), LAM_P);  %% unclipped

    lao_X = p2(G(6),3)/bus(Gbus(6), LAM_P);
    fro_X = p2(G(1),3)/bus(Gbus(1), LAM_P);
    lab_X = p2(L(3),2)/(bus(Lbus(3), LAM_P) + Qfudge(L(3),1));
    frb_X = p2(L(2),2)/(bus(Lbus(2), LAM_P) + Qfudge(L(2),1));

    t_is( lao_X, 1, 4, 'lao_X');
    t_is( fro_X, 1.1425, 4, 'fro_X');
    t_is( lab_X, 1.5813, 4, 'lab_X');
    t_is( frb_X, 0, 4, 'frb_X');

    t = 'gen 1 @ Pmin, marginal offer @ $60, auction_type = 1';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1110, 100, [], [], mpopt);
    cp1 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp5, 8, [t ' : prices'] );

    t = 'gen 1 @ Pmin, marginal offer @ $60, auction_type = 2';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1120, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp_lam*fro_X, 8, [t ' : prices'] );

    t = 'gen 1 @ Pmin, marginal offer @ $60, auction_type = 3';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1130, 100, [], [], mpopt);
    cp3 = cp;
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp_lam*lab_X, 8, [t ' : prices'] );

    t = 'gen 1 @ Pmin, marginal offer @ $60, auction_type = 4';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1140, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G,1), [65;40;42;44;46;60], 4, [t ' : gen prices'] );
    t_is( cp(L,:), cp_lam(L,:)*frb_X, 8, [t ' : prices'] );

    t = 'gen 1 @ Pmin, marginal offer @ $60, auction_type = 6';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1160, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp_lam*fro_X, 8, [t ' : prices'] );

    t = 'gen 1 @ Pmin, marginal offer @ $60, auction_type = 7';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1170, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, cp_lam * (lao_X+lab_X)/2, 8, [t ' : prices'] );
    t_is( cp, (cp_lam + cp3) / 2, 8, [t ' : prices'] );

    t = 'gen 1 @ Pmin, marginal offer @ $60, auction_type = 8';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1180, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp(G,:), cp5(G,:), 8, [t ' : prices'] );
    t_is( cp(L,:), cp3(L,:), 8, [t ' : prices'] );

    t = 'gen 1 @ Pmin, marginal offer @ $60, auction_type = 0';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p2, 1100, 100, [], [], mpopt);
    t_is( cq, cq5, 8, [t ' : quantities'] );
    t_is( cp, p2, 8, [t ' : prices'] );

    %%-----  gen 2 decommitted, one offer block marginal @ $60  -----
    p(G(2),:) = p(G(2),:) + 100;

    t = 'price of decommited gen, auction_type = 5';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1150, 200, [], [], mpopt);
    cp5 = cp;
    Qfudge =  zeros(size(p));
    Qfudge(L,:) = diag(gen(L,QG) ./ gen(L,PG) .* bus(Lbus, LAM_Q)) * ones(size(p(L,:)));
    t_is(sum(cq(2,:)), 0, 8, t);
    t_is(cp(2,1), 59.194, 3, t);

%     Xo = p(1:6, :) ./ (diag(bus(Gbus, LAM_P)) * ones(size(p(G,:))));
%     ao = (cq(1:6, :) ~= 0);
%     ro = (cq(1:6, :) == 0);
%     Xb = p(7:9, :) ./ (diag(bus(Lbus, LAM_P) + gen(L,QG) ./ gen(L,PG) .* bus(Lbus, LAM_Q)) * ones(size(p(L,:))));
%     ab = (cq(7:9, :) ~= 0);
%     rb = (cq(7:9, :) == 0);
%     aXo = ao .* Xo
%     rXo = ro .* Xo
%     aXb = ab .* Xb
%     rXb = rb .* Xb

    lao_X = p(G(6),3)/bus(Gbus(6), LAM_P);
    fro_X = p(G(1),3)/bus(Gbus(1), LAM_P);
    lab_X = p(L(1),2)/(bus(Lbus(1), LAM_P) + Qfudge(L(1),1));
    frb_X = p(L(1),3)/(bus(Lbus(1), LAM_P) + Qfudge(L(1),1));

    t_is( lao_X, 1, 4, 'lao_X');
    t_is( fro_X, 1.0212, 4, 'fro_X');
    t_is( lab_X, 1.1649, 4, 'lab_X');
    t_is( frb_X, 0.9985, 4, 'frb_X');

    t = 'price of decommited gen, auction_type = 1';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1110, 200, [], [], mpopt);
    t_is(cp(2,1), 59.194, 3, t);

    t = 'price of decommited gen, auction_type = 2';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1120, 200, [], [], mpopt);
    t_is(cp(2,1), cp5(2,1)*fro_X, 3, t);

    t = 'price of decommited gen, auction_type = 3';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1130, 200, [], [], mpopt);
    t_is(cp(2,1), cp5(2,1)*lab_X, 3, t);

    t = 'price of decommited gen, auction_type = 4';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1140, 200, [], [], mpopt);
    t_is(cp(2,1), cp5(2,1)*frb_X, 3, t);

    t = 'price of decommited gen, auction_type = 6';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1160, 200, [], [], mpopt);
    t_is(cp(2,1), cp5(2,1)*fro_X, 3, t);

    t = 'price of decommited gen, auction_type = 7';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1170, 200, [], [], mpopt);
    t_is(cp(2,1), cp5(2,1)*(lao_X+lab_X)/2, 3, t);

    t = 'price of decommited gen, auction_type = 0';
    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1100, 200, [], [], mpopt);
    t_is(cp(2,1), 120, 3, t);

    t = 'single block, marginal offer @ $50, auction_type = 5';
    q = [
        60;
        36;
        36;
        36;
        36;
        36;
        30;
        10;
        20;
    ];

    p = [
        50;
        40;
        42;
        44;
        46;
        48;
        100;
        100;
        100;
    ];

    [MVAbase, cq, cp, bus, gen, gencost, branch, f, dispatch, success, et] = ...
        runmkt('t_auction_case', q, p, 1150, 100, [], [], mpopt);
    t_is( cq(G(1)), 35.32, 2, t );
    t_is( cq(G(2:6)), q(G(2:6)), 8, [t ' : gen qtys'] );
    t_is( cp(G(1)), 50, 4, t );
    t_is( cq(L), q(L), 8, [t ' : load qtys'] );
    t_is( cp(L(2),:), 54.03, 2, t );
    t_is( cp(G), bus(Gbus, LAM_P), 8, [t ' : gen prices'] );
    Qfudge =  zeros(size(p));
    Qfudge(L,:) = diag(gen(L,QG) ./ gen(L,PG) .* bus(Lbus, LAM_Q)) * ones(size(p(L,:)));
    t_is( cp(L), bus(Lbus, LAM_P) + Qfudge(L,1), 8, [t ' : load prices'] );
end

t_end;