text:
function findBallFcn(greenBall1)

%% Find Green Object
% This script reads in an image file and then attempts to find a green
% object in the image. It is designed to find one green ball and highlight
% that ball on the original image

% Copyright 2013 The MathWorks, Inc.

%% Step 1: Read image into MATLAB
% First we read the specified image from the file and bring it into MATLAB
% as a variable. We also display the image to ensure it is correct.
% greenBall1 = imread('greenBall3.jpg');
% imtool(greenBall1);

%% Step 2: Identify Unique Characteristics of Object of Interest

%%
% Extract each color
% Next we using indexing to extract three 2D matrices from the 3D image
% data corresponding to the red, green, and blue components of the image.
r = greenBall1(:, :, 1);
g = greenBall1(:, :, 2);
b = greenBall1(:, :, 3);

%%
% View different color planes
% figure
% subplot(2,2,1),imagesc(r)
% subplot(2,2,2),imagesc(g)
% subplot(2,2,3),imagesc(b)

%%
% Calculate Green
% Then we perform an arithmetic operation on the matrices as a whole to try
% to create one matrix that represents an intensity of green.
justGreen = g - r/2 - b/2;
% colorPlanesPlot(r,g,b,justGreen);

%%
% close all

%% Step 3: Isolate Object of Interest

%%
% Threshold the image
% Now we can set a threshold to separate the parts of the image that we
% consider to be green from the rest.
bw = justGreen > 50;
% imagesc(bw);
% colormap(gray);

%%
% Remove small unwanted objects
% We can use special functions provided by the Image Processing toolbox to
% quickly perform common image processing tasks. Here we are using
% BWAREAOPEN to remove groups of pixels less than 30.
ball1 = bwareaopen(bw, 30);
% imagesc(ball1);

%% Step 4: Find center of green object
% Now we are using REGIONPROPS to extract the centroid of the group of
% pixels representing the ball.
s  = regionprops(ball1, {'centroid','area'});
if isempty(s)
  error('No ball found!');
else
  [~, id] = max([s.Area]);
  hold on, plot(s(id).Centroid(1),s(id).Centroid(2),'wp','MarkerSize',20,'MarkerFaceColor','r'), hold off
  disp(['Center location is (',num2str(s(id).Centroid(1),4),', ',num2str(s(id).Centroid(2),4),')'])
end

%% Step 5: Verify estimated location
% Finally we will plot the center on the original image to clearly evaluate
% how well we have found the center.
imagesc(greenBall1);
hold on, plot(s(id).Centroid(1),s(id).Centroid(2),'wp','MarkerSize',20,'MarkerFaceColor','r'), hold off


text:
filename='Gripping_quad_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'PROJECTED GRADIENT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.5;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ n_data, a, x, fx ] = gamma_inc_tricomi_values ( n_data )

%*****************************************************************************80
%
%% GAMMA_INC_TRICOMI_VALUES: values of Tricomi's incomplete Gamma function.
%
%  Discussion:
%
%    Tricomi's incomplete Gamma function is defined as:
%
%      1/Gamma(A) * 1/X^A * Integral ( 0 <= T <= X ) T^(A-1) * exp(-T) dT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 April 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real A, the parameter of the function.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  a_vec = [ ...
     0.10E+00, ...
     0.10E+00, ...
     0.10E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.50E+00, ...
     0.10E+01, ...
     0.10E+01, ...
     0.10E+01, ...
     0.11E+01, ...
     0.11E+01, ...
     0.11E+01, ...
     0.20E+01, ...
     0.20E+01, ...
     0.20E+01, ...
     0.60E+01, ...
     0.60E+01, ...
     0.11E+02, ...
     0.26E+02, ...
     0.41E+02  ];

  fx_vec = [ ...
    1.048292641463504E+00, ...
    1.024577737369574E+00, ...
    0.9493712443185374E+00, ...
    1.100793230316492E+00, ...
    0.8998911979655218E+00, ...
    0.5301656062431039E+00, ...
    0.9516258196404043E+00, ...
    0.6321205588285577E+00, ...
    0.1986524106001829E+00, ...
    0.9071784510537487E+00, ...
    0.5891809618706485E+00, ...
    0.1688269752193589E+00, ...
    0.4527034271637121E+00, ...
    0.1965220442795224E+00, ...
    0.02025928457705232E+00, ...
    0.0001721181724479739E+00, ...
    3.280858070850586E-07, ...
    5.244396471821590E-14, ...
    2.013462926183376E-37, ...
    1.230623887499875E-68 ];

  x_vec = [ ...
     0.30E-01, ...
     0.30E+00, ...
     0.15E+01, ...
     0.75E-01, ...
     0.75E+00, ...
     0.35E+01, ...
     0.10E+00, ...
     0.10E+01, ...
     0.50E+01, ...
     0.10E+00, ...
     0.10E+01, ...
     0.50E+01, ...
     0.15E+00, ...
     0.15E+01, ...
     0.70E+01, ...
     0.25E+01, ...
     0.12E+02, ...
     0.16E+02, ...
     0.25E+02, ...
     0.45E+02 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    a = a_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function out = NL_TSTL_ReturnTime(y,NNR,maxT,past,Nref,embedParams)
% NL_TSTL_ReturnTime    Analysis of the histogram of return times.
%
% Return times are the time taken for the time series to return to a similar
% location in phase space for a given reference point.
%
% Strong peaks in the histogram are indicative of periodicities in the data.
%
%---INPUTS:
%
% y, scalar time series as a column vector
% NNR, number of nearest neighbours
% maxT, maximum return time to consider
% past, Theiler window
% Nref, number of reference indicies
% embedParams, to feed into BF_Embed
%
%---OUTPUTS: include basic measures from the histogram, including the occurrence of
% peaks, spread, proportion of zeros, and the distributional entropy.
%
% Uses the code, return_time, from TSTOOL.
% TSTOOL: http://www.physik3.gwdg.de/tstool/

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check Inputs
% ------------------------------------------------------------------------------
N = length(y); % length of the input time series

% Number of nearest neighbours, NNR
if nargin < 2 || isempty(NNR)
    NNR = 5;
end
if (NNR > 0) && (NNR < 1) % specify a proportion of time series length
    NNR = floor(NNR*N); if NNR == 0, NNR = 1; end
end

% Maximum return time, maxT
if nargin < 3 || isempty(maxT)
    maxT = 0.1;
end
if (maxT > 0) && (maxT <= 1) % specify a proportion
    maxT = floor(N*maxT);
    if maxT == 0, maxT = 1; end
end

% Theiler window, past
if nargin < 4 || isempty(past)
    past = 10;
end
if (past > 0) && (past < 1) % specify a proportion
    past = floor(N*past);
    if past == 0, past = 1; end % round up from 0
end

% Number of reference points
if nargin < 5 || isempty(Nref)
    Nref = -1; % use all available points
end

% embed parameters
if nargin < 6 || isempty(embedParams)
    embedParams = {'ac','fnnmar'};
    fprintf(1,'Using default embedding using autocorrelation and cao\n');
end

doPlot = false; % plot outputs to figures

% ------------------------------------------------------------------------------
%% Embed the signal
% ------------------------------------------------------------------------------
s = BF_Embed(y,embedParams{1},embedParams{2},1,true);
if ~isa(s,'signal') && isnan(s); % embedding failed
    warning('Embedding failed');
    out = NaN; return
end
numPoints = size(data(s),1);
if numPoints < 10
    % Set heuristic minimum (10) on the number of points needed to perform a meaningful analysis
    warning('Time series not long enough for return time analysis')
    out = NaN; return
end

% ------------------------------------------------------------------------------
%% Run the code
% ------------------------------------------------------------------------------
try
    rs = return_time(s, NNR, maxT, past, Nref);
catch emsg
    if strcmp(emsg.message,'Index exceeds matrix dimensions.')
        fprintf(1,'Error evaluating return_time\n');
        out = NaN; return
    else
        error(emsg.message);
    end
end

Trett = data(rs);
NN = length(Trett);

% ------------------------------------------------------------------------------
%% Quantify structure in output
% ------------------------------------------------------------------------------
out.max = max(Trett);
out.std = std(Trett);
out.pzeros = sum(Trett == 0)/NN;
out.pg05 = sum(Trett>max(Trett)*0.5)/NN;
out.iqr = iqr(Trett);

% recurrent peaks:
icross05 = find((Trett(1:end-1)-0.5*max(Trett)).*(Trett(2:end)-0.5*max(Trett)) < 0);
if ~isempty(icross05) && length(icross05) > 2
    difficross05 = diff(icross05);
    difficross05 = difficross05(difficross05 > 0.4*max(difficross05)); % remove small entries, crossing peaks

    out.meanpeaksep = mean(difficross05)/NN;
    out.maxpeaksep = max(difficross05)/NN;
    out.minpeaksep = min(difficross05)/NN;
    out.rangepeaksep = range(difficross05)/NN;
    out.stdpeaksep = std(difficross05)/sqrt(NN);
else
    out.meanpeaksep = NaN;
    out.maxpeaksep = NaN;
    out.minpeaksep = NaN;
    out.rangepeaksep = NaN;
    out.stdpeaksep = NaN;
end

out.statrtys = std(Trett(1:floor(end/2)))/std(Trett(floor(end/2)+1:end));
out.statrtym = mean(Trett(1:floor(end/2)))/mean(Trett(floor(end/2)+1:end));

out.hhist = -sum(Trett(Trett>0).*log(Trett(Trett>0)));

% ------------------------------------------------------------------------------
%% Coarse-grain to 20 bins
% ------------------------------------------------------------------------------
numBins = 20;
cglav = zeros(numBins,1);
inds = round(linspace(0,NN,numBins+1));
for i = 1:numBins
    cglav(i) = sum(Trett(inds(i)+1:inds(i+1)));
end
if doPlot
    figure('color','w');
    box('on');
    plot(cglav,'k')
end
out.hcgdist = -sum(cglav(cglav > 0).*log(cglav(cglav > 0)));
out.rangecgdist = range(cglav);
out.pzeroscgdist = sum(cglav == 0)/numBins;

% ------------------------------------------------------------------------------
%% Get distribution of distribution of return times
% ------------------------------------------------------------------------------
[nhist, binEdges] = histcounts(Trett,'BinMethod','sqrt','Normalization','probability');
if doPlot
    binCenters = mean([binEdges(1:end-1); binEdges(2:end)]);
    figure('color','w');
    plot(binCenters,nhist,'o-k')
end
out.maxhisthist = max(nhist);
out.phisthistmin = nhist(1); % this is the same as maxhisthist
out.hhisthist = -sum(nhist(nhist > 0).*log(nhist(nhist > 0)));

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function so = tail(si)
% TAIL As streamCdr
so = streamCdr(si);

author: Sable

repo: mcbench-benchmarks


text:
function poly2 = resamplePolygon(poly, n)
%RESAMPLEPOLYGON  Distribute N points equally spaced on a polygon.
%
%   POLY2 = resamplePolygon(POLY, N)
%   Resample the input polygon POLY such that the resulting polygon POLY2
%   has N vertices. All points of POLY2 belong to the initial polygon, but
%   are not necessarily vertices of the original polygon.
%
%
%   Example
%     % creates a polygon from an ellipse
%     elli = [20 30 40 20 30];
%     poly = ellipseToPolygon(elli, 500);
%     figure; drawPolygon(poly, 'b');
%     % resample the polygon with a fixed number of vertices
%     poly2 = resamplePolygon(poly, 20);
%     drawPolygon(poly2, 'm');
%     drawPoint(poly2, 'mo');
%     axis equal; axis([-20 60 0 60]);
%
%   See also
%     polygons2d, resamplePolygonByLength, smoothPolygon, resamplePolyline
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2011-12-09, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

poly2 = resamplePolyline(poly([1:end 1],:), n+1);
poly2(end, :) = [];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function  Raw = simObservation(SimRob, SimSen, SimLmk, Opt)

% SIMOBSERVATION Observe simulated landmarks.
%   RAW = SIMOBSERVATION(SIMROB,SIMSEN,SIMLMK) returns the raw data
%   captured for the sensor.
%       SIMROB: simulated robot structure
%       SIMSEN: simulated sensor strucure
%       SIMLMK: simulated landmarks strucure
%
%   See also SIMMOTION PROJEUCPNTINTOPINHOLEONROB

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

Raw.type = 'simu';

switch SimSen.type

    case {'pinHole'}      % camera pinHole

        % Project virtual world's points
        [Raw.data.points.coord, s] = projEucPntIntoPinHoleOnRob(...
            SimRob.frame, ...
            SimSen.frame, ...
            SimSen.par.k, ...
            SimSen.par.d, ...
            SimLmk.points.coord);
        Raw.data.points.app  = SimLmk.points.id;

        % Add sensor noise
        Raw.data.points.coord = Raw.data.points.coord + ...
            SimSen.par.pixErr*randn(size(Raw.data.points.coord));

        % Remove non visible
        vis = isVisible(Raw.data.points.coord,s,SimSen.par.imSize);
        Raw.data.points.coord(:, ~vis)  = [];
        Raw.data.points.app(~vis) = [];

        % Project virtual world's segments
        [Raw.data.segments.coord, s] = projSegLinIntoPinHoleOnRob(...
            SimRob.frame, ...
            SimSen.frame, ...
            SimSen.par.k, ...
            SimLmk.segments.coord);
        Raw.data.segments.app  = SimLmk.segments.id;

        % Add sensor noise
        Raw.data.segments.coord = Raw.data.segments.coord + ...
            SimSen.par.cov*randn(size(Raw.data.points.coord));

        % Remove non visible
        [Raw.data.segments.coord,vis] = visibleSegment( ...
            Raw.data.segments.coord,...
            s,...
            SimSen.par.imSize,...
            0,...                      % N pixels margin
            Opt.obs.lines.minLength);  % min segment length
        Raw.data.segments.coord(:, ~vis)  = [];
        Raw.data.segments.app(~vis) = [];

    case {'pinHoleDepth'}      % camera pinHole

        % Project virtual world's points
        Raw.data.points.coord = projEucPntIntoPhdOnRob(...
            SimRob.frame, ...
            SimSen.frame, ...
            SimSen.par.k, ...
            SimSen.par.d, ...
            SimLmk.points.coord);
        Raw.data.points.app  = SimLmk.points.id;

        % Add sensor noise
        Raw.data.points.coord = Raw.data.points.coord + ...
            SimSen.par.cov*randn(size(Raw.data.points.coord));

        % Remove non visible
        vis = isVisible(Raw.data.points.coord(1:2,:),Raw.data.points.coord(3,:),SimSen.par.imSize);
        Raw.data.points.coord(:, ~vis)  = [];
        Raw.data.points.app(~vis) = [];

        % TODO: Lines not implemented for sensor type 'pihHoleDepth'

        %         % Project virtual world's segments
        Raw.data.segments.coord = [];
        %         [Raw.data.segments.coord, s] = projSegLinIntoPinHoleOnRob(...
        %             SimRob.frame, ...
        %             SimSen.frame, ...
        %             SimSen.par.k, ...
        %             SimLmk.segments.coord);
        %         Raw.data.segments.app  = SimLmk.segments.id;
        Raw.data.segments.app = [];
        %
        %         % Add sensor noise
        %         Raw.data.segments.coord = Raw.data.segments.coord + ...
        %             SimSen.par.pixErr*randn(size(Raw.data.segments.coord));
        %
        %         % Remove non visible
        %         [Raw.data.segments.coord,vis] = visibleSegment( ...
        %             Raw.data.segments.coord,...
        %             s,...
        %             SimSen.par.imSize,...
        %             0,...                      % N pixels margin
        %             Opt.obs.lines.minLength);  % min segment length
        %         Raw.data.segments.coord(:, ~vis)  = [];
        %         Raw.data.segments.app(~vis) = [];

    case {'omniCam'}      % Omnidirectional camera

        % Project virtual world's points
        [Raw.data.points.coord, s] = projEucPntIntoOmniCamOnRob(...
            SimRob.frame, ...
            SimSen.frame, ...
            SimSen.par.k, ...
            SimSen.par.d, ...
            SimLmk.points.coord);
        Raw.data.points.app  = SimLmk.points.id;

        % Add sensor noise
        Raw.data.points.coord = Raw.data.points.coord + ...
            SimSen.par.pixErr*randn(size(Raw.data.points.coord));

        % Remove non visible
        vis = isVisible(Raw.data.points.coord,s,SimSen.par.imSize);
        Raw.data.points.coord(:, ~vis)  = [];
        Raw.data.points.app(~vis) = [];

    otherwise
        % Print an error and exit
        error('??? Unknown sensor type ''%s''.',SimSen.type);

end % end of the "switch" on sensor type

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [ w, xyz ] = tetrahedron_unit_o01 ( )

%*****************************************************************************80
%
%% TETRAHEDRON_UNIT_O01 returns a 1 point quadrature rule for the unit tetrahedron.
%
%  Discussion:
%
%    The integration region is:
%
%      0 <= X
%      0 <= Y
%      0 <= Z
%      X + Y + Z <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carlos Felippa,
%    A compendium of FEM integration formulas for symbolic work,
%    Engineering Computation,
%    Volume 21, Number 8, 2004, pages 867-890.
%
%  Parameters:
%
%    Output, real W(1), the weights.
%
%    Output, real XYZ(3,1), the abscissas.
%
  w(1:1,1) = [ ...
    1.0 ];

  xyz(1:3,1:1) = [ ...
    0.25000000000000000000,  0.25000000000000000000,  0.25000000000000000000 ]';

  return
end


text:
function arctan2_values_test ( )

%*****************************************************************************80
%
%% ARCTAN2_VALUES_TEST demonstrates the use of ARCTAN2_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ARCTAN2_VALUES_TEST:\n' );
  fprintf ( 1, '  ARCTAN2_VALUES stores values of \n' );
  fprintf ( 1, '  the arc tangent function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           Y            F\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, y, f ] = arctan2_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %12f  %24.16f\n', x, y, f );

  end

  return
end


text:
classdef set

methods(Static)
    function subset = random_subset_1(set_size, subset_size)
        subset = randperm(set_size, subset_size);
    end

    function subset = random_subset_2(set_size, subset_size)
        flags = zeros(1, set_size);
        if (subset_size == 0)
            subset = [];
        end
        n = set_size;
        m = subset_size;
        % we will iterate m times
        for i=n-m+1 : n
            % pick an entry between 1:i to be selected.
            % note that i is never yet selected.
            j = ceil(i * rand());
            if (flags(j) == 0)
                % if the entry is not already selected, pick it
                flags(j) =1 ;
            else
                % otherwise select the i-th entry
                flags(i) = 1;
            end
        end
        subset = find (flags == 1);
    end

end
end

author: indigits

repo: sparse-plex


text:
function speedbump_sim(v,t,ym)
%Animation function for the motion of a 1 DOF auto suspension as it travels
%over a sinusoidal speed bump
%Written by T. Nordenholz, Fall 05
%To use, type speedbump_sim(v,t,x) where v is the horizontal speed (m/s) of the car,
%t is the time(sec)array, and x is the vertical (upwards) displacement (m)
%Geometrical and plotting parameters can be set within this program

%set geometric parameters
tsd=10; %slow down factor(ie animation will be tsd x slower than real time)
w=.8;%width of speedbump
h=.1;%height of speed bump
massw=.5;%length of mass
massh=.1;%height of mass
L0=.6-massh/2;%unstretched length of spring
Wsd=.1;
x=v*t;%horizontal distance travelled by car
%generate the base profile yb
for n=1:length(x)
    if x(n)<=0|x(n)>=w
        yb(n)=0;
    else
        yb(n)=h*sin(pi/w*x(n));
    end
end
massx=massw*[-.5,-.5,.5,.5,-.5];%data for drawing mass
massy=massh*[-.5,.5,.5,-.5,-.5];
y=ym+L0+massh/2; %vertical position of mass center

%set up figure and initialze plot/animation
Hf=figure('Position',[1,100,1000,500]);
%ground
Hp_ground=area(x,yb);axis([x(1)-massw,x(end)+massw,0,2*L0]),grid on,xlabel('s(m)'),ylabel('y (m)')
hold on
set(Hp_ground,'FaceColor','r');
%plot of x vs t
Hl_mass=line(x(1),y(1));
set(Hl_mass,'LineWidth',2)
%mass
Hp_mass=fill(x(1)+massx,y(1)+massy,'b');
%center of mass marker
Hl_cm=line(x(1),y(1),'Marker','O','MarkerSize',6,'MarkerFaceColor','k');
titlestr=['v = ',num2str(v),' m/s'];
title(titlestr)

% spring/damper
Hgt_springdamp=hgtransform;
Hl_Bend=line([0,0],[0,.1],'Color','k','Parent',Hgt_springdamp);
Hl_Tend=line([0,0],[.9,1],'Color','k','Parent',Hgt_springdamp);
Hl_Bbar=line(Wsd*[-1,1],[.1,.1],'Color','k','Parent',Hgt_springdamp);
Hl_Tbar=line(Wsd*[-1,1],[.9,.9],'Color','k','Parent',Hgt_springdamp);
Hl_spring=line(Wsd*[1,2,1,0,1,2,1,0,1],linspace(.1,.9,9),'Color','k','Parent',Hgt_springdamp);
Hl_dampB=line(Wsd*[-1,-1],[.1,.4],'Color','k','Parent',Hgt_springdamp);
Hl_dampBpist=line(Wsd*[-1.3,-.7],[.4,.4],'Color','k','Parent',Hgt_springdamp);
Hl_dampT=line(Wsd*[-1,-1],[.6,.9],'Color','k','Parent',Hgt_springdamp);
Hl_dampTcyl=line(Wsd*[-.5,-.5,-1.5,-1.5],[.55,.6,.6,.55],'Color','k','Parent',Hgt_springdamp);
set(Hgt_springdamp,'Matrix',[1,0,0,x(1);0,L0,0,yb(1);0,0,1,0;0,0,0,1]);

%draw all above objects and hold for 1 second
drawnow
tic,while toc<1,end
tic

% run animation
for n=1:length(t)
    L=L0+ym(n)-yb(n); %set spring length
    set(Hl_mass,'XData',x(1:n),'YData',y(1:n));
    set(Hp_mass,'XData',x(n)+massx,'YData',y(n)+massy);
    set(Hl_cm,'XData',x(n),'YData',y(n));
    set(Hgt_springdamp,'Matrix',[1,0,0,x(n);0,L,0,yb(n);0,0,1,0;0,0,0,1]);
    % time delay(run tsd x slower than real time)
    while toc<tsd*(t(n)-t(1)),end
    time(n)=toc;
    drawnow;
end


text:
function [dStates, contactForces] = dynamics_singleStanceTwo(States, Actuators, Parameters)
% function [dStates, Actuators] = dynamics_singleStanceTwo(States, Actuators, Parameters)
%
% Computer Generated File -- DO NOT EDIT
%
% This function was created by the function Write_ContinuousDynamics()
% 10-Dec-2013 19:40:47
%
% Dymanics Model: retractable double pendulum biped
% Motion Phase: Single Stance Two
%
% Matthew Kelly
% Cornell University
%

x = States(1,:); % (m) Foot One horizontal position
y = States(2,:); % (m) Foot One vertical position
th1 = States(3,:); % (rad) Leg One absolute angle
th2 = States(4,:); % (rad) Leg Two absolute angle
L1 = States(5,:); % (m) Leg One length
L2 = States(6,:); % (m) Leg Two length
dx = States(7,:); % (m/s) Foot One horizontal velocity
dy = States(8,:); % (m/s) Foot One vertical velocity
dth1 = States(9,:); % (rad/s) Leg One absolute angular rate
dth2 = States(10,:); % (rad/s) Leg Two absolute angular rate
dL1 = States(11,:); % (m/s) Leg One extension rate
dL2 = States(12,:); % (m/s) Leg Two extensioin rate

F1 = Actuators(1,:); % (N) Compresive axial force in Leg One
F2 = Actuators(2,:); % (N) Compresive axial force in Leg Two
T1 = 0; %Foot One not in contact with ground!
T2 = Actuators(4,:); % (Nm) External torque applied to Leg Two
Thip = Actuators(5,:); % (Nm) Torque acting on Leg Two from Leg One

m1 = Parameters.m1; % (kg) Foot One mass
m2 = Parameters.m2; % (kg) Foot Two mass
M = Parameters.M; % (kg) Hip mass
g = Parameters.g; % (m/s^2) Gravity

% Constraints for this phase:
H1 = 0; %(N) Foot One, horizontal contact force
V1 = 0; %(N) Foot One, vertical contact force

dStates = zeros(size(States));
dStates(1:6,:) = States((1+6):(6+6),:);
dStates(7,:) = (T1.*sin(th1) - Thip.*sin(th1) + H1.*L1 - F1.*L1.*cos(th1))./(L1.*m1);
dStates(8,:) = -(T1.*cos(th1) - L1.*V1 - Thip.*cos(th1) + L1.*g.*m1 + F1.*L1.*sin(th1))./(L1.*m1);
dStates(9,:) = -(L2.*Thip.*m1 - L2.*T1.*m1 + L1.*T2.*m1.*cos(th1 - th2) + L1.*Thip.*m1.*cos(th1 - th2) - L2.*M.*T1.*cos(th1).^2 + L2.*M.*Thip.*cos(th1).^2 - L2.*M.*T1.*sin(th1).^2 + L2.*M.*Thip.*sin(th1).^2 + L1.*L2.*M.*V1.*cos(th1) - H1.*L1.*L2.*M.*sin(th1) - F2.*L1.*L2.*m1.*sin(th1 - th2) + 2.*L1.*L2.*M.*dL1.*dth1.*m1)./(L1.^2.*L2.*M.*m1);
dStates(10,:) = (L2.*M.*T1.*sin(th1).*sin(th2) - L2.*T1.*m1.*cos(th1).*cos(th2) - L2.*M.*Thip.*sin(th1).*sin(th2) + L2.*Thip.*m1.*cos(th1).*cos(th2) - L2.*T1.*m1.*sin(th1).*sin(th2) + L2.*Thip.*m1.*sin(th1).*sin(th2) - L2.*M.*T1.*cos(th1).^3.*cos(th2) + L2.*M.*Thip.*cos(th1).^3.*cos(th2) - L1.*L2.*M.*V1.*cos(th2) + H1.*L1.*L2.*M.*sin(th2) - L2.*M.*T1.*sin(th1).^3.*sin(th2) + L2.*M.*Thip.*sin(th1).^3.*sin(th2) + L2.*M.*T1.*cos(th1).*cos(th2) - L2.*M.*Thip.*cos(th1).*cos(th2) + L1.*T2.*m1.*cos(th1 - th2).*cos(th1).*cos(th2) + L1.*Thip.*m1.*cos(th1 - th2).*cos(th1).*cos(th2) - L2.*M.*T1.*cos(th1).*cos(th2).*sin(th1).^2 + L2.*M.*Thip.*cos(th1).*cos(th2).*sin(th1).^2 + L1.*T2.*m1.*cos(th1 - th2).*sin(th1).*sin(th2) - L1.*T2.*m1.*sin(th1 - th2).*cos(th1).*sin(th2) + L1.*T2.*m1.*sin(th1 - th2).*cos(th2).*sin(th1) + L1.*Thip.*m1.*cos(th1 - th2).*sin(th1).*sin(th2) - L1.*Thip.*m1.*sin(th1 - th2).*cos(th1).*sin(th2) + L1.*Thip.*m1.*sin(th1 - th2).*cos(th2).*sin(th1) - F1.*L1.*L2.*M.*cos(th1).*sin(th2) + F1.*L1.*L2.*M.*cos(th2).*sin(th1) - L2.*M.*T1.*cos(th1).^2.*sin(th1).*sin(th2) + L2.*M.*Thip.*cos(th1).^2.*sin(th1).*sin(th2) + F1.*L1.*L2.*m1.*cos(th1).*sin(th2) - F1.*L1.*L2.*m1.*cos(th2).*sin(th1) - F1.*L1.*L2.*M.*cos(th2).*sin(th1).^3 + F1.*L1.*L2.*M.*cos(th1).^3.*sin(th2) + L1.*L2.*M.*V1.*cos(th1).^2.*cos(th2) - H1.*L1.*L2.*M.*cos(th1).^2.*sin(th2) + L1.*L2.*M.*V1.*cos(th2).*sin(th1).^2 - H1.*L1.*L2.*M.*sin(th1).^2.*sin(th2) + L1.*L2.*M.*g.*m1.*cos(th2) - F2.*L1.*L2.*m1.*cos(th1 - th2).*cos(th1).*sin(th2) + F2.*L1.*L2.*m1.*cos(th1 - th2).*cos(th2).*sin(th1) - F2.*L1.*L2.*m1.*sin(th1 - th2).*cos(th1).*cos(th2) - F1.*L1.*L2.*M.*cos(th1).^2.*cos(th2).*sin(th1) + F1.*L1.*L2.*M.*cos(th1).*sin(th1).^2.*sin(th2) - F2.*L1.*L2.*m1.*sin(th1 - th2).*sin(th1).*sin(th2) - 2.*L1.*L2.*M.*dL2.*dth2.*m1.*cos(th2).^2 - 2.*L1.*L2.*M.*dL2.*dth2.*m1.*sin(th2).^2)./(L1.*L2.^2.*M.*m1.*(cos(th2).^2 + sin(th2).^2));
dStates(11,:) = -(T2.*m1.*sin(th1 - th2) + Thip.*m1.*sin(th1 - th2) - F1.*L2.*m1 + H1.*L2.*M.*cos(th1) + L2.*M.*V1.*sin(th1) + F2.*L2.*m1.*cos(th1 - th2) - F1.*L2.*M.*cos(th1).^2 - F1.*L2.*M.*sin(th1).^2 - L1.*L2.*M.*dth1.^2.*m1)./(L2.*M.*m1);
dStates(12,:) = (L2.*M.*T1.*cos(th1).*sin(th2) - L2.*M.*T1.*cos(th2).*sin(th1) - L2.*M.*Thip.*cos(th1).*sin(th2) + L2.*M.*Thip.*cos(th2).*sin(th1) - L2.*T1.*m1.*cos(th1).*sin(th2) + L2.*T1.*m1.*cos(th2).*sin(th1) + L2.*Thip.*m1.*cos(th1).*sin(th2) - L2.*Thip.*m1.*cos(th2).*sin(th1) - H1.*L1.*L2.*M.*cos(th2) + L2.*M.*T1.*cos(th2).*sin(th1).^3 - L2.*M.*T1.*cos(th1).^3.*sin(th2) - L2.*M.*Thip.*cos(th2).*sin(th1).^3 + L2.*M.*Thip.*cos(th1).^3.*sin(th2) - L1.*L2.*M.*V1.*sin(th2) + L1.*L2.*M.*g.*m1.*sin(th2) + L1.*L2.^2.*M.*dth2.^2.*m1.*cos(th2).^2 + L1.*L2.^2.*M.*dth2.^2.*m1.*sin(th2).^2 + L1.*T2.*m1.*cos(th1 - th2).*cos(th1).*sin(th2) - L1.*T2.*m1.*cos(th1 - th2).*cos(th2).*sin(th1) + L1.*T2.*m1.*sin(th1 - th2).*cos(th1).*cos(th2) + L1.*Thip.*m1.*cos(th1 - th2).*cos(th1).*sin(th2) - L1.*Thip.*m1.*cos(th1 - th2).*cos(th2).*sin(th1) + L1.*Thip.*m1.*sin(th1 - th2).*cos(th1).*cos(th2) + F1.*L1.*L2.*M.*cos(th1).*cos(th2) + L2.*M.*T1.*cos(th1).^2.*cos(th2).*sin(th1) - L2.*M.*Thip.*cos(th1).^2.*cos(th2).*sin(th1) - L2.*M.*T1.*cos(th1).*sin(th1).^2.*sin(th2) + L2.*M.*Thip.*cos(th1).*sin(th1).^2.*sin(th2) + L1.*T2.*m1.*sin(th1 - th2).*sin(th1).*sin(th2) + L1.*Thip.*m1.*sin(th1 - th2).*sin(th1).*sin(th2) + F1.*L1.*L2.*M.*sin(th1).*sin(th2) - F1.*L1.*L2.*m1.*cos(th1).*cos(th2) - F1.*L1.*L2.*m1.*sin(th1).*sin(th2) - F1.*L1.*L2.*M.*cos(th1).^3.*cos(th2) + H1.*L1.*L2.*M.*cos(th1).^2.*cos(th2) + H1.*L1.*L2.*M.*cos(th2).*sin(th1).^2 - F1.*L1.*L2.*M.*sin(th1).^3.*sin(th2) + L1.*L2.*M.*V1.*cos(th1).^2.*sin(th2) + L1.*L2.*M.*V1.*sin(th1).^2.*sin(th2) + F2.*L1.*L2.*m1.*cos(th1 - th2).*cos(th1).*cos(th2) - F1.*L1.*L2.*M.*cos(th1).*cos(th2).*sin(th1).^2 + F2.*L1.*L2.*m1.*cos(th1 - th2).*sin(th1).*sin(th2) - F2.*L1.*L2.*m1.*sin(th1 - th2).*cos(th1).*sin(th2) + F2.*L1.*L2.*m1.*sin(th1 - th2).*cos(th2).*sin(th1) - F1.*L1.*L2.*M.*cos(th1).^2.*sin(th1).*sin(th2))./(L1.*L2.*M.*m1.*(cos(th2).^2 + sin(th2).^2));

% contactForces(1,:) == H1 == (N) Foot One, horizontal contact force
% contactForces(2,:) == V1 == (N) Foot One, vertical contact force
% contactForces(3,:) == H2 == (N) Foot Two, horizontal contact force
% contactForces(4,:) == V2 == (N) Foot Two, vertical contact force
contactForces = zeros(4,size(States,2));
contactForces(1,:) = zeros(1,size(States,2));
contactForces(2,:) = zeros(1,size(States,2));
contactForces(3,:) = (L2.*M.*T1.*m2.*sin(th1).^3 - L2.*M.*Thip.*m2.*sin(th1).^3 + L1.*M.*T2.*m1.*sin(th2) - L2.*M.*T1.*m2.*sin(th1) + L1.*M.*Thip.*m1.*sin(th2) + L2.*M.*Thip.*m2.*sin(th1) + L2.*T1.*m1.*m2.*sin(th1) - L2.*Thip.*m1.*m2.*sin(th1) - H1.*L1.*L2.*M.*m2 + L1.*T2.*m1.*m2.*cos(th1 - th2).^2.*sin(th2) + L1.*Thip.*m1.*m2.*cos(th1 - th2).^2.*sin(th2) - F1.*L1.*L2.*M.*m2.*cos(th1).^3 + H1.*L1.*L2.*M.*m2.*cos(th1).^2 + H1.*L1.*L2.*M.*m2.*cos(th2).^2 + L1.*T2.*m1.*m2.*sin(th1 - th2).^2.*sin(th2) + L1.*Thip.*m1.*m2.*sin(th1 - th2).^2.*sin(th2) + H1.*L1.*L2.*M.*m2.*sin(th1).^2 + H1.*L1.*L2.*M.*m2.*sin(th2).^2 - L1.*T2.*m1.*m2.*cos(th1 - th2).*sin(th1) + L1.*T2.*m1.*m2.*sin(th1 - th2).*cos(th1) - L2.*T1.*m1.*m2.*cos(th1 - th2).*sin(th2) - L2.*T1.*m1.*m2.*sin(th1 - th2).*cos(th2) - L1.*Thip.*m1.*m2.*cos(th1 - th2).*sin(th1) + L1.*Thip.*m1.*m2.*sin(th1 - th2).*cos(th1) + L2.*Thip.*m1.*m2.*cos(th1 - th2).*sin(th2) + L2.*Thip.*m1.*m2.*sin(th1 - th2).*cos(th2) + F1.*L1.*L2.*M.*m2.*cos(th1) - F2.*L1.*L2.*M.*m1.*cos(th2) + L2.*M.*T1.*m2.*cos(th1).^2.*sin(th1) + L2.*M.*T1.*m2.*cos(th2).^2.*sin(th1) - L2.*M.*Thip.*m2.*cos(th1).^2.*sin(th1) - L2.*M.*Thip.*m2.*cos(th2).^2.*sin(th1) + L2.*M.*T1.*m2.*sin(th1).*sin(th2).^2 - L2.*M.*Thip.*m2.*sin(th1).*sin(th2).^2 - F1.*L1.*L2.*m1.*m2.*cos(th1) - F1.*L1.*L2.*M.*m2.*cos(th1).*cos(th2).^2 - F1.*L1.*L2.*M.*m2.*cos(th1).*sin(th1).^2 - F1.*L1.*L2.*M.*m2.*cos(th1).*sin(th2).^2 - F1.*L1.*L2.*m1.*m2.*sin(th1 - th2).*sin(th2) + F2.*L1.*L2.*m1.*m2.*sin(th1 - th2).*sin(th1) - F2.*L1.*L2.*m1.*m2.*cos(th1 - th2).^2.*cos(th2) - F2.*L1.*L2.*m1.*m2.*sin(th1 - th2).^2.*cos(th2) - L2.*M.*T1.*m2.*cos(th1 - th2).*cos(th1).^2.*sin(th2) - L2.*M.*T1.*m2.*sin(th1 - th2).*cos(th1).^2.*cos(th2) + L2.*M.*Thip.*m2.*cos(th1 - th2).*cos(th1).^2.*sin(th2) + L2.*M.*Thip.*m2.*sin(th1 - th2).*cos(th1).^2.*cos(th2) - L2.*M.*T1.*m2.*cos(th1 - th2).*sin(th1).^2.*sin(th2) - L2.*M.*T1.*m2.*sin(th1 - th2).*cos(th2).*sin(th1).^2 + L2.*M.*Thip.*m2.*cos(th1 - th2).*sin(th1).^2.*sin(th2) + L2.*M.*Thip.*m2.*sin(th1 - th2).*cos(th2).*sin(th1).^2 + F1.*L1.*L2.*m1.*m2.*cos(th1 - th2).*cos(th2) + F2.*L1.*L2.*m1.*m2.*cos(th1 - th2).*cos(th1) + F1.*L1.*L2.*M.*m2.*cos(th1 - th2).*cos(th1).^2.*cos(th2) + F1.*L1.*L2.*M.*m2.*cos(th1 - th2).*cos(th2).*sin(th1).^2 - F1.*L1.*L2.*M.*m2.*sin(th1 - th2).*cos(th1).^2.*sin(th2) - F1.*L1.*L2.*M.*m2.*sin(th1 - th2).*sin(th1).^2.*sin(th2) - H1.*L1.*L2.*M.*m2.*cos(th1 - th2).*cos(th1).*cos(th2) + L1.*L2.*M.*V1.*m2.*cos(th1 - th2).*cos(th1).*sin(th2) - L1.*L2.*M.*V1.*m2.*cos(th1 - th2).*cos(th2).*sin(th1) + L1.*L2.*M.*V1.*m2.*sin(th1 - th2).*cos(th1).*cos(th2) - H1.*L1.*L2.*M.*m2.*cos(th1 - th2).*sin(th1).*sin(th2) + H1.*L1.*L2.*M.*m2.*sin(th1 - th2).*cos(th1).*sin(th2) - H1.*L1.*L2.*M.*m2.*sin(th1 - th2).*cos(th2).*sin(th1) + L1.*L2.*M.*V1.*m2.*sin(th1 - th2).*sin(th1).*sin(th2))./(L1.*L2.*M.*m1.*(cos(th2).^2 + sin(th2).^2));
contactForces(4,:) = -(L1.*L2.*M.*V1.*m2 + L1.*M.*T2.*m1.*cos(th2) - L2.*M.*T1.*m2.*cos(th1) + L1.*M.*Thip.*m1.*cos(th2) + L2.*M.*Thip.*m2.*cos(th1) + L2.*T1.*m1.*m2.*cos(th1) - L2.*Thip.*m1.*m2.*cos(th1) + L2.*M.*T1.*m2.*cos(th1).^3 - L2.*M.*Thip.*m2.*cos(th1).^3 + L1.*T2.*m1.*m2.*cos(th1 - th2).^2.*cos(th2) + L1.*Thip.*m1.*m2.*cos(th1 - th2).^2.*cos(th2) + L1.*T2.*m1.*m2.*sin(th1 - th2).^2.*cos(th2) + L1.*Thip.*m1.*m2.*sin(th1 - th2).^2.*cos(th2) + F1.*L1.*L2.*M.*m2.*sin(th1).^3 - L1.*L2.*M.*V1.*m2.*cos(th1).^2 - L1.*L2.*M.*V1.*m2.*cos(th2).^2 - L1.*L2.*M.*V1.*m2.*sin(th1).^2 - L1.*L2.*M.*V1.*m2.*sin(th2).^2 - L1.*L2.*M.*g.*m1.*m2 - L1.*T2.*m1.*m2.*cos(th1 - th2).*cos(th1) - L2.*T1.*m1.*m2.*cos(th1 - th2).*cos(th2) - L1.*Thip.*m1.*m2.*cos(th1 - th2).*cos(th1) + L2.*Thip.*m1.*m2.*cos(th1 - th2).*cos(th2) + L2.*M.*T1.*m2.*cos(th1).*cos(th2).^2 - L2.*M.*Thip.*m2.*cos(th1).*cos(th2).^2 + L2.*M.*T1.*m2.*cos(th1).*sin(th1).^2 + L2.*M.*T1.*m2.*cos(th1).*sin(th2).^2 - L2.*M.*Thip.*m2.*cos(th1).*sin(th1).^2 - L2.*M.*Thip.*m2.*cos(th1).*sin(th2).^2 - L1.*T2.*m1.*m2.*sin(th1 - th2).*sin(th1) + L2.*T1.*m1.*m2.*sin(th1 - th2).*sin(th2) - L1.*Thip.*m1.*m2.*sin(th1 - th2).*sin(th1) - L2.*Thip.*m1.*m2.*sin(th1 - th2).*sin(th2) - F1.*L1.*L2.*M.*m2.*sin(th1) + F2.*L1.*L2.*M.*m1.*sin(th2) + F1.*L1.*L2.*m1.*m2.*sin(th1) + L2.*M.*T1.*m2.*sin(th1 - th2).*sin(th1).^2.*sin(th2) - L2.*M.*Thip.*m2.*sin(th1 - th2).*sin(th1).^2.*sin(th2) - F1.*L1.*L2.*m1.*m2.*cos(th1 - th2).*sin(th2) - F1.*L1.*L2.*m1.*m2.*sin(th1 - th2).*cos(th2) - F2.*L1.*L2.*m1.*m2.*cos(th1 - th2).*sin(th1) + F2.*L1.*L2.*m1.*m2.*sin(th1 - th2).*cos(th1) + F1.*L1.*L2.*M.*m2.*cos(th1).^2.*sin(th1) + F1.*L1.*L2.*M.*m2.*cos(th2).^2.*sin(th1) + F1.*L1.*L2.*M.*m2.*sin(th1).*sin(th2).^2 + F2.*L1.*L2.*m1.*m2.*cos(th1 - th2).^2.*sin(th2) + F2.*L1.*L2.*m1.*m2.*sin(th1 - th2).^2.*sin(th2) - L2.*M.*T1.*m2.*cos(th1 - th2).*cos(th1).^2.*cos(th2) + L2.*M.*Thip.*m2.*cos(th1 - th2).*cos(th1).^2.*cos(th2) - L2.*M.*T1.*m2.*cos(th1 - th2).*cos(th2).*sin(th1).^2 + L2.*M.*Thip.*m2.*cos(th1 - th2).*cos(th2).*sin(th1).^2 + L2.*M.*T1.*m2.*sin(th1 - th2).*cos(th1).^2.*sin(th2) - L2.*M.*Thip.*m2.*sin(th1 - th2).*cos(th1).^2.*sin(th2) - F1.*L1.*L2.*M.*m2.*cos(th1 - th2).*cos(th1).^2.*sin(th2) - F1.*L1.*L2.*M.*m2.*sin(th1 - th2).*cos(th1).^2.*cos(th2) - F1.*L1.*L2.*M.*m2.*cos(th1 - th2).*sin(th1).^2.*sin(th2) - F1.*L1.*L2.*M.*m2.*sin(th1 - th2).*cos(th2).*sin(th1).^2 + L1.*L2.*M.*V1.*m2.*cos(th1 - th2).*cos(th1).*cos(th2) + H1.*L1.*L2.*M.*m2.*cos(th1 - th2).*cos(th1).*sin(th2) - H1.*L1.*L2.*M.*m2.*cos(th1 - th2).*cos(th2).*sin(th1) + H1.*L1.*L2.*M.*m2.*sin(th1 - th2).*cos(th1).*cos(th2) + L1.*L2.*M.*V1.*m2.*cos(th1 - th2).*sin(th1).*sin(th2) - L1.*L2.*M.*V1.*m2.*sin(th1 - th2).*cos(th1).*sin(th2) + L1.*L2.*M.*V1.*m2.*sin(th1 - th2).*cos(th2).*sin(th1) + H1.*L1.*L2.*M.*m2.*sin(th1 - th2).*sin(th1).*sin(th2))./(L1.*L2.*M.*m1.*(cos(th2).^2 + sin(th2).^2));

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function data = getImageBatch(imagePaths, varargin)
% GETIMAGEBATCH  Load and jitter a batch of images

opts.useGpu = false ;
opts.prefetch = false ;
opts.numThreads = 1 ;

opts.imageSize = [227, 227] ;
opts.cropSize = 227 / 256 ;
opts.keepAspect = true ;
opts.subtractAverage = [] ;

opts.jitterFlip = false ;
opts.jitterLocation = false ;
opts.jitterAspect = 1 ;
opts.jitterScale = 1 ;
opts.jitterBrightness = 0 ;
opts.jitterContrast = 0 ;
opts.jitterSaturation = 0 ;

opts = vl_argparse(opts, varargin);

args{1} = {imagePaths, ...
           'NumThreads', opts.numThreads, ...
           'Pack', ...
           'Interpolation', 'bicubic', ...
           'Resize', opts.imageSize(1:2), ...
           'CropSize', opts.cropSize * opts.jitterScale, ...
           'CropAnisotropy', opts.jitterAspect, ...
           'Brightness', opts.jitterBrightness, ...
           'Contrast', opts.jitterContrast, ...
           'Saturation', opts.jitterSaturation} ;

if ~opts.keepAspect
  % Squashign effect
  args{end+1} = {'CropAnisotropy', 0} ;
end

if opts.jitterFlip
  args{end+1} = {'Flip'} ;
end

if opts.jitterLocation
  args{end+1} = {'CropLocation', 'random'} ;
else
  args{end+1} = {'CropLocation', 'center'} ;
end

if opts.useGpu
  args{end+1} = {'Gpu'} ;
end

if ~isempty(opts.subtractAverage)
  args{end+1} = {'SubtractAverage', opts.subtractAverage} ;
end

args = horzcat(args{:}) ;

if opts.prefetch
  vl_imreadjpeg(args{:}, 'prefetch') ;
  data = [] ;
else
  data = vl_imreadjpeg(args{:}) ;
  data = data{1} ;
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function X_rec = recoverData(Z, U, K)
%RECOVERDATA Recovers an approximation of the original data when using the
%projected data
%   X_rec = RECOVERDATA(Z, U, K) recovers an approximation the
%   original data that has been reduced to K dimensions. It returns the
%   approximate reconstruction in X_rec.
%

% You need to return the following variables correctly.
X_rec = zeros(size(Z, 1), size(U, 1));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the approximation of the data by projecting back
%               onto the original space using the top K eigenvectors in U.
%
%               For the i-th example Z(i,:), the (approximate)
%               recovered data for dimension j is given as follows:
%                    v = Z(i, :)';
%                    recovered_j = v' * U(j, 1:K)';
%
%               Notice that U(j, 1:K) is a row vector.
%

for i = 1:size(Z,1)

	v = Z(i, :)';
	X_rec(i,:) = v' * U(:, 1:K)';

end

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

end

author: gopaczewski

repo: coursera-ml


text:
function [W,varW,X,varX] = gprotate2pca(W,varW,X,varX,weights)

[M,D] = size(W);
N = cols(X);

%WW = W'*W + diag(rowsum(varW));

% $$$ % Use ML zero mean (mu should be updated after this function!!) ..
% $$$ dmu = mean(X,2);

% $$$ % Move bias
% $$$ X = X - repmat(dmu,1,n);
% $$$ mu = mu + W*dmu;

Qx = 1;

%warning('DISCARDING VARIANCES!!')

% Whiten w.r.t. X
muX = mean(X,2);
X0 = bsxfun(@minus,X,muX);
XX = X0*X0' + diag(colsum(varX));
% $$$ if fixw
% $$$   [Vx,Dx,tmp] = svd(XX/(n-m)); % USE THIS IF FIXED w ??
% $$$ else
[Vx,Dx,tmp] = svd(XX/N);
% $$$ end
Qx = diag(1./sqrt(diag(Dx))) * Vx';
Qw = Vx*sqrt(Dx);
W = W * Qw;
for i=1:M
  varW(i,:) = diag(Qw'*diag(varW(i,:))*Qw);
end
X = Qx * X;
for j = 1:N
  varX(:,j) = diag(Qx*diag(varX(:,j))*Qx');
end

rotationX = Qx;

% $$$ % Check that XX is really whitened! (because of numerical issues!!)
% $$$ XX = X*X' + sum(CovX,3);
% $$$ if fixw
% $$$   [Vx,Dx,tmp] = svd(XX/(n-m)); % USE THIS IF FIXED w ??
% $$$ else
% $$$   [Vx,Dx,tmp] = svd(XX/n);
% $$$ end
% $$$ if Dx(1) > 1.1
% $$$   % Whiten w.r.t. X AGAIN
% $$$   Qx = diag(1./sqrt(diag(Dx))) * Vx';
% $$$   Qw = Vx*sqrt(Dx);
% $$$   W = W * Qw;
% $$$   for i=1:size(CovW,3)
% $$$     CovW(:,:,i) = Qw'*CovW(:,:,i)*Qw;
% $$$   end
% $$$   X = Qx * X;
% $$$   for j = 1:size(CovX,3)
% $$$     Sv{j} = Qx*Sv{j}*Qx';
% $$$     CovX(:,:,j) = Qx*CovX(:,:,j)*Qx';
% $$$   end
% $$$ end

Qx = 1;
% Diagonalize w.r.t. W
WW = W'*W + diag(rowsum(varW));
% Use weights!!
if nargin >= 5
  WW = W'*diag(weights(:))*W + diag(weights(:)' * varW);
%  WW = W'*diag(weights(:))*W + diag(rowsum(varW));
% $$$   w = sqrt(weights);
% $$$   WW = bsxfun(@times, WW, w(:));
% $$$   WW = bsxfun(@times, WW, w(:)');
end
%WW = W'*W + sum(CovW,3);
[Vw,Dw,tmp] = svd(WW);
%norms_of_W = sqrt(diag(Dw))
%[Dw,I] = sort(diag(Dw), 'descend');
%Vw = Vw(:,I);
Qx = Vw' * Qx;
Qw = Vw;
W = W * Qw;
for i=1:M
  varW(i,:) = diag(Qw'*diag(varW(i,:))*Qw);
end
X = Qx * X;
for j = 1:N
  varX(:,j) = diag(Qx*diag(varX(:,j))*Qx');
end

rotationX = Qx * rotationX;

% Rotate such that, the largest loadings are positive
for d=1:cols(W)
  if max(W(:,d)) < max(-W(:,d))
    W(:,d) = -W(:,d);
    X(d,:) = -X(d,:);
  end
end

% $$$ figure
% $$$ imagesc(abs(rotationX));
% $$$
% $$$ figure
% $$$ plot( sum(W.^2,1) )
% $$$
% $$$ figure
% $$$ plot( diag(Dw) )

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function i = end(x,k,n)
%END          Overloaded functions end, specifies last index
%

% written  10/03/02     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
%

  % call must be one-dimensional
  i = length(x.sign);

repo: OpenMAS


text:
function [nmps2] = ftps22nmps2(ftps2)
% Convert acceleration from feet per square-second to nanometers per second
% squared.
% Chad A. Greene 2012
nmps2 = ftps2*304800000;
end


text:
%Kalman Demo 1 (Sinc Regression)
%
%  The benefit from the state space formulation is that the
%  computational complexity is linear with respect to the number
%  of data points. Due to efficient matrix solvers in Matlab
%  (favoring the traditional GP solution over sequential looping),
%  the advantages in speed start to show in datasets with thousands
%  of data points.
%
%  The take-home message from this demo is that the GP can be
%  set up exactly as any other GP regression model in GPstuff,
%  and then solved by Kalman filtering methods by specifying the
%  'type' in the GP structure to be 'KALMAN'.
%
%  The implementation that is used below is primarily based on the
%  methods presented in the following publications:
%
%  [1] Simo Sarkka, Arno Solin, Jouni Hartikainen (2013).
%      Spatiotemporal learning via infinite-dimensional Bayesian
%      filtering and smoothing. IEEE Signal Processing Magazine,
%      30(4):51-61.
%
%  [2] Jouni Hartikainen and Simo Sarkka (2010). Kalman filtering and
%      smoothing solutions to temporal Gaussian process regression
%      models. Proceedings of IEEE International Workshop on Machine
%      Learning for Signal Processing (MLSP).
%
%  [3] Simo Sarkka (2013). Bayesian filtering and smoothing. Cambridge
%      University Press.
%
%  Copyright (c) 2014 Arno Solin and Jukka Koskenranta

%  This software is distributed under the GNU General Public
%  License (version 3 or later); please refer to the file
%  License.txt, included with the software, for details.

%% Generate data

  % Discretize x (measurement points)
  x = 10*rand(32,1)-5;

  % Test points (evaluation points)
  xt = linspace(-5,5,128)';

  % Measurement noise variance
  sigma2 = .01;

  % Simulate data
  y = sinc(x) + sqrt(sigma2)*randn(size(x));

%% Do the inference

  % Noise variance prior
  ps2 = prior_logunif();

  % The likelihood model
  lik = lik_gaussian('sigma2', 1, 'sigma2_prior', ps2);

  % Covariance function hyperparameter priors
  pl = prior_logunif();
  pm = prior_logunif();

  % The GP covariance function
  gpcf = gpcf_matern52('lengthScale', 1, ...
                       'magnSigma2', 1, ...
                       'lengthScale_prior', pl, ...
                       'magnSigma2_prior', pm);

  % Define Gaussian process model using type 'KALMAN'
  gp = gp_set('lik', lik, 'cf', gpcf, 'type', 'KALMAN');

  % Hyperparameter optimization (full model)
  gp_full = gp_optim(gp_set(gp,'type','FULL'), x, y);

  % Hyperparameter optimization (state space model)
  gp = gp_optim(gp, x, y);

  % Predict values at test inputs xt
  [Eft,Varft] = gp_pred(gp, x, y, 'xt', xt);

%% Compare against full GP solution (table)

  % Show table with comparison to full GP results
  fprintf('\n%12s | %8s | %8s \n', ...
      'Parameter','FULL', 'KALMAN')

  fprintf('-----------------------------------\n')

  fprintf('%12s | %8.4f | %8.4f \n', ...
      'magnSigma2',gp_full.cf{1}.magnSigma2,gp.cf{1}.magnSigma2)
  fprintf('%12s | %8.4f | %8.4f \n', ...
      'lengthScale',gp_full.cf{1}.lengthScale,gp.cf{1}.lengthScale)
  fprintf('%12s | %8.4f | %8.4f \n', ...
      'sigma2',gp_full.lik.sigma2,gp.lik.sigma2)

%% Show result

  % Plot the state space prediction
  figure(1); clf; hold on

    % Draw 95% uncertainty interval
    color=0.85*[1,1,1];
    p=patch([xt; flipud(xt)], ...
            [Eft+1.96*sqrt(Varft); flipud(Eft-1.96*sqrt(Varft))],color);
    set(p,'EdgeColor','none')

    % Show data and predicted mean
    h=plot(xt,sinc(xt),'-r', ...
           x,y,'+k', ...
           xt,Eft,'--','LineWidth', 1, 'MarkerSize',5);

    % Legend and labels
    xlabel('Input, x');
    ylabel('Output, y');
    title('Kalman Demo 1 (Sinc Regression)');
    legend([h;p], 'Sinc function', ...
                  'Noisy measurements', ...
                  'Estimated mean', ...
                  '95% confidence region');

    % Bring axes to front
    box on; set(gca,'Layer','top')

%% Compare the full posteriors

  % Return the posterior mean EFT and covariance COVFT of latent
  % variables:
  %     Eft =  E[f | xt,x,y,th]  = K_fy*(Kyy+s^2I)^(-1)*y
  %   Covft = Var[f | xt,x,y,th] = K_fy - K_fy*(Kyy+s^2I)^(-1)*K_yf.

  figure(2); clf

  % The full GP
  subplot(121)

    % Calculate
    [Eft,Covft] = gp_jpred(gp_full,x,y,xt);

    % Plot
    imagesc(Covft)
    axis equal tight

    % Title
    title('The full GP posterior covariance')

  % The state space model ('KALMAN')
  subplot(122)

    % Calculate
    [Eft,Covft] = gp_jpred(gp,x,y,xt);

    % Plot
    imagesc(Covft)
    axis equal tight

    % Title
    title('The state space GP posterior covariance')

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function y = mpower(x,d)
%MPOWER (overloaded)

if (numel(d)>1) || (size(x,1) ~= size(x,2))
   error('Inputs must be a scalar and a square matrix. To compute elementwise POWER, use POWER (.^) instead.');
end

if isa(d,'sdpvar')
    if numel(x) > 1
        if isa(x,'sdpvar')
            error('x^d not support SDPVARMATRIX^SDPVARSCALAR')
        else
            if isnumeric(x)
                [V,D] = eig(x);
                if ~isreal(D)
                    error('Matrix power x^d requires x to have real eigenvalues');
                end
                D = real(D);
                y = V*diag(diag(D).^d)*inv(V);
                y.extra.createTime = definecreationtime;
                return
            else
                error('Object class not support in x^d');
            end
        end
    end
    d = flush(d);d.conicinfo = [0 0];
    y = power_internal1(d,x);
    if isa(y,'sdpvar')
        y.extra.createTime = definecreationtime;
        y.extra.opname='';
    end
    return
end

x = flush(x);

% Trivial cases
if d==0
    y = eye(x.dim(1),x.dim(2))^0;
    return
end
if d==1
    y = x;
    return
end

% Check for special case norm(x)^2 which many users try to do
if d==2
    if length(x)==1
        base = getbase(x);
        if isequal(base,[0 1])
            if strcmp(x.extra.opname,'norm')
                model = yalmip('extstruct',getvariables(x));
                z = model.arg{1};
                if (isequal(model.arg{2},2) & min(size(z))==1) | isequal(model.arg{2},'fro')
                    z = reshape(model.arg{1},[],1);
                    y = real(z'*z);
                    y.extra.createTime = definecreationtime;
                     y.extra.opname='';
                    return
                end
            end
        end
    end
end

% Fractional and negative powers
if (ceil(d)-d>0) | (d<0)
    if x.dim(1)>1 | x.dim(2)>1
        error('Only scalars can have negative or non-integer powers');
    else
        base = getbase(x);
        if isequal(base,sparse([0 1])) % Simple unit scalar
            [mt,variabletype,hashM,hash] = yalmip('monomtable');
            var = getvariables(x);
          %  hash = randn(size(mt,2),1);
          %  hashM = mt*hash;
            hashV = (mt(var,:)*d)*hash;
            previous_var = find(abs(hashM - hashV) < 1e-20);
            if isempty(previous_var)
                newmt =  mt(getvariables(x),:)*d;
                mt(end+1,:) = newmt;
                yalmip('setmonomtable',mt,[variabletype newvariabletypegen(newmt)]);
                y = recover(size(mt,1));
            else
                y = recover(previous_var);
            end
        elseif  (size(base,2) == 2) & base(1)==0
            % Something like a*t^-d
            y = base(2)^d*recover(getvariables(x))^d;
        else
            % Bummer, something more complex, add an internal equality constraint
            y = (yalmip('define','mpower_internal',x))^d;
        end
    end
    y.extra.createTime = definecreationtime;
    y.extra.opname='';
    return
end

% Integer power of matrix
if x.dim(1)>1 | x.dim(2)>1
    switch d
        case 0
            y = 1;
        case 1
            y = x;
        otherwise
            y = x*mpower(x,d-1);
    end
else %Integer power of scalar

    base = x.basis;
    if isequal(base,[0 1]) % Unit scalar can be done fast
        [mt,variabletype,hashes,hash] = yalmip('monomtable');
        var = x.lmi_variables;
        newmt = mt(var,:)*d;
        newhash = newmt*hash;
        previous_var = findhash(hashes , newhash,length(hashes));
        if isempty(previous_var)
            mt = [mt;newmt];
            yalmip('setmonomtable',mt,[variabletype newvariabletypegen(newmt)]);
            y = x;
            y.lmi_variables = size(mt,1);
        else
            y = x;
            y.lmi_variables = previous_var;
        end

    else % General scalar
        switch d
            case 0
                y = 1;
            case 1
                y = x;
            case 2
                y = x*x;
            otherwise
                if even(d)
                    z = mpower(x,d/2);
                    y = z*z;
                else
                    y = x*mpower(x,d-1);
                end
        end
    end
    y.extra.createTime = definecreationtime;
    y.extra.opname='';
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [ x, w ] = rule01 ( n )

%*****************************************************************************80
%
%% RULE01 returns the rule of degree 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer N, the number of nodes.
%
%    Output, real X(2,N), the coordinates of the nodes.
%
%    Output, real W(N), the weights.
%
  xs = [ ...
       0.00000000000000000 ];
  ys = [ ...
       0.00000000000000000 ];
  ws = [ ...
       0.2828427124746189E+01 ];

  x(1,1:n) = xs(1:n);
  x(2,1:n) = ys(1:n);
  w(1:n) = ws(1:n);

  return
end


text:
function caffe_set_input(net, input)
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

input_blob_names = net.inputs;
input_size  = caffe_get_blobs_size(net, input_blob_names);
assert(numel(input_size) == numel(input));

caffe_set_blobs_data(net, input, input_blob_names);
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function w = map ( code, element_order )

%*****************************************************************************80
%
%% MAP returns the interpolation matrix for any available element.
%
%  Formula:
%
%    For an element of order ELEMENT_ORDER, we suppose we are given
%    ELEMENT_ORDER items of data Q associated with the nodes.
%
%   Let PHI(J)(R,S) be the Lagrange basis polynomial associated with
%   node J.  PHI(J)(R,S) is 1 at node J, and 0 at each of the other nodes.
%
%   Let P(R,S) be the polynomial of ELEMENT_ORDER terms which interpolates the
%   data Q, that is,
%
%      P(R(J),S(J)) = Q(J)
%
%   where the coordinates of node J are (R(J),S(J)).  Then we know
%   that we can write
%
%     P(R,S) = sum ( 1 <= J <= ELEMENT_ORDER ) Q(J) * PHI(J)(R,S)
%
%   But P(R,S) also has a standard representation as
%
%     P(R,S) = sum ( 1 <= I <= ELEMENT_ORDER ) A(I) * R**REXP(I) * S**SEXP(I)
%
%   where REXP(I) and SEXP(I) are the exponents of R and S and
%   the A(I) are the appropriate coefficients.
%
%   The interpolation matrix W allows us to immediately compute
%   the standard basis coefficients A from the data Q to be interpolated
%   using the formula:
%
%      A(I) = sum ( 1 <= J <= ELEMENT_ORDER ) W(I,J) * Q(J)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character CODE(*), identifies the element.
%    Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL',
%    'T3', 'T6' and 'T10'.
%
%    Input, integer ELEMENT_ORDER, the order associated with the code.
%
%    Output, real W(ELEMENT_ORDER,ELEMENT_ORDER), the interpolation matrix.
%

%
%  Get the (R,S) location of the nodes.
%
  [ r, s, area ] = node_reference ( code );
%
%  Get the associated monomials.
%
  [ rexp, sexp ] = poly ( code );
%
%  Set up the Vandermonde matrix.
%  Factors of the form 0**0 are to be understood as 1.
%
  for i = 1 : element_order
    for j = 1 : element_order

      if ( rexp(j) == 0 )
        rfact = 1.0;
      else
        rfact = r(i)^rexp(j);
      end

      if ( sexp(j) == 0 )
        sfact = 1.0;
      else
        sfact = s(i)^sexp(j);
      end

      w(i,j) = rfact * sfact;

    end
  end
%
%  Factor the Vandermonde matrix.
%
  [ w_lu, pivot, info ] = r8ge_fa ( element_order, w );

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MAP - Fatal error!\n' );
    fprintf ( 1, '  The Vandermonde matrix is singular.\n' );
    error ( 'MAP - Fatal error!' );
  end
%
%  Invert the Vandermonde matrix.
%
  w = r8ge_inverse ( element_order, w_lu, pivot );

  return
end


text:
% rf_blochSim.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% [mv,sc]=rf_blochSim(RF,tp,fspan,f0,peakB1,ph,M0);
%
% DESCRIPTION:
% Perform a bloch simulation of an RF pulse.  This code simply runs Martyn
% Klassen's excellent bloch equation simulator.  For more information, see
% help file for bes.m.  (~FID-A/rfPulseTools/mklassenTools/bes.m).
%
% INPUTS:
% RF        = RF pulse definition structure
% tp        = pulse duration in [ms]
% fspan     = Frequency span in [kHz] or, if the RF pulse includes a gradient
%             waveform (as indicated by having a 4th  column), then fspan
%             is the span of spatial positions in [cm] (optional.  Default=
%             10kHz or 10cm).
% f0        = Centre of frequnecy span [kHz] (optional.  Default=0)
% peakB1 	= Peak B1 amplitude in [kHz] (optional.  Default=RF.tw1/tp)
% ph        = Starting phase of the rf pulse [degrees] (optional. Default=0)
% M0        = Starting magnetization [units of M0] (optional. Default=[0,0,1])
%
% OUTPUTS:
% mv        = Simulated magnetization vector in three columns (x,y,z) as a
%             function of frequency.
% sc        = Frequency scale (in kHz), or if the pulse include a gradient
%             waveform, the position scale (in cm) corresponding to the
%             simulated mv vectors.

function [mv,sc]=rf_blochSim(RF,tp,fspan,f0,peakB1,ph,M0);

if nargin<7
    M0=[0,0,1];
    if nargin<6
        ph=0;
        if nargin<5
            peakB1=RF.tw1/tp;
            if nargin<4
                f0=0;
                if nargin<3
                    fspan=10;
                end
            end
        end
    end
end

[mv,sc]=bes(RF.waveform,tp,'f',peakB1,f0-fspan/2,f0+fspan/2,10000,ph,M0);

figure
subplot(4,1,1),plot(sc,mv(1,:),'LineWidth',1.2);
box off;
ylabel('M_x');

subplot(4,1,2),plot(sc,mv(2,:),'LineWidth',1.2);
box off;
ylabel('M_y');

subplot(4,1,3),plot(sc,sqrt(mv(2,:).^2 + mv(1,:).^2),'LineWidth',1.2);
box off;
ylabel('M_x_y');

subplot(4,1,4),plot(sc,mv(3,:),'LineWidth',1.2);
box off;
ylabel('M_z');

if size(RF.waveform,2)<4
    xlabel('Frequency (kHz)');
else
    xlabel('Position (cm)');
end

set(gcf,'color','w');

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function X = lyapunov_symmetric_eig(V, lambda, C, tol)
% Solves AX + XA = C when A = A', as a pseudo-inverse, given eig(A).
%
% function X = lyapunov_symmetric_eig(V, lambda, C)
% function X = lyapunov_symmetric_eig(V, lambda, C, tol)
%
% Same as lyapunov_symmetric(A, C, [tol]), where A is symmetric, its
% eigenvalue decomposition [V, lambda] = eig(A, 'vector') is provided as
% input directly, and C is a single matrix of the same size as A.
%
% See also: lyapunov_symmetric sylvester lyap sylvester_nocheck

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Aug. 31, 2018.
% Contributors:
% Change log:

    % AX + XA = C  is equivalent to DY + YD = M with
    % Y = V'XV, M = V'CV and D = diag(lambda).
    M = V'*C*V;

    % W(i, j) = lambda(i) + lambda(j)
    W = bsxfun(@plus, lambda, lambda');

    % Normally, the solution Y is simply this:
    Y = M ./ W;

    % But this may involve divisions by (almost) 0 in certain places.
    % Thus, we go for a pseudo-inverse.
    absW = abs(W);
    if ~exist('tol', 'var') || isempty(tol)
        tol = numel(C)*eps(max(absW(:))); % similar to pinv tolerance
    end
    Y(absW <= tol) = 0;

    % Undo the change of variable
    X = V*Y*V';

end

author: NicolasBoumal

repo: manopt


text:
function pyra = demo_deep_pyramid(im)
% Demonstrate basic usage and visualize features.
%
%
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the DeepPyramid code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

if exist('caffe') ~= 3
  error('You must add matcaffe to your path.');
end

if ~exist('data/caffe_nets/ilsvrc_2012_train_iter_310k')
  error(['You need the CNN model in %s. ' ...
         'You can get this model by following ' ...
         'the R-CNN installation instructions.'], ...
        'data/caffe_nets/ilsvrc_2012_train_iter_310k');
end

if 1
  % real use settings (compute features using the GPU)
  USE_GPU = true;
  USE_CACHE = false;
  USE_CAFFE = true;
else
  % fast demo settings
  USE_GPU = false;
  USE_CACHE = true;
  USE_CAFFE = false;
end

if ~exist('im', 'var') || isempty(im)
  im = imread('000084.jpg');
end
bbox = [263 145 381 225];

cnn = init_cnn_model('use_gpu', USE_GPU, 'use_caffe', USE_CAFFE);

if USE_CACHE
  cache_opts.cache_file = './cached_pyra.mat';
  cache_opts.debug = true;
  cache_opts.write_on_miss = true;
  th = tic;
  pyra = deep_pyramid_cache_wrapper(im, cnn, cache_opts);
  fprintf('deep_pyramid_cache_wrapper took %.3fs\n', toc(th));
else
  th = tic;
  pyra = deep_pyramid(im, cnn);
  fprintf('deep_pyramid took %.3fs\n', toc(th));
end
padx = 0;
pady = 0;

pyra = deep_pyramid_add_padding(pyra, padx, pady);

fprintf(['Press almost any key (with fig focused) to loop through ' ...
         'feature channels (or esc to exit).\n']);
for channel = 1:256
  vis_pyramid(im, pyra, bbox, channel);
  [~, ~, key_code] = ginput(1);
  if key_code == 27
    break;
  end
end

% ------------------------------------------------------------------------
function vis_pyramid(im, pyra, bbox, channel)
% ------------------------------------------------------------------------

pyra_boxes = im_to_pyra_coords(pyra, bbox);

clf;

rows = 2;
cols = 4;
subplot(rows, cols, 1);
imagesc(im);
axis image;
rectangle('Position', bbox_to_xywh(bbox), 'EdgeColor', 'g');
title(sprintf('input image feature %d', channel));

max_val = 0;
for level = 1:pyra.num_levels
  f = pyra.feat{level}(:,:,channel);
  max_val = max(max_val, max(f(:)));
end

ld = load('green_colormap');
colormap(ld.map); clear ld;

for level = 1:pyra.num_levels
  subplot(rows, cols, level+1);
  imagesc(pyra.feat{level}(:,:,channel), [0 max_val]);
  axis image;
  rectangle('Position', bbox_to_xywh(pyra_boxes{level}), 'EdgeColor', 'r');
  title(sprintf('level %d; scale = %.2fx', level, pyra.scales(level)));

  % project pyramid box back to image and display as red
  im_bbox = pyra_to_im_coords(pyra, [pyra_boxes{level} level]);
  subplot(rows, cols, 1);
  rectangle('Position', bbox_to_xywh(im_bbox), 'EdgeColor', 'r');
  %text(im_bbox(1), im_bbox(2), sprintf('%d', level));
end

function xywh = bbox_to_xywh(bbox)
xywh = [bbox(1) bbox(2) bbox(3)-bbox(1)+1 bbox(4)-bbox(2)+1];

author: rbgirshick

repo: DeepPyramid


text:
function sigma = snr2sigma(X,SNR)
%SNR2SIGMA Summary of this function goes here
%   Detailed explanation goes here
[N,B] = size(X);
p = mean(sum(X.^2,2));
sigma2 = p/(B*10^(SNR/10));
sigma = sqrt(sigma2);

end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function S1 = my_conv2(S1, sig, varargin)
% takes an extra argument which specifies which dimension to filter on
% extra argument can be a vector with all dimensions that need to be
% smoothed, in which case sig can also be a vector of different smoothing
% constants

if sig>.25
    idims = 2;
    if ~isempty(varargin)
        idims = varargin{1};
    end
    if numel(idims)>1 && numel(sig)>1
        sigall = sig;
    else
        sigall = repmat(sig, numel(idims), 1);
    end

    for i = 1:length(idims)
        sig = sigall(i);

        idim = idims(i);
        Nd = ndims(S1);

        S1 = permute(S1, [idim 1:idim-1 idim+1:Nd]);

        S1 = my_conv(S1, sig);
%         dsnew = size(S1);
%
%         S1 = reshape(S1, size(S1,1), []);
%         dsnew2 = size(S1);
%
%         tmax = ceil(4*sig);
%         dt = -tmax:1:tmax;
%         gaus = exp( - dt.^2/(2*sig^2));
%         gaus = gaus'/sum(gaus);
%
%         cNorm = filter(gaus, 1, cat(1, ones(dsnew2(1), 1), zeros(tmax,1)));
%         cNorm = cNorm(1+tmax:end, :);
%
%         S1 = filter(gaus, 1, cat(1, S1, zeros([tmax, dsnew2(2)])));
%         S1(1:tmax, :) = [];
%         S1 = reshape(S1, dsnew);
%
%         S1 = bsxfun(@rdivide, S1, cNorm);

        S1 = permute(S1, [2:idim 1 idim+1:Nd]);
    end
end

author: cortex-lab


text:
function [theta] = trainLinearReg(X, y, lambda)
%TRAINLINEARREG Trains linear regression given a dataset (X, y) and a
%regularization parameter lambda
%   [theta] = TRAINLINEARREG (X, y, lambda) trains linear regression using
%   the dataset (X, y) and regularization parameter lambda. Returns the
%   trained parameters theta.
%

% Initialize Theta
initial_theta = zeros(size(X, 2), 1);
% initial_theta = ones(size(X, 2), 1);

% Create "short hand" for the cost function to be minimized
costFunction = @(t) linearRegCostFunction(X, y, t, lambda);

% Now, costFunction is a function that takes in only one argument
options = optimset('MaxIter', 200, 'GradObj', 'on');

% Minimize using fmincg
theta = fmincg(costFunction, initial_theta, options);

end

repo: machine-learning-coursera


text:
function [rxnInLoops, N, loopInfo] = findMinNull(model, formulation, varargin)
% Find a minimal null-space for all internal cycles by solving a MILP,
% as proposed in Chan et al., 2017.
%
% USAGE:
%    [rxnInLoops, N, loopInfo] = findMinNull(model, formulation, parameters)
%
% INPUT:
%    model:               COBRA model structure
%
% OPTIONAL INPUTS:
%    formulation:         1: solve the MILP problem for a minimal null-space basis
%                            in an interlaced fashion by presolving some relaxed LPs to simultaneously
%                            determine the directions of reactions that participate in internal cycles
%                            (default, the quickest way as of 2017 Nov)
%                         2: directly solve the MILP. Then determine the directions
%                            of reactions that participate in internal cycles
%    parameters:          solver-specific parameter structure or name-value pair argument for solverCobraMILP
%
% OUTPUTS:
%    rxnInLoops:          #rxns-by-2 matrix. rxnInLoops(j, 1) = true => reverse direction of rxn j in loops
%                         rxnInLoops(j, 2) = true => forward direction of rxn j in loops
%    N:                   Minimal feasible null-space matrix for internal cycles
%    loopInfo:            structure containing the following parameters/information:
%                         *.M:          the bound for minimum/maximum flux
%                         *.minFlux:    minimum flux required for a reaction to be active
%                         *.ignoreRxns: rxns with small coefficients that are prechecked
%                                       before solving the MILP to avoid numerical issues
%                         *.nsTime:      wall time for finding the null-space
%                         *.nsCPU:      CPU time for finding the null-space
%                         *.loopPreprocessCPU:  CPU time for finding the directions of reactions participating in loops
%                         *.loopPreprocessTime: wall time for finding the directions of reactions participating in loops

if nargin < 2 || isempty(formulation)
	formulation = 1;
else
	if ~isscalar(formulation)
        % first varargin recognized as formulation
		varargin = [{formulation}; varargin(:)];
		formulation = 1;
	end
end

% filter out reactions with very small stoichiometric coefficents. They are
% usually not in loops. This can accelerate MILP solution time.
% (part of nullspace preprocessing)
cpu0 = cputime;
t = tic;

bigM = 100;
minFlux = 0.1;
rxnIn = sum(model.S ~= 0, 1) > 1;
rxnInIDs = find(rxnIn);
nR = sum(rxnIn);
metIn = any(model.S(:, rxnIn), 2);
nM = sum(metIn);
feasTol = getCobraSolverParams('LP', 'feasTol');

% rxnInLoops(j, 1) = true => reverse direction in loops
% rxnInLoops(j, 2) = true => forward direction in loops
rxnInLoops = false(numel(model.rxns), 2);

nsCPU = cputime - cpu0;
nsTime = toc(t);

if formulation == 1
    % Formulation 1: the procedure stated in SI Methods section 7 in Chan et al. (2017)
    % Determine reactions in loops and their feasible directions in loops
    % in an interlaced fashion

    % find maximum number of active reactions
    % (part of nullspace preprocessing)
    cpu0 = cputime;
    t = tic;

    MILP.A = [model.S(metIn, rxnIn),  sparse(nM, nR * 2); ... Sv = 0
             -speye(nR),             -bigM * speye(nR)    sparse(nR, nR); ...  -v - M z_pos <= -minFlux
              speye(nR),              sparse(nR, nR),    -bigM * speye(nR)];  % v - M z_neg <= -minFlux
    MILP.b = [zeros(nM, 1); -minFlux * ones(nR * 2, 1)];
    MILP.c = [zeros(nR, 1); ones(nR * 2, 1)];

    % Step 2 in SI Methods section 7.
    % Find all irreversible reactions that are in loops
    % fix z_pos = 1 if lb < 0, i.e., v may be negative
    % fix z_neg = 1 if ub > 0, i.e., v may be positive
    % (=> z_pos = z_neg = 1 for truly reversible reactions)
    MILP.lb = [-bigM * (model.lb(rxnIn) < 0); model.lb(rxnIn) < 0; model.ub(rxnIn) > 0];
    MILP.ub = [bigM * (model.ub(rxnIn) > 0); ones(nR * 2, 1)];
    MILP.osense = 1;
    MILP.csense = char(['E' * ones(1, nM), 'L' * ones(1, nR * 2)]);
    sol.irr = solveCobraLP(MILP, varargin{:});

    % Step 3 in SI Methods section 7.
    % Record identified reactions
    rxnInLoops(rxnIn, 1) = rxnInLoops(rxnIn, 1) | sol.irr.full(1:nR) < -1e-3;
    rxnInLoops(rxnIn, 2) = rxnInLoops(rxnIn, 2) | sol.irr.full(1:nR) > 1e-3;

    % Step 4 in SI Methods section 7.
    % Check reversible reactions whose forward direction not yet found to be in loops
    rxnCk = ~rxnInLoops(rxnIn, 2) & (model.lb(rxnIn) < 0 & model.ub(rxnIn) >0);
    % fix irreversible reactions known to be not in loops
    MILP.lb((model.lb(rxnIn) >= 0 | model.ub(rxnIn) <= 0) & ~any(rxnInLoops(rxnIn, :), 2)) = 0;
    MILP.ub((model.lb(rxnIn) >= 0 | model.ub(rxnIn) <= 0) & ~any(rxnInLoops(rxnIn, :), 2)) = 0;
    % fix z_pos = 1 for reactions not in rxnCk
    MILP.lb((nR + 1):end) = [~rxnCk; ones(nR, 1)];
    sol.rev1 = solveCobraLP(MILP, varargin{:});
    % Record identified reactions
    rxnInLoops(rxnIn, 1) = rxnInLoops(rxnIn, 1) | sol.rev1.full(1:nR) < -1e-3;
    rxnInLoops(rxnIn, 2) = rxnInLoops(rxnIn, 2) | sol.rev1.full(1:nR) > 1e-3;

    % Step 5 in SI Methods section 7.
    % Check reversible reactions whose reverse direction not yet found to be in loops
    rxnCk = ~rxnInLoops(rxnIn, 1) & (model.lb(rxnIn) < 0 & model.ub(rxnIn) >0);
    % fix z_neg = 1 for reactions not in rxnCk
    MILP.lb((nR + 1):end) = [ones(nR, 1); ~rxnCk];
    sol.rev2 = solveCobraLP(MILP, varargin{:});
    % Record identified reactions
    rxnInLoops(rxnIn, 1) = rxnInLoops(rxnIn, 1) | sol.rev2.full(1:nR) < -1e-3;
    rxnInLoops(rxnIn, 2) = rxnInLoops(rxnIn, 2) | sol.rev2.full(1:nR) > 1e-3;

    % This step is not mentioned in SI.
    % Filter reactions with small coefficients before solving MILP.
    % They are usually not in loops. This can accelerate MILP solution time.
    % reactions to be checked
    rxnCk = find(any(abs(model.S(:, rxnIn)) > 0 & abs(model.S(:, rxnIn)) < 1e-3, 1)' ...
        & model.lb(rxnIn) < 0 & model.ub(rxnIn) > 0 & ~any(rxnInLoops(rxnIn, :), 2));
    % fix all z_pos and z_neg
    MILP.lb((nR + 1):end) = 1;
    MILP.c(:) = 0;
    for j = 1:numel(rxnCk)
        % check feasibility in the forward direction
        if ~rxnInLoops(rxnInIDs(rxnCk(j)), 2) && MILP.ub(rxnCk(j)) > 0
            bd0 = MILP.lb(rxnCk(j));
            MILP.lb(rxnCk(j)) = minFlux / 10;
            solJ = solveCobraLP(MILP, varargin{:});
            if checkSolFeas(MILP, solJ) <= feasTol
                rxnInLoops(rxnIn, 1) = rxnInLoops(rxnIn, 1) | solJ.full(1:nR) < -1e-3;
                rxnInLoops(rxnIn, 2) = rxnInLoops(rxnIn, 2) | solJ.full(1:nR) > 1e-3;
            end
            MILP.lb(rxnCk(j)) = bd0;
        end
        % check feasibility in the reverse direction
        if ~rxnInLoops(rxnInIDs(rxnCk(j)), 1) && MILP.lb(rxnCk(j)) < 0
            bd0 = MILP.ub(rxnCk(j));
            MILP.ub(rxnCk(j)) = -minFlux / 10;
            solJ = solveCobraLP(MILP, varargin{:});
            if checkSolFeas(MILP, solJ) <= feasTol
                rxnInLoops(rxnIn, 1) = rxnInLoops(rxnIn, 1) | solJ.full(1:nR) < -1e-3;
                rxnInLoops(rxnIn, 2) = rxnInLoops(rxnIn, 2) | solJ.full(1:nR) > 1e-3;
            end
            MILP.ub(rxnCk(j)) = bd0;
        end
    end
    ignoreRxns = rxnCk(~any(rxnInLoops(rxnInIDs(rxnCk), :), 2));
    MILP.lb(ignoreRxns) = 0;
    MILP.ub(ignoreRxns) = 0;
    ignoreRxns = rxnInIDs(ignoreRxns);

    % Step 6 in SI Methods section 7.
    % Solve MILP for the forward directions
    rxnCk = ~rxnInLoops(rxnIn, 2) & (model.lb(rxnIn) < 0 & model.ub(rxnIn) >0);
    % fix z_pos for those not rxnCk. Fix all z_neg
    MILP.lb((nR + 1):end) = [~rxnCk; ones(nR, 1)];
    MILP.vartype = char(['C' * ones(1, nR), 'B' * ones(1, nR * 2)]);
    MILP.x0 = [];
    % minimize z_pos + z_neg
    MILP.c((nR + 1):end) = 1;
    sol.rev3 = solveCobraMILP(MILP, varargin{:});
    % Record identified reactions
    rxnInLoops(rxnIn, 1) = rxnInLoops(rxnIn, 1) | sol.rev3.full(1:nR) < -1e-3;
    rxnInLoops(rxnIn, 2) = rxnInLoops(rxnIn, 2) | sol.rev3.full(1:nR) > 1e-3;

    % end of the first part of null-space preprocessing
    nsCPU = nsCPU + (cputime - cpu0);
    nsTime = nsTime + toc(t);

    % Step 7 in SI Methods section 7.
    % Check whether the reactions with negative fluxes in sol.rev3.x can
    % actual have positive fluxes (part of looping reaction preprocessing)
    cpu0 = cputime;
    t = tic;

    rxnCk = find(sol.rev3.full(1:nR) < -1e-3 & rxnCk);
    if ~isempty(rxnCk)
        % fix all z when solving LP
        MILP.lb((nR + 1):end) = 1;
        MILP.c(:) = 0;
        for j = 1:numel(rxnCk)
            if ~rxnInLoops(rxnInIDs(rxnCk(j)), 2)
                bd0 = MILP.lb(rxnCk(j));
                MILP.lb(rxnCk(j)) = minFlux / 10;
                solJ = solveCobraLP(MILP, varargin{:});
                if checkSolFeas(MILP, solJ) <= feasTol
                    rxnInLoops(rxnIn, 1) = rxnInLoops(rxnIn, 1) | solJ.full(1:nR) < -1e-3;
                    rxnInLoops(rxnIn, 2) = rxnInLoops(rxnIn, 2) | solJ.full(1:nR) > 1e-3;
                end
                MILP.lb(rxnCk(j)) = bd0;
            end
        end
    end

    % end of the first part of loop preprocessing time
    loopPreprocessCPU = cputime - cpu0;
    loopPreprocessTime = toc(t);

    % Step 8 in SI Methods section 7.
    % Solve MILP for the reverse direction. (2nd part of nullspace preprocessing)
    cpu0 = cputime;
    t = tic;

    rxnCk = ~rxnInLoops(rxnIn, 1) & (model.lb(rxnIn) < 0 & model.ub(rxnIn) >0);
    % fix z_neg for those not rxnCk. Fix all z_pos
    MILP.lb((nR + 1):end) = [ones(nR, 1); ~rxnCk];
    MILP.c((nR + 1):end) = 1;
    sol.rev4 = solveCobraMILP(MILP, varargin{:});

    rxnInLoops(rxnIn, 1) = rxnInLoops(rxnIn, 1) | sol.rev4.full(1:nR) < -1e-3;
    rxnInLoops(rxnIn, 2) = rxnInLoops(rxnIn, 2) | sol.rev4.full(1:nR) > 1e-3;

    nsCPU = nsCPU + (cputime - cpu0);
    nsTime = nsTime + toc(t);

    % Step 9 in SI Methods section 7.
    % Check whether the reactions with positive fluxes in sol.rev4.cont can
    % actual have negative fluxes (part of looping reaction preprocessing)
    cpu0 = cputime;
    t = tic;

    rxnCk = find(sol.rev4.full(1:nR) > 1e-3 & rxnCk);
    if ~isempty(rxnCk)
        % fix all z when solving LP
        MILP.lb((nR + 1):end) = 1;
        MILP.c(:) = 0;
        for j = 1:numel(rxnCk)
            if ~rxnInLoops(rxnInIDs(rxnCk(j)), 1)
                bd0 = MILP.ub(rxnCk(j));
                MILP.ub(rxnCk(j)) = -minFlux / 10;
                solJ = solveCobraLP(MILP, varargin{:});
                if checkSolFeas(MILP, solJ) <= feasTol
                    rxnInLoops(rxnIn, 1) = rxnInLoops(rxnIn, 1) | solJ.full(1:nR) < -1e-3;
                    rxnInLoops(rxnIn, 2) = rxnInLoops(rxnIn, 2) | solJ.full(1:nR) > 1e-3;
                end
                MILP.ub(rxnCk(j)) = bd0;
            end
        end
    end
    loopPreprocessCPU = loopPreprocessCPU + (cputime - cpu0);
    loopPreprocessTime = loopPreprocessTime + toc(t);
else
    % Formulation 2: directly solve the MILP and find the feasible reaction
    % direction participating in loops for each reaction

    % find maximum number of active reactions (nullspace preprocessing)
    cpu0 = cputime;
    t = tic;

    % filter reactions with small coefficients before solving MILP
    % reactions to be checked
    LP.A = model.S(metIn, rxnIn);
    LP.b = zeros(nM, 1);
    LP.c = zeros(nR, 1);
    LP.lb = -bigM * (model.lb(rxnIn) < 0);
    LP.ub = bigM * (model.ub(rxnIn) > 0);
    LP.csense = char('E' * ones(1, nM));
    [lb0, ub0] = deal(LP.lb, LP.ub);
    % reactions to be checked
    rxnCk = find(any(abs(model.S(:, rxnIn)) > 0 & abs(model.S(:, rxnIn)) < 1e-3, 1));
    for j = 1:numel(rxnCk)
        % check feasibility in the forward direction
        if LP.ub(rxnCk(j)) > 0
            LP.lb(rxnCk(j)) = minFlux / 10;
            solJ = solveCobraLP(LP, varargin{:});
            if checkSolFeas(LP, solJ) <= feasTol
                rxnInLoops(rxnInIDs(rxnCk(j)), 2) = true;
            end
            LP.lb(rxnCk(j)) = lb0(rxnCk(j));
        end
        % check feasibility in the reverse direction
        if ~rxnInLoops(rxnInIDs(rxnCk(j)), 2) && LP.lb(rxnCk(j)) < 0
            LP.ub(rxnCk(j)) = -minFlux / 10;
            solJ = solveCobraLP(LP, varargin{:});
            if checkSolFeas(LP, solJ) <= feasTol
                rxnInLoops(rxnInIDs(rxnCk(j)), 1) = true;
            end
            LP.ub(rxnCk(j)) = ub0(rxnCk(j));
        end
    end
    clear LP
    ignoreRxns = rxnInIDs(rxnCk(~any(rxnInLoops(rxnInIDs(rxnCk), :), 2)));
    rxnIn(ignoreRxns) = false;
    nR = sum(rxnIn);
    metIn = any(model.S(:, rxnIn), 2);
    nM = sum(metIn);
    rxnInIDs = find(rxnIn);

    % solve MILP
    MILP.A = [model.S(metIn, rxnIn),  sparse(nM, nR * 2); ... Sv = 0
             -speye(nR),             -bigM * speye(nR)    sparse(nR, nR); ...  -v - M z_pos <= -1
              speye(nR),              sparse(nR, nR),    -bigM * speye(nR); ...   % v - M z_neg <= -1
              sparse(nR, nR),         speye(nR),          speye(nR)];   % z_pos + z_neg >= 1
    MILP.b = [zeros(nM, 1); -minFlux * ones(nR * 2, 1); ones(nR, 1)];
    MILP.c = [zeros(nR, 1); ones(nR * 2, 1)];
    % fix z_pos = 1 for reactions with ub <= 0. Fix z_neg = 1 for reactions with lb >= 0
    MILP.lb = [-bigM * (model.lb(rxnIn) < 0); model.ub(rxnIn) <= 0; model.lb(rxnIn) >= 0];
    MILP.ub = [bigM * (model.ub(rxnIn) > 0); ones(nR * 2, 1)];
    MILP.osense = 1;
    MILP.csense = char(['E' * ones(1, nM), 'L' * ones(1, nR * 2), 'G' * ones(1, nR)]);
    MILP.vartype = char(['C' * ones(1, nR), 'B' * ones(1, nR * 2)]);
    % z_pos for reactions with lb >= 0 and z_neg for reactions with ub <= 0 can be treated as continuous
    MILP.vartype(nR + find(model.lb(rxnIn) >= 0)) = 'C';
	MILP.vartype(nR * 2 + find(model.ub(rxnIn) <= 0)) = 'C';
    MILP.x0 = [zeros(nR, 1); ones(nR * 2, 1)];
    sol = solveCobraMILP(MILP, varargin{:});
    if isempty(sol.full)
		varargin2 = [{'feasTol'; 1e-7}; varargin(:)];
		sol = solveCobraMILP(MILP, varargin2{:});
    end

    nsCPU = nsCPU + (cputime - cpu0);
    nsTime = nsTime + toc(t);

    cpu0 = cputime;
    t = tic;

    % Check whether the reversible reactions with nonzero fluxes in sol.cont can
    % have fluxes in the opposite sign (part of looping reaction preprocessing)

    rxnInLoops(rxnIn, 1) = sol.full(1:nR) < -1e-3;
    rxnInLoops(rxnIn, 2) = sol.full(1:nR) > 1e-3;

    rxnCk = find(sol.full(1:nR) < -1e-3 | sol.full(1:nR) > 1e-3 & model.lb(rxnIn) < 0 & model.ub(rxnIn) > 0);

    LP.A = model.S(metIn, rxnIn);
    LP.b = zeros(nM, 1);
    LP.c = zeros(nR, 1);
    LP.lb = MILP.lb(1:nR) .* any(rxnInLoops(rxnIn, :), 2);
    LP.ub = MILP.ub(1:nR) .* any(rxnInLoops(rxnIn, :), 2);
    LP.csense = MILP.csense(1:nM);
    clear MILP

    for j = 1:numel(rxnCk)
        if ~rxnInLoops(rxnInIDs(rxnCk(j)), 2)
			bd0 = LP.lb(rxnCk(j));
            LP.lb(rxnCk(j)) = minFlux / 10;
            solJ = solveCobraLP(LP, varargin{:});
            if checkSolFeas(LP, solJ) <= feasTol
                rxnInLoops(rxnIn, 1) = rxnInLoops(rxnIn, 1) | solJ.full(1:nR) < -1e-3;
				rxnInLoops(rxnIn, 2) = rxnInLoops(rxnIn, 2) | solJ.full(1:nR) > 1e-3;
            end
            LP.lb(rxnCk(j)) = bd0;
        end
        if ~rxnInLoops(rxnInIDs(rxnCk(j)), 1)
			bd0 = LP.ub(rxnCk(j));
            LP.ub(rxnCk(j)) = -minFlux / 10;
            solJ = solveCobraLP(LP, varargin{:});
            if checkSolFeas(LP, solJ) <= feasTol
                rxnInLoops(rxnIn, 1) = rxnInLoops(rxnIn, 1) | solJ.full(1:nR) < -1e-3;
				rxnInLoops(rxnIn, 2) = rxnInLoops(rxnIn, 2) | solJ.full(1:nR) > 1e-3;
            end
            LP.ub(rxnCk(j)) = bd0;
        end
    end

    loopPreprocessCPU = cputime - cpu0;
    loopPreprocessTime = toc(t);
end

% get the nullspace matrix (part of nullspace preprocessing)
cpu0 = cputime;
t = tic;

rxnCyc = any(rxnInLoops, 2);
N = sparseNull(model.S(any(model.S(:, rxnCyc), 2), rxnCyc));
[row, col, entry] = find(N);
rxnCyc = find(rxnCyc);
row = rxnCyc(row);
N = sparse(row, col, entry, size(model.S, 2), size(N, 2));

nsCPU = nsCPU + (cputime - cpu0);
nsTime = nsTime + toc(t);

loopInfo.M = bigM;
loopInfo.minFlux = minFlux;
loopInfo.nsCPU = nsCPU;
loopInfo.nsTime = nsTime;
loopInfo.loopPreprocessCPU = loopPreprocessCPU;
loopInfo.loopPreprocessTime = loopPreprocessTime;
loopInfo.ignoreRxns = model.rxns(ignoreRxns);
end

author: opencobra

repo: cobratoolbox


text:
%CODEGENERATOR.FINSERTFRONT Insert a string at the beginning of a textfile.
%
% cGen.finsertfront(FNAME,NTEXT)
% FNAME is the full or relative path to the text file.
% NTEXT is the string containing the text to be inserted at the beginning
% of the file.
%
% Notes::
%  MatLab ships with functions for reading, overwriting and appending text
%  to files. This function is used to insert text at the beginning of a
%  text file.
%
% Authors::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
%  See also CodeGenerator.ffindreplace.

% Copyright (C) 1993-2012, by Peter I. Corke
% Copyright (C) 2012-2013, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function [ ] = finsertfront(CGen, fName,nText)

%% analyse the file
fid = fopen(fName,'r'); % open the file
% determine number of lines in the file
nLines = getNLines(fid);
% read all lines to cell
oldLines = cell(1,nLines);
for iLines = 1:nLines
    oldLines{iLines} = fgets(fid);
end
% close the file again
fclose(fid);

%% rewrite the file
fid = fopen(fName,'w');
fprintf(fid,'%s',nText);

for iLines = 1:nLines
   fprintf(fid,'%s',oldLines{iLines});
end
fclose(fid);

end

%% Determine the number of lines in a file
function [nLines] = getNLines(fid)
nLines = 0;
while (~feof(fid))      % go through each line until eof is reached
    fgets(fid);
    nLines = nLines +1; % increment line counter
end
frewind(fid);           % set file indicator back to the beginning of the file
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function ui = mrViewRemoveOverlay(ui,o);
%
% ui = mrViewRemoveOverlay([ui],[o]);
%
% Close/remove an overlay from a mrViewer UI.
% o is the index into the overlay; if omitted
% it removes the last overlay.
%
% ras 07/05.
if ~exist('ui','var') | isempty(ui), ui = mrViewGet;            end
if ~exist('o','var') | isempty(o), o = length(ui.overlays);     end

% new overlay order
newInd = setdiff(1:length(ui.overlays),o);

% close the UI panel if it exists
if checkfields(ui,'panels','overlays') & ishandle(ui.panels.overlays(o))
    % first check if it's in its own window: if it is, close the
    % window too
    pos = get(ui.panels.overlays(o),'Position');
    if isequal(pos,[0 0 1 1])
        par = get(ui.panels.overlays(o),'Parent');
        delete(par);
    else
        mrvPanelToggle(ui.panels.overlays(o),'off');
        delete(ui.panels.overlays(o));
    end
end

% remove the menu option to toggle this overlay
if checkfields(ui,'menus','overlayToggle') & ...
   ishandle(ui.menus.overlayToggle(o))
    delete(ui.menus.overlayToggle(o));
    ui.menus.overlayToggle = ui.menus.overlayToggle(newInd);
end

% remove the overlay field
ui.overlays = ui.overlays(newInd);

ui = mrViewRefresh(ui);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%islmatrix - Test if parameter is a logical matrix (>= 2 columns).
%
%  USAGE
%
%    test = islmatrix(x)
%
%    x              parameter to test
%
%  NOTE
%
%    To be considered logical, the matrix should contain only values 0 and 1, but it
%    does not need to actually be of class 'logical' (class(x) could be e.g. 'double').
%
%  SEE ALSO
%
%    See also islscalar, islvector, isdmatrix, isdvector, isdscalar, isivector, isiscalar,
%    isstring.
%

% Copyright (C) 2010-2011 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.

function test = islmatrix(x)

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help islmatrix'' for details).');
end

% Test: logical, two dimensions, two or more columns?
test = islogical(x) & length(size(x)) == 2 & size(x,2) >= 2;

function test = islogical(x)

test = builtin('islogical',x) | all(x(:)==0|x(:)==1);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%  Computes the hessian by finite differences
%
%  ::
%
%     H = finite_differences(Objective,params)
%     H = finite_differences(Objective,params,varargin)
%
%  Args:
%
%     - **Objective** [char|function handle]: function to differentiate
%     - **params** [vector]: point at which the differentiation is taken
%     - **varargin** : optional/further arguments of the objective function
%
%  Returns:
%     :
%
%     - **H** [matrix]: Hessian matrix
%
%  See also:
%
%     - utils.hessian.outer_product
%
%

author: jmaih


text:
function varargout = voxel(i,d,c,alpha)
% VOXEL Draw a 3-D voxel.
%    VOXEL(S) Draws a voxel centered at the specified point S. S is a
%       coordinate point in the form [ x y z ].
%    VOXEL(S,EDGE) Draws a voxel using the specified EDGE size. If no
%       EDGE is provided, it is set by default to 1. EDGE is a three
%       element vector [dx,dy,dz].
%    VOXEL(S,EDGE,C) Uses the specified colour C to draw the faces of the
%       voxel. C is a character string to specify color (type 'help plot'
%       to see list of valid colors). If no C is provided, it is set by
%       default to blue.
%    VOXEL(S,EDGE,C,ALPHA) Uses the given ALPHA to define the transparency
%       level (1 for opaque, 0 for transparent). If no ALPHA is given, it
%       is set to 1.
%    H = VOXEL(...) Return a vector of handles for the voxel drawn.
%
%
% See also PLOT, PATCH

%   Suresh Joel Apr 15,2003
%           Updated Feb 25, 2004

switch(nargin),
case 0
    disp('Too few arguements for voxel');
    return;
case 1
    d=[ 1 1 1 ]; %default length of side of voxel is 1
    c='b';       %default color of voxel is blue
    alpha=1;
case 2
    c='b';
    alpha=1;
case 3
    alpha=1;
end;

x=[i(1)+[0 0 0 0 d(1) d(1) d(1) d(1)]; ...
        i(2)+[0 0 d(2) d(2) 0 0 d(2) d(2)]; ...
        i(3)+[0 d(3) 0 d(3) 0 d(3) 0 d(3)]]';
h = [];
for n=1:3,
    if n==3,
        x=sortrows(x,[n,1]);
    else
        x=sortrows(x,[n n+1]);
    end;
    temp=x(3,:);
    x(3,:)=x(4,:);
    x(4,:)=temp;
    h1=patch(x(1:4,1),x(1:4,2),x(1:4,3),c);
    set(h1,'FaceAlpha',alpha);
    h = vertcat(h1,h);
    temp=x(7,:);
    x(7,:)=x(8,:);
    x(8,:)=temp;
    h1=patch(x(5:8,1),x(5:8,2),x(5:8,3),c);
    set(h1,'FaceAlpha',alpha);
    h = vertcat(h1,h);
end;

if nargout>0
  varargout{1} = h;
end

author: zouchuhang

repo: LayoutNet


text:
function agm_pi_test ( )

%*****************************************************************************80
%
%% AGM_PI_TEST tests the AGM_PI function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2012
%
%  Author:
%
%    John Burkardt
%
  clear all

  fprintf ( 1, '\n' );
  fprintf ( 1, 'AGM_PI_TEST\n' );
  fprintf ( 1, '  AGM_PI computes the value of pi using the \n' );
  fprintf ( 1, '  arithmetic-geometric-mean approach.\n' );
%
%  Get pi to 10, 20, 40, 80 decimal digits:
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Request PI to 10, 20, 40, 80 digits.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Digits   PI:\n' );
  fprintf ( 1, '\n' );

  d = 10;

  for i = 1 : 4

    p = agm_pi ( d );
    s = [ '  ', num2str(d,'%6d'), '  ', p ];
%
%  Trying to format output involving a symbolic variable seems to
%  be rather unsatisfactory.  FPRINTF won't work.  DISP can be used,
%  but you have to concatenate items first, and the concatenation operator
%  goes nuts when you have unexpected spaces in a phrase like "num2str ( d, '%6d' )"
%  and then DISP automatically prints the item PLUS an apparently unavoidable blank line
%  while suppressing all but a single blank when you inserted blanks so stuff
%  would line up, AND you also get the extraneous square brackets in your output.
%
    disp ( s )
    d = d * 2;

  end
%
%  Timings.
%  Apparently, D = 1 is illegal input to AGM_PI!
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  How long does it take to compute PI to D digits?\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '          Digits    Seconds\n' );
  fprintf ( 1, '\n' );
  d = 2;
  for i = 1 : 20
    tic
    p = agm_pi ( d );
    t = toc;
    fprintf ( 1, '  %2d  %10d  %14g\n', i, d, t );
    d = d * 2;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [f] = elec_fit_ellipse_optim(r, X, Y, Z, xo, yo, zo)

% elec_fit_ellipse_optim - Optimization for elec_fit_ellipse.m
%
% Called from elec_fit_ellipse.m
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:55 $

% Licence:  GNU GPL, no implied or express warranties
% History:  02/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% r is a 3x1 vector of radius values for each x,y,z axis component of ellipse
%
% equation of ellipsoid with center (xo,yo,zo) and radius for each axis (x,y,z) = (a,b,c):
% (( x - xo )^2 / a^2) + (( y - yo )^2 / b^2) + (( z - zo )^2 / c^2) = 1
%
% This function below creates a scalar value to
% return to the fminsearch function in elec_fit_sphere.

E = ( (X-xo).^2 )/r(1).^2 + ((Y-yo).^2)/r(2).^2 + ((Z-zo).^2)/r(3).^2  - 1;

f = sum( E .* E );  % sum of squares returned

return

author: PatternRecognition

repo: OpenBMI


text:
function gp = gplite_clean(gp)
%GPLITE_CLEAN Remove auxiliary info from lite GP struct (less memory usage).
%   GP = GPLITE_CLEAN(GP) removes auxiliary computational structs from
%   the GP. These can be reconstructed via a call to GPLITE_POST.
%
%   See also GPLITE_POST.

if ~isempty(gp) && isfield(gp,'post')
    copyfields = {'hyp'};
    emptyfields = {'alpha','sW','L','sn2_mult','Lchol'};
    checkfields = {'intmean'};
    for ff = copyfields; post0.(ff{:}) = []; end
    for ff = emptyfields; post0.(ff{:}) = []; end
    for ff = checkfields
        if isfield(gp.post(1),ff{:}); post0.(ff{:}) = []; end
    end

    for iG = 1:numel(gp)
        Ns = numel(gp(iG).post);
        postnew = post0;
        for iS = 1:Ns
            post_tmp = post0;
            for ff = copyfields
                post_tmp.(ff{:}) = gp(iG).post(iS).(ff{:});
            end
            postnew(iS) = post_tmp;
        end
        gp(iG).post = postnew;
    end
end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc