text:
function step_distance(nlp, bounds)
    % constraints for step length and step width

    domain = nlp.Plant;
    x = domain.States.x;

    [right_foot_frame] = sys.frames.RightFoot(domain);
    [left_foot_frame] = sys.frames.LeftFoot(domain);
    right_foot = getCartesianPosition(domain, right_foot_frame);
    left_foot = getCartesianPosition(domain, left_foot_frame);
    constraint = tomatrix(left_foot(1:2) - right_foot(1:2));
    constraint_func = SymFunction(['step_distance_',domain.Name], constraint, {x});
    lb = [bounds.constrBounds.stepWidth.lb
        bounds.constrBounds.stepLength.lb];
    ub = [bounds.constrBounds.stepWidth.ub
        bounds.constrBounds.stepLength.ub];
    addNodeConstraint(nlp, constraint_func, {'x'}, 'first', lb, ub, 'NonLinear');

    %     wt = 0.3662;
    %     lb = [-wt-0.015
    %         0-0.03];
    %     ub = [-wt+0.015
    %         0+0.03];
    %     addNodeConstraint(nlp, constraint_func, {'x'}, floor(nlp.NumNode/2)+1, lb, ub, 'NonLinear');
end


text:
% [PYR, INDICES, STEERMTX, HARMONICS] = buildFullSFpyr2(IM, HEIGHT, ORDER, TWIDTH)
%
% Construct a steerable pyramid on matrix IM, in the Fourier domain.
% Unlike the standard transform, subdivides the highpass band into
% orientations.

function [pyr,pind,steermtx,harmonics] = buildFullSFpyr2(im, ht, order, twidth)

%-----------------------------------------------------------------
%% DEFAULTS:

max_ht = floor(log2(min(size(im)))+1);

if (exist('ht') ~= 1)
  ht = max_ht;
else
  if (ht > max_ht)
    error(sprintf('Cannot build pyramid higher than %d levels.',max_ht));
  end
end

if (exist('order') ~= 1)
  order = 3;
elseif ((order > 15)  | (order < 0))
  fprintf(1,'Warning: ORDER must be an integer in the range [0,15]. Truncating.\n');
  order = min(max(order,0),15);
else
  order = round(order);
end
nbands = order+1;

if (exist('twidth') ~= 1)
  twidth = 1;
elseif (twidth <= 0)
  fprintf(1,'Warning: TWIDTH must be positive.  Setting to 1.\n');
  twidth = 1;
end

%-----------------------------------------------------------------
%% Steering stuff:

if (mod((nbands),2) == 0)
  harmonics = [0:(nbands/2)-1]'*2 + 1;
else
  harmonics = [0:(nbands-1)/2]'*2;
end

steermtx = steer2HarmMtx(harmonics, pi*[0:nbands-1]/nbands, 'even');

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

dims = size(im);
ctr = ceil((dims+0.5)/2);

[xramp,yramp] = meshgrid( ([1:dims(2)]-ctr(2))./(dims(2)/2), ...
    ([1:dims(1)]-ctr(1))./(dims(1)/2) );
angle = atan2(yramp,xramp);
log_rad = sqrt(xramp.^2 + yramp.^2);
log_rad(ctr(1),ctr(2)) =  log_rad(ctr(1),ctr(2)-1);
log_rad  = log2(log_rad);

%% Radial transition function (a raised cosine in log-frequency):
[Xrcos,Yrcos] = rcosFn(twidth,(-twidth/2),[0 1]);
Yrcos = sqrt(Yrcos);

YIrcos = sqrt(1.0 - Yrcos.^2);
lo0mask = pointOp(log_rad, YIrcos, Xrcos(1), Xrcos(2)-Xrcos(1), 0);
imdft = fftshift(fft2(im));
lo0dft =  imdft .* lo0mask;

[pyr,pind] = buildSFpyrLevs(lo0dft, log_rad, Xrcos, Yrcos, angle, ht, nbands);

%% Split the highpass band into orientations

hi0mask = pointOp(log_rad, Yrcos, Xrcos(1), Xrcos(2)-Xrcos(1), 0);

lutsize = 1024;
Xcosn = pi*[-(2*lutsize+1):(lutsize+1)]/lutsize;  % [-2*pi:pi]
order = nbands-1;
const = (2^(2*order))*(factorial(order)^2)/(nbands*factorial(2*order));
Ycosn = sqrt(const) * (cos(Xcosn)).^order;

bands = zeros(prod(size(imdft)), nbands);
bind = zeros(nbands,2);

for b = 1:nbands
    anglemask = pointOp(angle, Ycosn, Xcosn(1)+pi*(b-1)/nbands, Xcosn(2)-Xcosn(1));
    Mask = ((-sqrt(-1))^(nbands-1))*anglemask.*hi0mask;
    % make real the contents in the HF cross (to avoid information loss in these freqs.)
    % It distributes evenly these contents among the nbands orientations
    Mask(1,:) = ones(1,size(im,2))/sqrt(nbands);
    Mask(2:size(im,1),1) = ones(size(im,1)-1,1)/sqrt(nbands);

    banddft =  imdft .* Mask;
    band = real(ifft2(fftshift(banddft)));
    bands(:,b) = real(band(:));
    bind(b,:)  = size(band);
end

pyr = [bands(:); pyr];
pind = [bind; pind];

pind = [ [0 0]; pind];  %% Dummy highpass

author: ricedsp


text:
% GET THE STANDARD 2D TRAJECTORY
function [currentFigure,figureHandle] = GetFigure_plan(SIM,objectIndex,DATA,currentFigure)

% CONFIGURE THE PLOT ATTRIBUTES
figureHandle = figure('Name','OpenMAS top-down image');
set(figureHandle,'Position',DATA.figureProperties.windowSettings);         % [x y width height]
set(figureHandle,'Color',DATA.figureProperties.figureColor);               % Background colour
ax = axes(figureHandle);
% setappdata(figureHandle, 'SubplotDefaultAxesLocation', [0.08, 0.08, 0.90, 0.88]); % MAXIMISE GRAPH SIZE IN WINDOW

legendCounter = 1; legendEntries = cell(DATA.totalObjects,1);

hold on; grid on; box on; grid minor;
for ID1 = 1:DATA.totalObjects
    % GET OBJECT OVERVIEW DATA
    legendString = sprintf('[ID-%s] %s',num2str(SIM.OBJECTS(ID1).objectID),SIM.OBJECTS(ID1).name);
    legendEntries(legendCounter) = {legendString};
    % THE OBJECT HANDLE
    objectHandle = objectIndex{SIM.OBJECTS(ID1).objectID == SIM.globalIDvector};
    % EXTRACT FINAL POSITION DATA FROM THE TRAJECTORY MATRIX
    [finalStates] = OMAS_getTrajectoryData_mex(...
        DATA.globalTrajectories,...
        SIM.globalIDvector,...
        SIM.OBJECTS(ID1).objectID,...
        SIM.TIME.endStep);
    finalPosition = finalStates(1:3,:);
    % LOCAL FIXED TO GLOBAL ROTATED
    R_final = OMAS_geometry.quaternionToRotationMatrix(finalStates(7:10));
    % DISPLAY THE OBJECT
    if numel(objectHandle.GEOMETRY.vertices) > 0
        patch(ax,'Vertices',objectHandle.GEOMETRY.vertices*R_final + finalPosition',...
            'Faces',objectHandle.GEOMETRY.faces,...
            'FaceColor',SIM.OBJECTS(ID1).colour,...
            'EdgeColor',DATA.figureProperties.edgeColor,...
            'EdgeAlpha',DATA.figureProperties.edgeAlpha,...
            'FaceLighting',DATA.figureProperties.faceLighting,...
            'FaceAlpha',DATA.figureProperties.faceAlpha,...
            'LineWidth',DATA.figureProperties.patchLineWidth);             % Patch properties
    else
        % PLOT THE TERMINAL POSITIONS
        plot3(ax,finalPosition(1),finalPosition(2),finalPosition(3),...
              'Marker',SIM.OBJECTS(ID1).symbol,...
              'MarkerSize',DATA.figureProperties.markerSize,...
              'MarkerFaceColor',SIM.OBJECTS(ID1).colour,...
              'MarkerEdgeColor',DATA.figureProperties.markerEdgeColor,...
              'LineWidth',DATA.figureProperties.lineWidth,...
              'LineStyle',DATA.figureProperties.lineStyle,...
              'Color',SIM.OBJECTS(ID1).colour);
    end
    legendCounter = legendCounter + 1;
end
legend('off')
hold on;
for ID1 = 1:DATA.totalObjects
    % EXTRACT TIME-STATE DATA FROM THE TRAJECTORY MATRIX
    idleFlag = NaN('double');
    [objectStates] = OMAS_getTrajectoryData_mex(DATA.globalTrajectories,SIM.globalIDvector,SIM.OBJECTS(ID1).objectID,idleFlag);
    positions = objectStates(1:3,:);
    plot3(positions(1,:),positions(2,:),positions(3,:),...
          'LineStyle',DATA.figureProperties.lineStyle,...
          'LineWidth',DATA.figureProperties.lineWidth,...
          'Color',SIM.OBJECTS(ID1).colour);

    % ADD THE INTIAL POSITION REFERENCE POSITION
    if SIM.OBJECTS(ID1).type == OMAS_objectType.agent
        initialPlotPosition = positions(1:2,1) - SIM.OBJECTS(ID1).radius; % The first position of the object
        rectangle('Position',[initialPlotPosition(1) initialPlotPosition(2) (2*SIM.OBJECTS(ID1).radius) (2*SIM.OBJECTS(ID1).radius)],...
                  'Curvature',[1 1],...
                  'FaceColor',SIM.OBJECTS(ID1).colour,...
                  'EdgeColor',DATA.figureProperties.edgeColor,...
                  'LineWidth',DATA.figureProperties.patchLineWidth)
    end
end

% Title
title(ax,...
    sprintf('Object trajectories over a period of %ss',num2str(SIM.TIME.endTime)),...
    'Interpreter',DATA.figureProperties.interpreter,...
    'fontname',DATA.figureProperties.fontName,...
    'Fontweight',DATA.figureProperties.fontWeight,...
    'FontSize',DATA.figureProperties.titleFontSize,...
    'FontSmoothing','on');
% X-Label
xlabel(ax,'x(m)',...
    'Interpreter',DATA.figureProperties.interpreter,...
    'fontname',DATA.figureProperties.fontName,...
    'Fontweight',DATA.figureProperties.fontWeight,...
    'FontSize',DATA.figureProperties.axisFontSize,...
    'FontSmoothing','on');
% Y-Label
ylabel(ax,'y(m)',...
    'Interpreter',DATA.figureProperties.interpreter,...
    'fontname',DATA.figureProperties.fontName,...
    'Fontweight',DATA.figureProperties.fontWeight,...
    'FontSize',DATA.figureProperties.axisFontSize,...
    'FontSmoothing','on');
% Axes
set(ax,...
    'TickLabelInterpreter',DATA.figureProperties.interpreter,...
    'fontname',DATA.figureProperties.fontName,...
    'Fontweight',DATA.figureProperties.fontWeight,...
    'FontSize',DATA.figureProperties.axisFontSize,...
    'FontSmoothing','on',...
    'Color',DATA.figureProperties.axesColor,...
    'GridLineStyle','--',...
    'GridAlpha',0.25,...
    'GridColor','k');
% Legend
legend(legendEntries,...
    'location','northeast',...
    'interpreter',DATA.figureProperties.interpreter,...
    'fontname',DATA.figureProperties.fontName);
axis square equal;

% axis square;
view([0 90]);
% set(ax,'outerposition',[0.05 0.15 1 0.68]);                               % Set the axes offset position in the figure window
hold off;

figurePath = strcat(SIM.outputPath,'plan');
% SAVE THE OUTPUT FIGURE
savefig(figureHandle,figurePath);

% Publish as .pdf if requested
if DATA.figureProperties.publish
	GetFigurePDF(figureHandle,figurePath);
end

% FIGURE COMPLETE
DATA.figureProperties.alignment = DATA.figureProperties.alignment + DATA.figureProperties.spacing;
currentFigure = currentFigure + 1;
end

repo: OpenMAS


text:
%% Distributed Marginal Value-at-Risk Simulation
% This example uses the Parallel Computing Toolbox(TM) to perform a Monte
% Carlo simulation of a number of stocks in a portfolio. At a given
% confidence level, we predict the value at risk (VaR) of the portfolio as
% well as the marginal value at risk (mVaR) of each of the stocks in the
% portfolio.  We also provide confidence intervals for our estimates.
%
%   Copyright 2007-2012 The MathWorks, Inc.
%   Adapted and simplified in 2013 by Michael Weidman, Quantitative Support
%   Services, Ltd.

%% 1. Open a pool of MATLAB workers and load input data

if matlabpool('size') == 0
    matlabpool open
end

load pctdemo_data_mvar

%% 2. Define Parameters

nSims           = 1e5;
relativeWeights = [1 1 2 2 1 3 2]/12;
time            = 50;%(10:1:50)';
confLevel       = 95;

nTimes          = length(time);
nAssets         = size(stock,2);

%% 3. Estimate asset moments

returns         = (stock(2:end,:) - stock(1:end-1,:)) ./ stock(1:end-1,:);
expReturns      = mean(returns);
expCovariances  =  cov(returns);

%% 4. Simulate asset prices

tic

simPrices   = zeros(nTimes, nAssets, nSims);
portPrices  = zeros(nTimes,       1, nSims);
price0      = stock(end,:);
portPrice0  = price0 * relativeWeights';

parfor iSim = 1:nSims
    simPrices( :,:,iSim) = simulatePrices(expReturns, expCovariances, ...
        price0, time);
    portPrices(:,:,iSim) = sum(bsxfun(@times, simPrices(:,:,iSim), ...
        relativeWeights), 2);
end

toc

%% 5. Calculate mVaR and VaR

simRelativeLosses = bsxfun(@rdivide, ...
    bsxfun(@minus, price0, simPrices), price0);
simRelativePortLosses = bsxfun(@minus, portPrice0, portPrices) ...
    ./ portPrice0;

mVaR = prctile(simRelativeLosses    , confLevel, 3);
VaR  = prctile(simRelativePortLosses, confLevel, 3);

%% 6. Plot the Results

plotMVAR(VaR, mVaR, time, names);

author: Sable

repo: mcbench-benchmarks


text:
function qq = ppdiff(pp,j)
%PPDIFF Differentiate piecewise polynomial.
%   QQ = PPDIFF(PP,J) returns the J:th derivative of a piecewise
%   polynomial PP. PP must be on the form evaluated by PPVAL. QQ is a
%   piecewise polynomial on the same form. Default value for J is 1.
%
%   Example:
%       x = linspace(-pi,pi,9);
%       y = sin(x);
%       pp = spline(x,y);
%       qq = ppdiff(pp);
%       xx = linspace(-pi,pi,201);
%       plot(xx,cos(xx),'b',xx,ppval(qq,xx),'r')
%
%   See also PPVAL, SPLINE, SPLINEFIT, PPINT

%   Author: Jonas Lundgren <splinefit@gmail.com> 2009

if nargin < 1, help ppdiff, return, end
if nargin < 2, j = 1; end

% Check diff order
if ~isreal(j) || mod(j,1) || j < 0
    msgid = 'PPDIFF:DiffOrder';
    message = 'Order of derivative must be a non-negative integer!';
    error(msgid,message)
end

% Get coefficients
coefs = pp.coefs;
[m n] = size(coefs);

if j == 0
    % Do nothing
elseif j < n
    % Derivative of order J
    D = [n-j:-1:1; ones(j-1,n-j)];
    D = cumsum(D,1);
    D = prod(D,1);
    coefs = coefs(:,1:n-j);
    for k = 1:n-j
        coefs(:,k) = D(k)*coefs(:,k);
    end
else
    % Derivative kills PP
    coefs = zeros(m,1);
end

% Set output
qq = pp;
qq.coefs = coefs;
qq.order = size(coefs,2);


text:
function [m,n] = size(a,dim)
%SIZE         Implements  size(a)  for intervals
%
%   [m,n] = size(a,dim)
%
% functionality as Matlab function size, second argument dim optional
%

% written  10/16/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if a.complex
    if nargout==2
      [m n] = size(a.mid);
    else
      if nargin==1
        m = size(a.mid);
      else
        m = size(a.mid,dim);
      end
    end
  else
    if nargout==2
      [m n] = size(a.inf);
    else
      if nargin==1
        m = size(a.inf);
      else
        m = size(a.inf,dim);
      end
    end
  end


text:
% Copyright (C) 1994-2015 John W. Eaton
%
% This file is part of Octave.
%
% Octave 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.
%
% Octave 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 Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

% -*- texinfo -*-
% @deftypefn  {Function File} {[@var{h}, @var{w}] =} freqz (@var{b}, @var{a}, @var{n}, "whole")
% @deftypefnx {Function File} {[@var{h}, @var{w}] =} freqz (@var{b})
% @deftypefnx {Function File} {[@var{h}, @var{w}] =} freqz (@var{b}, @var{a})
% @deftypefnx {Function File} {[@var{h}, @var{w}] =} freqz (@var{b}, @var{a}, @var{n})
% @deftypefnx {Function File} {@var{h} =} freqz (@var{b}, @var{a}, @var{w})
% @deftypefnx {Function File} {[@var{h}, @var{w}] =} freqz (@dots{}, @var{Fs})
% @deftypefnx {Function File} {} freqz (@dots{})
%
% Return the complex frequency response @var{h} of the rational IIR filter
% whose numerator and denominator coefficients are @var{b} and @var{a},
% respectively.
%
% The response is evaluated at @var{n} angular frequencies between 0 and
% @ifnottex
% 2*pi.
% @end ifnottex
% @tex
% $2\pi$.
% @end tex
%
% @noindent
% The output value @var{w} is a vector of the frequencies.
%
% If @var{a} is omitted, the denominator is assumed to be 1 (this
% corresponds to a simple FIR filter).
%
% If @var{n} is omitted, a value of 512 is assumed.  For fastest computation,
% @var{n} should factor into a small number of small primes.
%
% If the fourth argument, @qcode{'whole'}, is omitted the response is
% evaluated at frequencies between 0 and
% @ifnottex
% pi.
% @end ifnottex
% @tex
% $\pi$.
% @end tex
%
% @code{freqz (@var{b}, @var{a}, @var{w})}
%
% Evaluate the response at the specific frequencies in the vector @var{w}.
% The values for @var{w} are measured in radians.
%
% @code{[@dots{}] = freqz (@dots{}, @var{Fs})}
%
% Return frequencies in Hz instead of radians assuming a sampling rate
% @var{Fs}.  If you are evaluating the response at specific frequencies
% @var{w}, those frequencies should be requested in Hz rather than radians.
%
% @code{freqz (@dots{})}
%
% Plot the magnitude and phase response of @var{h} rather than returning them.
%
% @seealso{freqz_plot}
% @end deftypefn

% Author: jwe ???

function [h_r, f_r] = freqz (b, a, n, region, Fs)

  if (nargin < 1 || nargin > 5)
    print_usage ();
  elseif (nargin == 1)
    % Response of an FIR filter.
    a = [];
    n = [];
    region = [];
    Fs = [];
  elseif (nargin == 2)
    % Response of an IIR filter
    n = [];
    region = [];
    Fs = [];
  elseif (nargin == 3)
    region = [];
    Fs = [];
  elseif (nargin == 4)
    Fs = [];
    if (~ischar (region) && ~isempty (region))
      Fs = region;
      region = [];
    end
  end

  if (isempty (b))
    b = 1;
  elseif (~isvector (b))
    error ('freqz: B must be a vector');
  end
  if (isempty (a))
    a = 1;
  elseif (~isvector (a))
    error ('freqz: A must be a vector');
  end
  if (isempty (n))
    n = 512;
  elseif (isscalar (n) && n < 1)
    error ('freqz: N must be a positive integer');
  end
  if (isempty (region))
    if (isreal (b) && isreal (a))
      region = 'half';
    else
      region = 'whole';
    end
  end
  if (isempty (Fs))
    freq_norm = true;
    if (nargout == 0)
      Fs = 2;
    else
      Fs = 2*pi;
    end
  else
    freq_norm = false;
  end
  plot_output = (nargout == 0);
  whole_region = strcmp (region, 'whole');

  a = a(:);
  b = b(:);

  if (~isscalar (n))
    % Explicit frequency vector given
    w = n;
    f = n;
    if (nargin == 4)
      % Sampling rate Fs was specified
      w = 2*pi*f/Fs;
    end
    k = max (length (b), length (a));
    hb = polyval (oc_postpad (b, k), exp (j*w));
    ha = polyval (oc_postpad (a, k), exp (j*w));
  else
    % polyval(fliplr(P),exp(jw)) is O(p n) and fft(x) is O(n log(n)),
    % where p is the order of the polynomial P.  For small p it
    % would be faster to use polyval but in practice the overhead for
    % polyval is much higher and the little bit of time saved isn't
    % worth the extra code.
    k = max (length (b), length (a));
    if (k > n/2 && nargout == 0)
      % Ensure a causal phase response.
      n = n * 2 .^ ceil (log2 (2*k/n));
    end

    if (whole_region)
      N = n;
      if (plot_output)
        f = Fs * (0:n).' / N;    % do 1 more for the plot
      else
        f = Fs * (0:n-1).' / N;
      end
    else
      N = 2*n;
      if (plot_output)
        n = n+1;
      end
      f = Fs * (0:n-1).' / N;
    end

    pad_sz = N*ceil (k/N);
    b = oc_postpad (b, pad_sz);
    a = oc_postpad (a, pad_sz);

    hb = zeros (n, 1);
    ha = zeros (n, 1);

    for i = 1:N:pad_sz
      fftb = fft(oc_postpad (b(i:i+N-1), N));
      ffta = fft(oc_postpad (a(i:i+N-1), N));
      hb = hb + fftb(1:n);
      ha = ha + ffta(1:n);
    end

  end

  h = hb ./ ha;

  if (plot_output)
    % Plot and don't return values.
    if (whole_region)
      h(end+1) = h(1); % Solution is periodic.  Copy first value to end.
    end
    freqz_plot (f, h, freq_norm);
  else
    % Return values and don't plot.
    h_r = h;
    f_r = f;
  end

end

author: brainstorm-tools


text:
function points = intersectLineSphere(line, sphere, varargin)
    %INTERSECTLINESPHERE Return intersection points between a line and a sphere
    %
    %   PTS = intersectLineSphere(LINE, SPHERE);
    %   Returns the two points which are the intersection of the given line and
    %   sphere.
    %   LINE   : [x0 y0 z0  dx dy dz]
    %   SPHERE : [xc yc zc  R]
    %   PTS     : [x1 y1 z1 ; x2 y2 z2]
    %   If there is no intersection between the line and the sphere, return a
    %   2-by-3 array containing only NaN.
    %
    %   Example
    %     % draw the intersection between a sphere and a collection of parallel
    %     % lines
    %     sphere = [50.12 50.23 50.34 40];
    %     [x, y] = meshgrid(10:10:90, 10:10:90);
    %     n = numel(x);
    %     lines = [x(:) y(:) zeros(n,1) zeros(n,2) ones(n,1)];
    %     figure; hold on; axis equal;
    %     axis([0 100 0 100 0 100]); view(3);
    %     drawSphere(sphere);
    %     drawLine3d(lines);
    %     pts = intersectLineSphere(lines, sphere);
    %     drawPoint3d(pts, 'ro');
    %
    %   See also
    %   spheres, circles3d, intersectPlaneSphere
    %
    %   ---------
    %   author : David Legland
    %   INRA - TPV URPOI - BIA IMASTE
    %   created the 18/02/2005.
    %
    %   HISTORY
    %   2011-06-21 bug for tangent lines, add tolerance
    % Process input arguments
    % check if user-defined tolerance is given
    tol = 1e-14;
    if ~isempty(varargin)
        tol = varargin{1};
    end
    % difference between centers
    dc = bsxfun(@minus, line(:, 1:3), sphere(:, 1:3));
    % equation coefficients
    a = sum(line(:, 4:6) .* line(:, 4:6), 2);
    b = 2 * sum(bsxfun(@times, dc, line(:, 4:6)), 2);
    c = sum(dc.*dc, 2) - sphere(:,4).*sphere(:,4);
    % solve equation
    delta = b.*b - 4*a.*c;
    % initialize empty results
    points = NaN * ones(2 * size(delta, 1), 3);
    % process couples with two intersection points
    % proces couples with two intersection points
    inds = find(delta > tol);
    if ~isempty(inds)
        % delta positive: find two roots of second order equation
        u1 = (-b(inds) -sqrt(delta(inds))) / 2 ./ a(inds);
        u2 = (-b(inds) +sqrt(delta(inds))) / 2 ./ a(inds);

        % convert into 3D coordinate
        points(inds, :) = line(inds, 1:3) + bsxfun(@times, u1, line(inds, 4:6));
        points(inds+length(delta),:) = line(inds, 1:3) + bsxfun(@times, u2, line(inds, 4:6));
    end
    % process couples with one intersection point
    % proces couples with two intersection points
    inds = find(abs(delta) < tol);
    if ~isempty(inds)
        % delta around zero: find unique root, and convert to 3D coord.
        u = -b(inds) / 2 ./ a(inds);

        % convert into 3D coordinate
        pts = line(inds, 1:3) + bsxfun(@times, u, line(inds, 4:6));
        points(inds, :) = pts;
        points(inds+length(delta),:) = pts;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout = drawEllipseCylinder(cyl, varargin)
%DRAWELLIPSECYLINDER Draw a cylinder with ellipse cross-section.
%
%   drawEllipseCylinder(CYL)
%   draws the cylinder CYL on the current axis.
%   CYL is a cylinder defined by [x1 y1 z1 x2 y2 z2 r1 r2 roll], with:
%   * [x1 y2 z1] are coordinates of starting point,
%   * [x2 y2 z2] are coordinates of ending point,
%   * R1 and R2 are the lengths of the ellipse semi axes, and
%   * ROLL is the rotation of the cylinder around its main axis (in
%      degrees)
%
%   drawEllipseCylinder(CYL, N)
%   uses N points for discretisation of angle. Default value is 32.
%
%   drawEllipseCylinder(..., OPT)
%   with OPT = 'open' (default) or 'closed', specify if bases of the
%   cylinder should be drawn.
%
%   drawEllipseCylinder(..., 'FaceColor', COLOR)
%   Specifies the color of the cylinder. Any couple of parameters name and
%   value can be given as argument, and will be transfered to the 'surf'
%   matlab function
%
%   H = drawEllipseCylinder(...)
%   returns a handle to the patch representing the cylinder.
%
%
%   Example:
%   figure;drawEllipseCylinder([0 0 0 10 20 30 5]);
%
%   figure;drawEllipseCylinder([0 0 0 10 20 30 5], 'open');
%
%   figure;drawEllipseCylinder([0 0 0 10 20 30 5], 'FaceColor', 'r');
%
%   figure;
%   h = drawEllipseCylinder([0 0 0 10 20 30 5]);
%   set(h, 'facecolor', 'b');
%
%   % Draw three mutually intersecting elliptic cylinders
%     p1 = [30 0 0];
%     p2 = [0 30 0];
%     p3 = [0 0 30];
%     radii = [20 10];
%     figure;
%     drawEllipseCylinder([-p1 p1 radii 0], 'FaceColor', 'r');
%     hold on
%     drawEllipseCylinder([-p2 p2 radii 90], 'FaceColor', 'g');
%     drawEllipseCylinder([-p3 p3 radii 90], 'FaceColor', 'b');
%     axis equal
%     set(gcf, 'renderer', 'opengl')
%     view([60 30]); light;
%
%   See Also:
%   drawCylinder, drawSphere, cylinderMesh, drawLine3d, surf

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 27/02/2014

%   HISTORY

%% Input argument processing

if iscell(cyl)
    res = zeros(length(cyl), 1);
    for i = 1:length(cyl)
        res(i) = drawEllipseCylinder(cyl{i}, varargin{:});
    end

    if nargout > 0
        varargout{1} = res;
    end
    return;
end

% default values
N = 32;
closed = true;

% check number of discretization steps
if ~isempty(varargin)
    var = varargin{1};
    if isnumeric(var)
        N = var;
        varargin = varargin(2:end);
    end
end

% check if cylinder must be closed or open
if ~isempty(varargin)
    var = varargin{1};
    if ischar(var)
        if strncmpi(var, 'open', 4)
            closed = false;
            varargin = varargin(2:end);
        elseif strncmpi(var, 'closed', 5)
            closed = true;
            varargin = varargin(2:end);
        end
    end
end

%% Computation of mesh coordinates

% extreme points of cylinder
p1 = cyl(1:3);
p2 = cyl(4:6);

% radius of cylinder
r1 = cyl(7);
r2 = cyl(8);
roll = 0;
if size(cyl, 2) > 8
    roll = cyl(9);
end

% compute orientation angle of cylinder (in degrees)
[theta, phi, rho] = cart2sph2d(p2 - p1);
dphi = linspace(0, 2*pi, N+1);

% generate a cylinder oriented upwards
x = repmat(cos(dphi) * r1, [2 1]);
y = repmat(sin(dphi) * r2, [2 1]);
z = repmat([0 ; rho], [1 length(dphi)]);

% transform points
trans   = localToGlobal3d(p1, theta, phi, roll);
pts     = transformPoint3d([x(:) y(:) z(:)], trans);

% reshape transformed points
x2 = reshape(pts(:,1), size(x));
y2 = reshape(pts(:,2), size(x));
z2 = reshape(pts(:,3), size(x));

%% Display cylinder mesh

% add default drawing options
varargin = [{'FaceColor', 'g', 'edgeColor', 'none'} varargin];

% plot the cylinder as a surface
hSurf = surf(x2, y2, z2, varargin{:});

% eventually plot the ends of the cylinder
if closed
    ind = find(strcmpi(varargin, 'facecolor'), 1, 'last');
    if isempty(ind)
        color = 'k';
    else
        color = varargin{ind+1};
    end

    patch(x2(1,:)', y2(1,:)', z2(1,:)', color, 'edgeColor', 'none');
    patch(x2(2,:)', y2(2,:)', z2(2,:)', color, 'edgeColor', 'none');
end

% format ouptut
if nargout == 1
    varargout{1} = hSurf;
end


text:
function meanF1 = MAMUDA_enterFunc(inputPath,ite1,ite2,itermediateD,reducedD,sharedD)
    load(inputPath);
    [result,F1,meanF1,R,Rt,Q]=MAMUDA(testfea,testlabel,trainFea,trainLabel,unlabelFea,unlabelLabel,ite1,ite2,0.0,itermediateD,reducedD,sharedD,2);
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function varargout = subsref(f, index)
%SUBSREF   DELTAFUN subsref.
% ( )
%   F(x) tries to evaluate the distribution F at the points given in x.
%   Mathematically, this doesn't make sense when x coincides with a point having
%   a non-trivial delta function.
%
%   F(f) [TODO]: tries to compose the distribution F with the linear chebfun f.
%
% .
%   F.PROP returns the property PROP of F as defined by GET(F, 'PROP').
%
% {}
%   F{S1, S2} restricts F to the subdomain [S1, S2]. See DELTAFUN/RESTRICT for
%   further details. Note that F{[S1, S2]} is not supported due to the behaviour
%   of the MATLAB subsref() command.
%
% See also FEVAL, COMPOSE, GET.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

idx = index(1).subs;
switch index(1).type

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% FEVAL / COMPOSE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case '()'

        % Where to evaluate:
        x = idx{1};
        varin = {};

        % Deal with additional arguments:
        if ( length(idx) == 2 )
            varin = {idx(2)};
        elseif ( length(idx) > 2 )
            error('CHEBFUN:DELTAFUN:subsref:dimensions', ...
                'Index exceeds chebfun dimensions.')
        end

        % Compute the output:
        if ( isnumeric(x) )
            % Call FEVAL():
            out = feval(f, x, varin{:});

        elseif ( isa(x, 'chebfun') )
            % Call COMPOSE():
            % TODO: write compose and check for linearity of f?
            out = compose(x, f);
        else
            error('CHEBFUN:DELTAFUN:subsref:nonnumeric', ...
              'Cannot evaluate chebfun for non-numeric type.')
        end

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GET %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case '.'

        % Call GET() for .PROP access.
        out = get(f, idx);
        if ( numel(index) > 1 )
            % Recurse on SUBSREF():
            index(1) = [];
            out = subsref(out, index);
        end

%% %%%%%%%%%%%%%%%%%%%%%%%%% ACTION of a DISTRIBUTION %%%%%%%%%%%%%%%%%%%%%%%%
    case '{}'
        if ( length(idx) == 1 )
            if ( isequal(idx{1}, ':') )
                % F{:} returns F:
                out = f;
            else
                error('CHEBFUN:DELTAFUN:subsref:badDomain', 'Invalid domain syntax.')
            end

        elseif ( size(idx, 1) == 1 )
            % F{s1,s2,...,sk} returns RESTRICT(F, [s1,s2,...,sk]):
            x = cat(2, idx{:});
            out = restrict(f, x);
        else
            error('CHEBFUN:DELTAFUN:subsref:dimensions', ...
                'Index exceeds chebfun dimensions.')
        end

    otherwise

        error('CHEBFUN:DELTAFUN:subsref:unexpectedType',...
            ['??? Unexpected index.type of ', index(1).type]);
end

% Convert to a cell:
varargout = {out};

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function D = Floyd_Steinberg_Dithering(G)
% ================================================================
% FUNCTION Floyd_Steinberg Dithering Algorithm
%
% Input: G = a 8-bit grayscale / color image
% Output: D = dithered image of G's format with only values 0 and 255 of the same
% ----------------------------------------------------------------
% Demo:
%        G = imread('peppers.png');
%        D = Floyd_Steinberg_Dithering(G);
%        figure('position',[50,50,600,900]),subplot(211),imshow(G),title('Original Image');
%        subplot(212),imshow(D),title('Dithered Image');
% ----------------------------------------------------------------
% For more details about Floyd Steinberg Dithering Algorithm
% Please check
% http://en.wikipedia.org/wiki/Floyd%E2%80%93Steinberg_dithering
% ----------------------------------------------------------------
% Oct. 18th 2011
% By Yue Wu,
% Department of Electrical and Computer Engineering
% Tufts University,
% Medford, MA 02155
% ================================================================

switch size(G,3)
    case 1
    G = double(G); % convert the original image from unit8 to double
    D = zeros(size(G)); % initialize dithered image D
    [M,N] = size(G); % extract size information of the original image
    for r = 1:M
        % applying '2'-like scanning order
        % ->>>>>>>>>>>>>>>>>>>>>>>-
        %                        |
        % -<<<<<<<<<<<<<<<<<<<<<<<-
        % |
        % ->>>>>>>>>>>>>>>>>>>>>>>-
        if mod(r,2) == 0 % scan pixel from left to right
            cOrder = 1:N;
            direction = 'l2r';
        else % scan pixel from right to left
            cOrder = N:-1:1;
            direction = 'r2l';
        end
        for c = cOrder
            tP = G(r,c); % current pixel intensity
            % pick nearest intensity scale two options 0 or 255
            if tP>=128 % close to 255
                D(r,c) = 255; % pick 255
            else % close to 0
                D(r,c) = 0; % pick 0
            end
            eP = tP-D(r,c); % difference before and after selection
            % diffuse difference eP to neighbor pixels
            if r~=M % deal with none bottom rows
                switch direction
                    case 'l2r'
                        if c == 1 % left-most pixel case
                            G(r,c+1) = G(r,c+1)+eP*7/13;
                            G(r+1,c) = G(r+1,c)+eP*5/13;
                            G(r+1,c+1) = G(r+1,c+1)+eP*1/13;
                        elseif c == N % deal with the right-most pixel case
                            G(r+1,c) = G(r+1,c)+eP*5/8;
                            G(r+1,c-1) = G(r+1,c-1)+eP*3/8;
                        else % the normal case
                            G(r,c+1) = G(r,c+1)+eP*7/16;
                            G(r+1,c) = G(r+1,c)+eP*5/16;
                            G(r+1,c+1) = G(r+1,c+1)+eP*1/16;
                            G(r+1,c-1) = G(r+1,c-1)+eP*3/16;
                        end
                    case 'r2l'
                        if c == N % right-most pixel case
                            G(r,c-1) = G(r,c-1)+eP/2;
                            G(r+1,c) = G(r+1,c)+eP*3/8;
                            G(r+1,c-1) = G(r+1,c-1)+eP*1/8;
                        elseif c == 1 % left-most pixel case
                            G(r+1,c) = G(r+1,c)+eP*5/8;
                            G(r+1,c+1) = G(r+1,c+1)+eP*3/8;
                        else % normal case
                            G(r,c-1) = G(r,c-1)+eP*7/16;
                            G(r+1,c) = G(r+1,c)+eP*5/16;
                            G(r+1,c-1) = G(r+1,c-1)+eP*1/16;
                            G(r+1,c+1) = G(r+1,c+1)+eP*3/16;
                        end
                end
            else % deal with the bottom row
                switch direction
                    case 'l2r'
                        if c ~=N % normal case
                            G(r,c+1) = G(r,c+1)+eP;
                        end
                    case 'r2l'
                        if c ~= 1 % normal case
                            G(r,c-1) = G(r,c-1)+eP;
                        end
                end
            end
        end
    end
    otherwise % if the input image is not one-layer gray-scale image, then apply algorithm with respect to each layer
        for i = 1:size(G,3)
            tD = Floyd_Steinberg_Dithering(G(:,:,i));
            D(:,:,i) = tD;
        end
end

D = uint8(D); % convert double D to uint8


text:
classdef MimLevelSlider < MimGuiPluginSlider
    % MimLevelSlider. Gui Plugin for changing overlay opacity
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimLevelSlider is a Gui Plugin for changing the value of the
    %     window used to display the greyscale image
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Level'
        SelectedText = 'Level'
        ToolTip = 'Change the level'
        Category = 'Window / Level'
        Visibility = 'Dataset'
        Mode = 'Toolbar'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Location = 26

        MinValue = 0
        MaxValue = 100
        SmallStep = 0.01
        LargeStep = 0.1
        DefaultValue = 50

        EditBoxPosition = 50
        EditBoxWidth = 40

        StackVertically = true
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && ~strcmp(gui_app.ImagePanel.Mode, MimModes.View3DMode);
        end

        function is_selected = IsSelected(gui_app)
            is_selected = true;
        end

        function [value_instance_handle, value_property_name, limits_instance_handle, limits_property_name] = GetHandleAndProperty(gui_app)
            value_instance_handle = gui_app.ImagePanel.GetBackgroundImageDisplayParameters;
            value_property_name = 'Level';
            limits_instance_handle = gui_app.ImagePanel;
            limits_property_name = 'LevelLimits';
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%%%  ModeSwitch_FLY
clear
path( './icon/',path);
load vehicle_local_position.mat;
%Constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
%throttle when UAV is hovering
THR_HOVER = 0.609;

%% control parameter
%Attitude PID parameters
Kp_PITCH_ANGLE = 6.5;
Kp_PITCH_AngleRate = 0.1;
Ki_PITCH_AngleRate = 0.02;
Kd_PITCH_AngleRate = 0.001;
Kp_ROLL_ANGLE = 6.5;
Kp_ROLL_AngleRate = 0.1;
Ki_ROLL_AngleRate = 0.02;
Kd_ROLL_AngleRate = 0.001;

Kp_YAW_ANGLE = 3;
Kp_YAW_AngleRate = 0.3;
Ki_YAW_AngleRate = 0.001;
Kd_YAW_AngleRate = 0.00;
%Position PID parameters
Kpxp = 1.0;
Kpyp = 1.0;
Kpzp = 4.0;
Kvxp = 2.5; Kvxi = 0.4; Kvxd = 0.01;
Kvyp = 2.5; Kvyi = 0.4; Kvyd = 0.01;
Kvzp = 0.45; Kvzi = 0.01; Kvzd = 0.005;
%integral saturation
Saturation_I_RP_Max = 0.3;
Saturation_I_RP_Min = -0.3;
Saturation_I_Y_Max = 0.2;
Saturation_I_Y_Min = -0.2;
Saturation_I_ah = 3.43;
Saturation_I_az = 5;

%max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
MAX_CONTROL_ANGLE_NAV_ROLL = 15;
MAX_CONTROL_ANGLE_NAV_PITCH  = 15;
%max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
%Maximum control speed, m/s
MAX_CONTROL_VELOCITY_XY = 5;
MAX_CONTROL_VELOCITY_Z = 3;
%Throttle amplitude
MAX_MAN_THR = 0.9;
MIN_MAN_THR = 0.05;
%% run simulink model
ModeSwitch_FLY

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function [epoch, jd] = mydatemjdi (mjd)
    jd = mjd + 2400000.5;
    epoch = mydatejdi (jd);
end

%!test
%! num = mydatenum([1858 11 17  0 0 0]);
%! num2 = mydatemjdi(0);
%! myassert (num2, num)


text:
function [Fitness,I,C] = CalFitness(Population,kappa)
% Calculate the fitness of each solution

%------------------------------- 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 Wenhua Li

    PopObj = Population.objs;
    N = size(PopObj,1);
    PopObj = (PopObj-repmat(min(PopObj),N,1))./(repmat(max(PopObj)-min(PopObj),N,1));
    I      = zeros(N);
    for i = 1 : N
        for j = 1 : N
            I(i,j) = max(PopObj(i,:)-PopObj(j,:));
        end
    end
    C = max(abs(I));
    Fitness = sum(-exp(-I./repmat(C,N,1)/kappa)) + 1;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [retdat,algo] =  training(algo,dat)

retdat=dat;

 disp(['training ' get_name(algo) '.... '])
 switch algo.optimizer
  case {'svmtorch'}
  %%<<----------------svmtorch optimizer-------------------->>

    multi = 0;
    regression = 1;
    degree = 1;
    gamma = 1;
    eps = algo.epsilon;
    C = algo.C;
    if strcmp(algo.child.ker,'linear')
       kernelType = 0;
    elseif strcmp(algo.child.ker,'poly')
       kernelType = 1;
       degree = algo.child.kerparam;
    elseif strcmp(algo.child.ker,'rbf')
    	kernelType = 2;
	gamma = algo.child.kerparam;
    end;
    [x y] = get_xy(dat);
    numEx = get_dim(dat);
    [alpha,threshold0,xSV] = SVMTorch(x,y,regression,multi,kernelType,degree,gamma,C,eps);
    if isempty(alpha)
    	alpha = zeros(numEx,1);
    end;

 %<<--------------sparse svr------------------>>%

   case{'sparse'}
   	kern1=get_kernel(algo.child,dat,dat);
	y=get_y(dat);
	kNum=size(kern1,1);
	c=ones(kNum*4,1);
	c(2*kNum+1:4*kNum)=algo.C;
	kern2 = [-kern1 kern1 -1*eye(kNum) zeros(kNum); kern1 -kern1 zeros(kNum) -1*eye(kNum); -1*eye(4*kNum)];
	cst = zeros(6*kNum,1);
	cst(1:kNum) = (algo.epsilon+algo.b0) - y;
	cst(kNum+1:2*kNum) = (algo.epsilon-algo.b0) + y;
	opts= optimset('display','off','MaxIter',10000,'LargeScale','off');
	[alphaTemp,fval,exit,out,lambda] = linprog(c,kern2,cst,[],[],[],[],[],opts);
	threshold0=lambda.ineqlin(1);
	alpha=alphaTemp(1:kNum)-alphaTemp(kNum+1:2*kNum);

  %<<------------andre optimizer--------------------->>

  case {'andre'}
   kern1=get_kernel(algo.child,dat,dat);   %% <--- calculate kernel
   y=get_y(dat);
   yLen = length(y(:,1));
   xSV = get_x(dat);
   kern2 = [kern1 , -kern1 ; -kern1 , kern1];
   cst = ones(2*yLen,1);
   cst(yLen+1:2*yLen) = -1;
   if algo.nu ==0,
       c = zeros(2*yLen,1);
       c(1:yLen) = algo.epsilon*ones(yLen,1) - y;
       c(yLen+1:2*yLen) = algo.epsilon*ones(yLen,1) + y;
       [alphaTemp,threshold] = quadsolve(kern2,c,cst',0,algo.C);
       alpha=alphaTemp(1:yLen)-alphaTemp(yLen+1:2*yLen);
       threshold0 = -threshold;
   else
     if algo.C==Inf,
         algo.C=10000;
     end;

        c = zeros(2*yLen+1,1);
        c(1:yLen) = -y;
        c(yLen+1:2*yLen) = y;
        cst2 = ones(1,2*yLen+1)/(yLen*algo.nu);
        cst2(2*yLen+1)=-1;
        kern2 = [kern2,zeros(2*yLen,1);zeros(1,2*yLen+1)];
        cst=[cst',0;cst2];
        [alphaTemp,threshold] = quadsolve(kern2,c,cst,[0;0],algo.C); threshold0 = -threshold(1);
        alpha=alphaTemp(1:yLen)-alphaTemp(yLen+1:2*yLen);
        epsilon = -threshold(2);
    end

  %<<------------quadprog optimizer--------------------->>

  case {'quadprog'}

   kern1=get_kernel(algo.child,dat,[]);   %% calculate kernel
   kNum=size(kern1,1);  y=get_y(dat);
   kern2 = [kern1 , -kern1 ; -kern1 , kern1];
   c = zeros(2*kNum,1);
   c(1:kNum) = algo.epsilon*ones(kNum,1) - y;
   c(kNum+1:2*kNum) = algo.epsilon*ones(kNum,1) + y;
   cst(kNum+1:2*kNum) = -1;
   cst = ones(2*kNum,1);
   opts= optimset('display','off','MaxIter',10000,'LargeScale','off');
   [alphaTemp,fval,exit,out,lambda] = quadprog(kern2,c,[],[],cst',0,...
				       zeros(2*kNum,1),algo.C*ones(2*kNum,1),[],opts);
   threshold0=lambda.eqlin(1);
   alpha=alphaTemp(1:kNum)-alphaTemp(kNum+1:2*kNum);

%<<------------libsvm optimizer--------------------->>
case {'libsvm'}

        %
        x=[];
        y=[];
        svm_type=3;
        kernelType=0;
        degree=3;
        gamma=0;
        coef0=0;

        nu=algo.nu;
        if(nu>0)
            svm_type=4;
        end

        cachesize=40;
        C=algo.C;
        eps=algo.epsilon;
        p=0.05;
        shrinking=1;

        weight_label=[];
        weight=[];
        nr_weight=0;

        if strcmp(algo.child.ker,'linear')
            kernelType = 0;
        end;
        if strcmp(algo.child.ker,'poly')
            kernelType = 1;
            degree = algo.child.kerparam;
            coef0 = 1;
            gamma = 1;
        end;
        if strcmp(algo.child.ker,'rbf'),
            kernelType = 2;
            sigma = algo.child.kerparam;
            gamma = 1/(2*sigma^2);
        end;

        y=get_y(dat);
        x=get_x(dat);

        if strcmp(algo.child.ker,'custom'),
          kernelType = 4;
          K= algo.child.kerparam;
          l = get_dim( retDat);
          x = get_index( retDat);
          x = [ reshape( x, l, 1) [ 1:l]']; % using x to pass indices in Matrix and real indices
        end;

        s=whos('libsvm_cachesize','global');

        if (length(s)>0)
            global libsvm_cachesize;
            cachesize=libsvm_cachesize;
        else
            cachesize=40;
        end
        if algo.algorithm.verbosity>1
         fprintf('Using %d MB Cache for Libsvm\n',cachesize)
        end

         if( kernelType == 4)
          [alpha,xSV,bias0]=libsvm_regressor_spider(x,y,svm_type,kernelType,...
                     degree,gamma,coef0,nu,cachesize,C,eps,p,weight_label,weight,nr_weight,K);
          algo.Xsv=get(retDat, xSV( :, 2));
        else
          [alpha,xSV,bias0]=libsvm_regressor_spider(x,y,svm_type,kernelType,...
                     degree,gamma,coef0,nu,cachesize,C,eps,p,weight_label,weight,nr_weight);
          algo.Xsv=data(xSV);
        end

        threshold0 = bias0 * y(1);

        algo.b0=bias0;

        algo.epsilon = eps;

        algo.Xsv = data(xSV);
        algo.alpha=alpha;

        if algo.algorithm.do_not_evaluate_training_error==1
            retdat=set_x(dat,get_y(dat));
        else
            retdat=test(algo,dat);
        end

        return

        %         fin=find(abs(alpha)>algo.alpha_cutoff);

% case {'libsvm'}
%
%    if algo.nu ==0,
%         svm_type = 3;
%         C = algo.C;
%         epsilon = algo.epsilon;
%         nu=0;
%     else
%         svm_type = 4;
%         nu = algo.nu;
%         C = algo.C;
%         epsilon = -1;
%     end;
%     %% default values for libsvm
%     cacheSize = 40;
%     eps = 0.001;
%     shrinking=1;
%     nrWeight = 0;
%     weightLabel =0;
%     weight = 1;
%     gamma=1;
%     deg = 0;
%     coef0 = 0;
%     kerTmp = algo.child;
%     if strcmp(kerTmp.ker,'linear')
%       kernelType = 0;
%     end;
%     if strcmp(kerTmp.ker,'poly')
%       kernelType = 1;
%       ptmp = kerTmp.kerparam;
%       deg = ptmp;
%       coef0 = 1;
%     end;
%    if strcmp(kerTmp.ker,'rbf'),
%          kernelType = 2;
%          ptmp = kerTmp.kerparam;
%          gamma = 1/(2*ptmp^2);
%    end;
% %    if algo.balanced_ridge~=0,
% %         disp('Warning: balanced ridge not implemented for libsvm.');
% %    end;
%    y=get_y(dat);
%    x=get_x(dat);
%    [alpha,threshold0,xSV,eps,CC] = svmlibtrain(x,y,svm_type,kernelType,deg,gamma,coef0,nu,cacheSize,C,eps,epsilon,...
%        shrinking,nrWeight,weightLabel,weight,0);
%    threshold0=-threshold0;% from libsvm
%    % alpha is reordered in order to have the same xsp for all runs (important for one_vs_rest)
%    alphaTemp = zeros(size(xSV,1),1);
%    indTemp = find(xSV(:,size(xSV,2))~=0);
%    indTemp2 = xSV(indTemp,size(xSV,2));
%    alphaTemp(indTemp2) = alpha(indTemp);
%    alpha = alphaTemp;
%    epsilon=eps;
%  end
end

 if algo.nu~=0,
   algo.epsilon = epsilon;
 end;

 algo.b0=threshold0;
 fin=find(abs(alpha)>algo.alpha_cutoff);
 algo.alpha=alpha(fin);
 algo.Xsv = get(dat,fin);

 if algo.algorithm.do_not_evaluate_training_error==1
   retdat=set_x(dat,get_y(dat));
 else
   retdat=test(algo,dat);
 end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function docking = ma_createDocking(name, oScId, undockTime, lineColor, lineStyle, lineWidth, massLoss)
%ma_createDocking Summary of this function goes here
%   Detailed explanation goes here

    docking = struct();

    docking.id          = rand(1);
    docking.type        = 'Docking';
    docking.name        = name;
    docking.oScId       = oScId;
    docking.undockTime	= undockTime;
    docking.lineColor   = lineColor;
    docking.lineStyle   = lineStyle;
    docking.lineWidth   = lineWidth;
    docking.massloss	= massLoss;
end


text:
function q = quaternion(sym)

q = quaternion(sym.properGroup.rot);

end

author: mtex-toolbox

repo: mtex

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


text:
function plterrel(xo,yo,C,s,ltype)
% PLTERREL  Plots an error ellipse on screen.
%   Note 1: x & y represent north & south (opposite of
%   normal MatLab convention).  Note 2: use a square
%   aspect ratio for plot; i.e., use axis('square')
%   command.  Non-vectorized.
% Version: 18 Jan 96
% Useage:  plterrel(xo,yo,C,s,ltype)
% Input:   xo - x (north) origin of ellipse
%          yo - y (east) origin of ellipse
%          C  - covariance matrix (2x2)
%          s  - ellipse scale factor (default=1)
%          ltype - line type for error ellipse
%                  (default=solid,red)
% Output:  Plotted error ellipse centred at (xo,yo)

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

if nargin<3
  error('Too few input arguements');
end
if nargin<4
  s=1;             % Define default scale factor
end
if nargin<5
  ltype='-r';      % Define default color red
end
dt=0.1;            % Angular resolution
t=[(0:dt:2*pi)';0];
[V,D]=eig(C);      % Eigenvalues & vectors
r=sqrt(diag(D));   % Length of axes
x=s*r(1)*cos(t);
y=s*r(2)*sin(t);
xx=xo+[x y]*V(1,:)';
yy=yo+[x y]*V(2,:)';
plot(yy,xx,ltype);  % Switch x (north) & y (east)


text:
function test_suite = test_projPointOnLine
%TESTPROJPOINTONLINE  One-line description here, please.
%   output = testProjPointOnLine(input)
%
%   Example
%   testProjPointOnLine
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testHorizontal(testCase) %#ok<*DEFNU>
point = [0 0];
line = [1 0 0 1];
testCase.assertEqual(projPointOnLine(point, line), [1 0], 'AbsTol', .01);

function testDiagonal(testCase)
point = [0 0];
line = [2 0 1 1];
testCase.assertEqual(projPointOnLine(point, line), [1 -1], 'AbsTol', .01);

function testBigDerivative(testCase)
point = [0 0];
line = [2 0 1000 1000];
testCase.assertEqual(projPointOnLine(point, line), [1 -1], 'AbsTol', .01);

function testDiagonal2(testCase)
point = [2 3];
line = [-2 -4 6 4];
testCase.assertEqual(projPointOnLine(point, line), [4 0], 'AbsTol', .01);

function test_SingleMulti(testCase)

point = [4 2];
line1 = [0 0 1 0];
line2 = [0 0 0 1];
line3 = [0 2 2 2];

res = projPointOnLine(point, [line1;line2;line3]);
exp = [4 0;0 2;2 4];
testCase.assertEqual(exp, res, 'AbsTol', .01);

function test_MultiSingle(testCase)

line = [0 2 4 2];
p1 = [3 1];
p2 = [5 2];
p3 = [3 6];

res = projPointOnLine([p1;p2;p3], line);
exp = [2 3;4 4;4 4];
testCase.assertEqual(exp, res, 'AbsTol', .01);

function test_MultiMulti(testCase)

line1 = [0 0 1 0];
line2 = [0 0 0 1];
line3 = [0 2 2 2];
p1 = [3 1];
p2 = [2 3];
p3 = [1 5];
res = projPointOnLine([p1;p2;p3], [line1;line2;line3]);
exp = [3 0;0 3;2 4];
testCase.assertEqual(exp, res, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% -------------------------------------------------------------------------
clear all; close all;

%connection settings
number_of_retries = 1; % set to -1 for infinite
port=2057; % This code use for communication "port" and "port+1"
host='localhost';
samplingrate=1000; %in ms

while(true)
    data_to_send=['|',regexprep(num2str(randn(1),6),'\.',','),'|',regexprep(num2str(randn(1),6),'\.',',')];
    data_received=exchangeData(port,host,number_of_retries,samplingrate,data_to_send)
    if(numel(data_received)>0 )
        % Do your stuff with the data
    end;
end


text:
%TestCaseTest Unit tests for the TestCase class

%   Steven L. Eddins
%   Copyright The MathWorks 2008

classdef TestCaseTest < TestCaseInDir

    methods
        function self = TestCaseTest(name)
            self = self@TestCaseInDir(name, ...
                fullfile(fileparts(which(mfilename)), 'helper_classes'));
        end

        function testConstructor(self)
            % Exercise the constructor.  Verify that the Name and Location
            % properties are set correctly.
            tc = TwoPassingTests('testMethod1');
            assertEqual(tc.Name, 'testMethod1');
            assertEqual(tc.Location, which('TwoPassingTests'));
        end

        function testPassingTests(self)
            % Verify that the expected observer notifications are received in
            % the proper order.
            logger = TestRunLogger();
            TestSuite('TwoPassingTests').run(logger);
            assertTrue(isequal(logger.Log, ...
                {'TestRunStarted', 'TestComponentStarted', ...
                'TestComponentStarted', 'TestComponentFinished', ...
                'TestComponentStarted', 'TestComponentFinished', ...
                'TestComponentFinished', 'TestRunFinished'}));
        end

        function testFixtureCalls(self)
            % Verify that fixture calls are made in the proper order.
            tc = LoggingTestCase('testMethod');
            tc.run(TestRunLogger());
            assertTrue(isequal(tc.log, {'setUp', 'testMethod', 'tearDown'}));
        end

        function testTestFailure(self)
            % Verify that a test failure is recorded.
            logger = TestRunLogger();
            TestSuite('FailingTestCase').run(logger);
            assertTrue(isequal(logger.NumFailures, 1));
        end

        function testTestError(self)
            % Verify that a test error is recorded.
            logger = TestRunLogger();
            TestSuite('BadFixture').run(logger);
            assertTrue(isequal(logger.NumErrors, 1));
        end

    end

end


text:
function comp_unrank_grlex_test ( )

%*****************************************************************************80
%
%% COMP_UNRANK_GRLEX_TEST tests COMP_UNRANK_GRLEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 October 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMP_UNRANK_GRLEX_TEST\n' );
  fprintf ( 1, '  A COMP is a composition of an integer N into K parts.\n' );
  fprintf ( 1, '  Each part is nonnegative.  The order matters.\n' );
  fprintf ( 1, '  COMP_UNRANK_GRLEX determines the parts\n' );
  fprintf ( 1, '  of a COMP from its rank.\n' );

  kc = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Rank: ->  NC       COMP\n' );
  fprintf ( 1, '  ----:     --   ------------\n' );

  for rank1 = 1 : 71

    xc = comp_unrank_grlex ( kc, rank1 );
    nc = sum ( xc(1:kc) );

    fprintf ( 1, '   %3d: ', rank1 );
    fprintf ( 1, '    %2d = ', nc );
    for j = 1 : kc - 1
      fprintf ( 1, '%2d + ', xc(j) );
    end
    fprintf ( 1, '%2d\n', xc(kc) );
%
%  When XC(1) == NC, we have completed the compositions associated with
%  a particular integer, and are about to advance to the next integer.
%
    if ( xc(1) == nc )
      fprintf ( 1, '  ----:     --   ------------\n' );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mode = addReservedColors(mode)
%
% mode = addReservedColors(mode)
%
% Adds reserved colors to a displayMode and its underlying colormap.
%
% djh 1/98

mode.white=256;
mode.cmap(256,:) = [1,1,1]; 		% white

mode.black=255;
mode.cmap(255,:) = [0,0,0]; 		% black

mode.red=254;
mode.cmap(254,:) = [1,0,0]; 		% red

mode.green=253;
mode.cmap(253,:) = [0,1,0]; 		% green

mode.blue=252;
mode.cmap(252,:) = [0,0,1]; 		% blue

mode.cyan=251;
mode.cmap(251,:) = [0,1,1]; 		% cyan

mode.magenta=250;
mode.cmap(250,:) = [1,0,1]; 		% magenta

mode.yellow=249;
mode.cmap(249,:) = [1,1,0]; 		% yellow

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%Detection the blueprint for Detection objects in GISMO
% An Detection object is a container for sta/lta detection metadata
% See also Arrival, Catalog
classdef Detection
    properties
        channelinfo
        time
        state
        filterString
        signal2noise
        traveltime
    end
    properties(Dependent)
        numel
    end
    methods
        function obj = Detection(sta, chan, time, state, filterString, signal2noise)

            % Blank constructor
            if ~exist('sta','var')
                return
            end

            % Parse required, optional and param-value pair arguments,
            % set default values, and add validation conditions
            p = inputParser;
            p.addRequired('sta', @iscell);
            p.addRequired('chan', @iscell);
            %p.addRequired('time', @(t) t>0 & t<now+1);
            p.addRequired('time', @isnumeric);
            p.addOptional('state', {}, @iscell);
            p.addOptional('filterString', {}, @iscell);
            p.addOptional('signal2noise', [], @isnumeric);

            % Missed several properties out here just because of laziness.
            % Add them as needed.
            p.parse(sta, chan, time, state, filterString, signal2noise);
            ctag = ChannelTag.array('',p.Results.sta,'',p.Results.chan)';
            obj.channelinfo = ctag.string();
            obj.time = p.Results.time;
            obj.state = p.Results.state;
            obj.filterString = p.Results.filterString;
            obj.signal2noise = p.Results.signal2noise;
            N = numel(obj.time);
            fprintf('\nGot %d detections\n',N);
            if numel(obj.state) == 0
                obj.state = repmat({''},1,N);
            end
            if numel(obj.filterString) == 0
                obj.filterString = repmat({''},1,N);
            end
            if numel(obj.signal2noise) == 0
                obj.signal2noise = repmat({''},1,N);
            end
            if numel(obj.traveltime) == 0
                obj.traveltime = repmat(NaN,1,N);
            end

        end

        function val = get.time(obj)
            val = obj.time;
        end

        function val = get.numel(obj)
            val = numel(obj.time);
        end

%         function val = get.channelinfo(obj)
%             val = obj.channelinfo;
%         end
%
%         function val = get.iphase(obj)
%             val = obj.iphase;
%         end
%
%         function val = get.amp(obj)
%             val = obj.amp;
%         end
%
%         function val = get.signal2noise(obj)
%             val = obj.signal2noise;
%         end
%
%         function obj = set.amp(obj, amp)
%             obj.amp = amp;
%         end

        function plot(obj)
            ctaguniq = unique(obj.channelinfo)
            N = numel(ctaguniq);
            hf1=figure;
            suptitle('Cumulative detections vs time')
            hf2=figure;
            suptitle('Percentage of detections captured by signal2noise')
            disp(sprintf('Detections: %d',numel(obj.time)))
            for c=1:N
               indexes = find(strcmp(obj.channelinfo,ctaguniq{c})==1);
%                size(indexes)
%                indexes(1:10)
               figure(hf1)
               %subplot(N,1,c);
               hold on
               t = obj.time(indexes);
               y = cumsum(ones(size(t)));
plot(t,y,'LineWidth',5);
               %ylabel(ctaguniq{c});
               ylabel(sprintf('Cumulative #\nDetections'))
               xlabel('Date/Time')
               datetick('x')
               set(gca,'XLim',[min(obj.time) max(obj.time)]);
               disp(sprintf('- %s: %d',ctaguniq{c},length(indexes==1)));

               figure(hf2)
               %subplot(N,1,c);
               hold on
               s = obj.signal2noise(indexes);
               s(s>101)=101;
               [n x]=hist(s, 1:0.1:100);
plot(x,100-cumsum(n)/sum(n)*100,'LineWidth',5);
               xlabel('signal2noise');
               %ylabel(ctaguniq{c});
               ylabel('%age')
               set(gca, 'XLim', [2.4 20]);
            end
            figure(hf1)
            legend(ctaguniq)
            figure(hf2)
            legend(ctaguniq)
        end

        function summary(obj, showall)
        % DETECTION.SUMMARY Summarise Detection object
            for c=1:numel(obj)
                obj(c)
                numrows = numel(obj(c).time);
                fprintf('Number of detections: %d\n',numrows);
                if numrows > 0
                    if ~exist('showall','var')
                            showall = false;
                    end
                    if numel(obj) == 1
                        if numrows <= 50 || showall
                            for rownum=1:numrows
                                summarize_row(obj, rownum);
                            end
                        else
                            for rownum=1:50
                                summarize_row(obj, rownum);
                            end
                            disp('* Only showing first 50 rows/arrivals - to see all rows/arrivals use:')
                            disp('*      arrivalObject.disp(true)')
                        end
                    end
                end
            end
        end

        function summarize_row(self, rownum)
            fprintf('%s\t%s\t%s\t%e\n', ...
                self.channelinfo(rownum), ...
                datestr(self.time(rownum)), ...
                self.state(rownum), ...
                self.signal2noise(rownum));
        end

        function self2 = subset(self, columnname, findval)
            self2 = self;
            N = numel(self.time);
            indexes = [];
            if ~exist('findval','var')
                % assume columnname is actually row numbers
                indexes = columnname;
            else

                for c=1:N
                    gotval = eval(sprintf('self.%s(c);',columnname));
                    if isa(gotval,'cell')
                        gotval = cell2mat(gotval);
                    end
                    if isnumeric(gotval)
                        if gotval==findval
                            indexes = [indexes c];
                        end
                    else
                        if strcmp(gotval,findval)
                            indexes = [indexes c];
                        end
                    end
                end
            end
%             if isa(self.channelinfo,'cell')
%                 self2.channelinfo = self.channelinfo{indexes};
%             else
                 self2.channelinfo = self.channelinfo(indexes);
%             end
            self2.time = self.time(indexes);
            self2.state = self.state(indexes);
            if numel(self.filterString)==N
                self2.filterString = self.filterString(indexes);
            end
            if numel(self.signal2noise)==N
                self2.signal2noise = self.signal2noise(indexes);
            end
            if numel(self.traveltime)==N
                self2.traveltime = self.traveltime(indexes);
            end
        end

        function self = append(self1, self2)
            disp('Appending...')
            [newtime, indices] = sort([self1.time self2.time]);
            size(cellstr(self1.channelinfo))
            size(cellstr(self2.channelinfo))
            newchannelinfo = [cellstr(self1.channelinfo); cellstr(self2.channelinfo)];
            newchannelinfo = newchannelinfo(indices);
            newstate = [cellstr(self1.state); cellstr(self2.state)];
            newstate = newstate(indices);
            newfs = [cellstr(self1.filterString); cellstr(self2.filterString)];
            newfs = newfs(indices);
            newsnr = [self1.signal2noise self2.signal2noise];
            newsnr = newsnr(indices);
            ctags = ChannelTag(newchannelinfo);
            %self = Detection(cellstr([get(ctags,'station')]), cellstr([get(ctags,'channel')]), newtime, cellstr([newstate]), cellstr([newfs]), newsnr)
            self = Detection([get(ctags,'station')], ...
                [get(ctags,'channel')], ...
                newtime, ...
                newstate, ...
                [newfs], ...
                newsnr)
        end

        % prototypes
         catalogobj = associate(self, maxTimeDiff, sites, source)
%         write(detectionobj, FORMAT, PATH);
    end
    methods(Static)
        function self = retrieve(dbname, subset_expr)
        %DETECTION.RETRIEVE Read detections from an Antelope CSS3.0 table
        %
        % Usage:
        %       detectionObj = DETECTION.RETRIEVE(dbname, subset_expr)

        % Author: Glenn Thompson (glennthompson1971@gmail.com)

            debug.printfunctionstack('>')
            self = [];
            if ~(antelope.dbtable_present(dbname, 'detection'))
                fprintf('No detection table belonging to %s\n',dbname);
                return
            end

            fprintf('Loading detections from %s\n',dbname);

            % Open database
            db = dbopen(dbname,'r');
            disp('- database opened');

            % Apply subset expression
            db = dblookup_table(db,'detection');
            disp('- detection table opened');
            if exist('subset_expr','var')
                db = dbsubset(db,subset_expr);
                disp('- subsetted database')
            end

            nrows = dbnrecs(db);
            if nrows > 0

                % Sort by arrival time
                db = dbsort(db,'time');
                disp('- sorted detection table')

                % Get the values
                fprintf('- reading %d rows\n',nrows);
                [sta,chan,time,state,filterString,signal2noise] = dbgetv(db,'sta','chan','time','state', 'filter','snr');

                % Close database link
                dbclose(db);
                disp('- database closed')

                % Create detection object
                disp('- creating detection object')
                self = Detection(cellstr(sta), cellstr(chan), epoch2datenum(time), cellstr(state), cellstr(filterString), signal2noise);

                disp('- complete!')
            else
                fprintf('No detections found matching request\n');
            end

            debug.printfunctionstack('<')
        end

        %cookbook()
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function [inerr,dw,inLL,right ]=batch_equal_nomask_lstm(W,dataall  )
% gradient threshold for bptt need done   ,inLL,right
global   mzeros convert in_size gate_size out_size share_size usegpu...
    share_size2  in ingate cellstate cells outgate globalData    ...
    node_outgateInit cellinInit node_cellbiasInit delta_outInit
if nargin==1
    data=globalData;
else
    data=dataall;
end
if usegpu
    fun=@sigmoidnGpu;
    delta_fun=@delta_sigmoidnGpu;
    activation=@activationGpu;
    deactivation=@deactivationGpu;
else
    fun=@sigmoidnCpu;
    delta_fun=@delta_sigmoidnCpu;
    activation=@activationCpu;
    deactivation=@deactivationCpu;
end

actNum1=convert(1);
actNum2=2*actNum1;
numsamples=size(data,1) ;

dw= zeros(1 ,2*share_size2 +share_size2*numMmcell + (gate_size*numMmcell+1)*out_size);
[ Wingate,Wcell , Woutgate, Wout]=unpack(W);

[dwingate,dwcell,dwoutgate,dwout] = unpack2(dw);
dwingate=convert(dwingate);
dwcell=convert(dwcell);
dwoutgate=convert(dwoutgate);
dwout=convert(dwout);

right=zeros(numsamples,1,'int32');
inLL = 0;
inerr = zeros( numsamples,out_size);

timespan = size( data  , 3 );
node_tmpindex=cell(1,timespan);
node_tmpindex2=cell(1,timespan);
node_tmpindex3=cell(1,timespan);

node_outgate= node_outgateInit;
cellin=cellinInit;
node_cellbias=node_cellbiasInit;
delta_out=delta_outInit;
cellstatus= cellstatusInit;

node_ingate = node_outgate;
node_cell = cellin;
Y_cellout = cellin;
delta_outgate=node_outgate;

errorstate=cellin;
delta_cellin=cellin;
delta_ingate=node_outgate;

Wingate_in=convert( Wingate(in,:) );
Wingate_ingate= convert( Wingate(ingate,:) );
Wingate_cellstate=convert( Wingate(cellstate,:) );
Wingate_cell=convert( Wingate(cells,:) );
Wingate_outgate=convert( Wingate(outgate,:) );

% to cellin
Wcell_in=convert( Wcell(in,:) ) ;
Wcell_ingate=convert( Wcell(ingate,:) );
Wcell_cell=convert( Wcell(cells,:) );
Wcell_outgate=convert( Wcell(outgate,:) );

% to outgate
Woutgate_in=convert( Woutgate(in,:) );
Woutgate_ingate=convert( Woutgate(ingate,:) );
Woutgate_cellstate=convert( Woutgate(cellstate,:) );
Woutgate_cell = convert( Woutgate(cells,:) );
Woutgate_outgate=convert( Woutgate(outgate,:) );

t=2;
while t<=timespan

    output = data(:,1:out_size,t ) ;

    % forward pass
    node_in = data(:,1:in_size,t-1);
    node_cellbias{t } = data(:,in_size+1,t-1);

    % to ingate
    tmpp=node_in*Wingate_in + node_outgate{t-1}*Wingate_outgate +...
        cellstatus{t-1} * Wingate_cellstate + node_cell{t-1}*Wingate_cell +  node_ingate{t-1} * Wingate_ingate;
    node_ingate{t}=fun( tmpp ,signum);

    tmpp=node_in*Wcell_in +  node_ingate{t-1}  * Wcell_ingate  +...
        + node_cell{t-1}*Wcell_cell +node_outgate{t-1} * Wcell_outgate ;
    cellin{t} = activation( tmpp,actNum2 ) ;

    cellstatus{t}  =  cellstatus{t-1}  +    cellin{t }  .* repmat(node_ingate{t},1,numMmcell )  ;

    Y_cellout{t}= activation( cellstatus{t } ,actNum1 ) ;

    tmpp=node_in*Woutgate_in +   node_outgate{t-1} *Woutgate_outgate +...
        cellstatus{t}  * Woutgate_cellstate + node_cell{t-1}*Woutgate_cell +  node_ingate{t-1} * Woutgate_ingate;
    node_outgate{t }  = fun(tmpp  ,signum);

    node_cell{t }  = repmat(node_outgate{t},1,numMmcell )   .* Y_cellout{t }  ;

    node_tmpindex{t } = [ node_in  node_ingate{t-1} cellstatus{t-1} node_cell{t-1} node_outgate{t-1} ];
    node_tmpindex2{t} = [ node_in  node_ingate{t-1} cellstatus{t} node_cell{t-1} node_outgate{t-1} ];
    node_tmpindex3{t } = [ node_in  node_ingate{t-1}  node_cell{t-1} node_outgate{t-1} ];

    node_out  = fun(  [ node_cell{t}     node_cellbias{t}   ] * Wout    ,signum ) ;

    delta_out{t} =  ( - output + node_out  )  .* delta_fun(node_out,signum);

    inerr = inerr +     (  (output  - node_out  )).^2    ;
    %       right = right +   rightfun(masko,output,node_out) ;
    t=t+1;
end
t=t-1;
% just  cell
tmp_outgate=repmat(node_outgate{t},1,numMmcell);
tmp_ingate=repmat(node_ingate{t},1,numMmcell);
errorcell = delta_out{t} * Wout(1:end-1,:)';

% just cell to outgate
delta_outgate{t} =  squeezing(delta_fun( tmp_outgate,signum) .* errorcell  .* Y_cellout{t})  ;  %  .* nodenew_outgate.* ( 1 - nodenew_outgate  )  ;
% peephole
errorstate{t} = errorcell .* tmp_outgate .* deactivation(Y_cellout{t},actNum1)+...
    delta_outgate{t} * Woutgate_cellstate' ;
delta_cellin{t} =tmp_ingate .* errorstate{t}.* deactivation(cellin{t},actNum2);
delta_ingate{t} =  squeezing(cellin{t} .* errorstate{t} .* delta_fun(tmp_ingate,signum));
for t = timespan-1:-1:2
    % just  cell
    tmp_outgate=repmat(node_outgate{t},1,numMmcell);
    tmp_ingate=repmat(node_ingate{t},1,numMmcell);
    errorcell = delta_out{t } * Wout(1:end-1,:)' + delta_outgate{t+1} * Woutgate_cell' +...
        delta_cellin{t+1} * Wcell_cell' + delta_ingate{t+1} * Wingate_cell'; % W( out , %cell);

    % just cell to outgate
    delta_outgate{t} = delta_fun(  node_outgate{t},signum )  .* ( squeezing(errorcell  .* Y_cellout{t})  +...
        delta_ingate{t+1} * Wingate_outgate'+ delta_cellin{t+1} * Wcell_outgate'+...
        delta_outgate{t+1} *Woutgate_outgate') ;  %  .* nodenew_outgate.* ( 1 - nodenew_outgate  )  ;
    % peephole
    errorstate{t} = errorcell .* tmp_outgate .* deactivation( Y_cellout{t} ,actNum1 )  + ...
        errorstate{t+1} + delta_ingate{t+1}* Wingate_cellstate' +...
        delta_outgate{t} * Woutgate_cellstate' ;

    delta_cellin{t} = tmp_ingate .* deactivation( cellin{t } ,actNum2) .* errorstate{t} ;

    delta_ingate{t} = delta_fun( node_ingate{t} ,signum).* ( squeezing(cellin{t} .* errorstate{t})  +...
        + delta_ingate{t+1} * Wingate_ingate' +delta_cellin{t+1} * Wcell_ingate' + ...
        delta_outgate{t+1} * Woutgate_ingate' );
end

for t = 2:timespan
    dwout = dwout + [node_cell{t} node_cellbias{t} ]' *  delta_out{t} ;
    dwoutgate =dwoutgate + node_tmpindex2{t}' * delta_outgate{t} ;
    dwcell=dwcell+ node_tmpindex3{t}'*delta_cellin{t};
    dwingate = dwingate + node_tmpindex{t}'  * delta_ingate{t}   ;
end

inerr = gather(  1/2*    sum(    inerr)   /numsamples)  ;%  1/2*  for gradient checking
%right=gather(sum(right));

dw=pack2(dwingate,dwcell ,dwoutgate,dwout);
dw = gather(dw / numsamples );

    function [Wingate,Wcell ,Woutgate,Wout] = unpack(W)
        Wingate= reshape( W( 1:share_size2 )   , share_size ,gate_size );
        Wcell  = reshape( W( share_size2 +1 :   share_size2 *(1+numMmcell) ),  share_size, numMmcell*gate_size);
        Woutgate = reshape( W( (1+numMmcell) * share_size2 +1 : (2+numMmcell)* share_size2  ), share_size , gate_size);
        Wout = reshape( W( (2+numMmcell) * share_size2 +1 : end ) ,gate_size*numMmcell+1 ,  out_size);
    end
    function [Wingate,Wcell ,Woutgate,Wout] = unpack2(W)
        Wingate= reshape( W( 1:share_size2 )   , share_size ,gate_size );
        % Wcell  = reshape( W( share_size2 +1 : 2* share_size2  ),  share_size,gate_size);
        Wcell  = reshape( W( share_size2 +1 :   share_size2 *(1+numMmcell) ),  share_size, numMmcell*gate_size);
        %   share_size = in_size + gate_size * (2+2*problem.numMmcell) ;
        Wcell =  Wcell([1:in_size+gate_size (in_size+gate_size+gate_size*numMmcell+1):end],:) ;
        Woutgate = reshape( W( (1+numMmcell) * share_size2 +1 : (2+numMmcell)* share_size2  ), share_size , gate_size);
        Wout = reshape( W( (2+numMmcell) * share_size2 +1 : end ) ,gate_size*numMmcell+1 ,  out_size);
    end
    function [W]=pack2(Wingate,Wcell , Woutgate,Wout)
        tmp=mzeros(share_size,gate_size*numMmcell);
        tmp([1:in_size+gate_size (in_size+gate_size+gate_size*numMmcell+1):end],:)=Wcell;
        Wcell=tmp;
        W = [Wingate(:);Wcell(:) ;Woutgate(:);Wout(:)];
    end

    function y=sigmoidnGpu(x,num)
        y=arrayfun(@(x,num)1./(1+exp(- num*x)),x,num);
    end
    function y=delta_sigmoidnGpu(x,num)
        y=arrayfun(@(x,num)num*x.*(1-x),x,num);
    end
    function y=sigmoidnCpu(x,num)
        y=1./(1+exp(- num*x));
    end
    function y=delta_sigmoidnCpu(x,num)
        y=  num*x.*(1-x);
    end
    function y =  activationCpu(x,num)
        y=num*2./(1+exp(-x))-num;
    end
    function y =  deactivationCpu(x,num)
        y=0.5/num*(num+x).*(num-x);
    end
    function y =  activationGpu(x,num)
        y= arrayfun(@(x,num)num*2./(1+exp(-x))-num,x,num);
    end
    function y =  deactivationGpu(x,num)
        y=arrayfun(@(x,num)0.5/num*(num+x).*(num-x),x,num);
    end
    function y = squeezing(x)
        y=zeros(size(x,1),gate_size);
        for i = 1:numMmcell
            y = y + x(:,1+(i-1)*gate_size: i*gate_size);
        end
        %         x=mat2cell(x,size(x,1),gate_size*ones(1,numMmcell));
        %
        %         for i=2:numMmcell
        %             x{1}=x{1}+x{i};
        %         end
        %         y=x{1};
        %    y= cell2mat(cellfun(@plus  ,x, 'UniformOutput',false));
    end

end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function clrIX = MapXto1Dcolormap(X,Xrange,numC)
% note: Xrange only saturates high and low values, does not remove data
% points in X.

if nargin < 2
    Xrange = [min(X),max(X)];
end
if nargin < 3
    numC = 64;
end

X(X<Xrange(1)) = Xrange(1);
X(X>Xrange(2)) = Xrange(2);

clrIX = round((X-Xrange(1))/(Xrange(2)-Xrange(1))*(numC-1))+1;
if size(clrIX,1)<size(clrIX,2)
    clrIX = clrIX';
end
end

author: xiuyechen

repo: FishExplorer


text:
function varargout = norm_nuclear(varargin)
% NORM_NUCLEAR Returns sum of singular values
%
% s = SUMK(X,k)
%
% For a vector X, NORM_NUCLEAR returns the sum of singular values
%
% For a matrix X, NORM_NUCLEAR returns the sum of absolute values.
%
% See also SUMABSK

switch class(varargin{1})

    case 'double' % What is the numerical value of this argument (needed for displays etc)
        varargout{1} = sum(svd(varargin{1}));

    case 'sdpvar'
        varargout{1} = yalmip('define',mfilename,varargin{:});

    case 'char' % YALMIP sends 'model' when it wants the epigraph or hypograph
        if isequal(varargin{1},'graph')
            t = varargin{2};
            X = varargin{3};
            [n,m] = size(X);
            if is(X,'real')
                U = sdpvar(m);
                V = sdpvar(n);
            else
                U = sdpvar(m,m,'hermitian','complex');
                V = sdpvar(n,n,'hermitian','complex');
            end
            F = [trace(U)+trace(V) <= 2*t, [U X';X V]>=0];
            varargout{1} = F;
            varargout{2} = struct('convexity','convex','monotonicity','none','definiteness','positive','model','graph');
            varargout{3} = X;
        else
            varargout{1} = [];
            varargout{2} = [];
            varargout{3} = [];
        end
    otherwise
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP