text:
%Vehicle Car-like vehicle class
%
% This class models the kinematics of a car-like vehicle (bicycle model) on
% a plane that moves in SE(2).  For given steering and velocity inputs it
% updates the true vehicle state and returns noise-corrupted odometry
% readings.
%
% Methods::
%   init         initialize vehicle state
%   f            predict next state based on odometry
%   step         move one time step and return noisy odometry
%   control      generate the control inputs for the vehicle
%   update       update the vehicle state
%   run          run for multiple time steps
%   Fx           Jacobian of f wrt x
%   Fv           Jacobian of f wrt odometry noise
%   gstep        like step() but displays vehicle
%   plot         plot/animate vehicle on current figure
%   plot_xy      plot the true path of the vehicle
%   add_driver   attach a driver object to this vehicle
%   display      display state/parameters in human readable form
%   char         convert to string
%
% Class methods::
%   plotv        plot/animate a pose on current figure
%
% Properties (read/write)::
%   x               true vehicle state: x, y, theta (3x1)
%   V               odometry covariance (2x2)
%   odometry        distance moved in the last interval (2x1)
%   rdim             dimension of the robot (for drawing)
%   L               length of the vehicle (wheelbase)
%   alphalim        steering wheel limit
%   maxspeed        maximum vehicle speed
%   T               sample interval
%   verbose         verbosity
%   x_hist          history of true vehicle state (Nx3)
%   driver          reference to the driver object
%   x0              initial state, restored on init()
%
% Examples::
%
% Create a vehicle with odometry covariance
%       v = Vehicle( diag([0.1 0.01].^2 );
% and display its initial state
%       v
% now apply a speed (0.2m/s) and steer angle (0.1rad) for 1 time step
%       odo = v.update([0.2, 0.1])
% where odo is the noisy odometry estimate, and the new true vehicle state
%       v
%
% We can add a driver object
%      v.add_driver( RandomPath(10) )
% which will move the vehicle within the region -10<x<10, -10<y<10 which we
% can see by
%      v.run(1000)
% which shows an animation of the vehicle moving for 1000 time steps
% between randomly selected wayoints.
%
% Notes::
% - Subclasses the MATLAB handle class which means that pass by reference semantics
%   apply.
%
% Reference::
%
%   Robotics, Vision & Control, Chap 6
%   Peter Corke,
%   Springer 2011
%
% See also RandomPath, EKF.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% 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

classdef Vehicle < handle

    properties
        % state
        x           % true state (x,y,theta)
        x_hist      % x history

        % parameters
        L           % length of vehicle
        alphalim    % steering wheel limit
        maxspeed    % maximum speed
        dim         % dimension of the world -dim -> +dim in x and y
        rdim    % dimension of the robot
        dt           % sample interval
        V           % odometry covariance

        odometry    % distance moved in last interval

        verbose

        driver      % driver object
        x0          % initial state
    end

    methods

        function veh = Vehicle(V, varargin)
        %Vehicle Vehicle object constructor
        %
        % V = Vehicle(V_ACT, OPTIONS)  creates a Vehicle object with actual odometry
        % covariance V_ACT (2x2) matrix corresponding to the odometry vector [dx dtheta].
        %
        % Options::
        % 'stlim',A       Steering angle limited to -A to +A (default 0.5 rad)
        % 'vmax',S        Maximum speed (default 5m/s)
        % 'L',L           Wheel base (default 1m)
        % 'x0',x0         Initial state (default (0,0,0) )
        % 'dt',T          Time interval
        % 'rdim',R        Robot size as fraction of plot window (default 0.2)
        % 'verbose'       Be verbose
        %
        % Notes::
        % - Subclasses the MATLAB handle class which means that pass by reference semantics
        %   apply.

            if ~isnumeric(V)
                error('first arg is V');
            end
            veh.x = zeros(3,1);
            if nargin < 1
                V = zeros(2,2);
            end

            opt.stlim = 0.5;
            opt.vmax = 5;
            opt.L = 1;
            opt.rdim = 0.2;
            opt.dt = 0.1;
            opt.x0 = zeros(3,1);
            opt = tb_optparse(opt, varargin);

            veh.V = V;

            veh.dt = opt.dt;
            veh.alphalim = opt.stlim;
            veh.maxspeed = opt.vmax;
            veh.L = opt.L;
            veh.x0 = opt.x0(:);
            veh.rdim = opt.rdim;
            veh.verbose = opt.verbose;

            veh.x_hist = [];
        end

        function init(veh, x0)
            %Vehicle.init Reset state of vehicle object
            %
            % V.init() sets the state V.x := V.x0, initializes the driver
            % object (if attached) and clears the history.
            %
            % V.init(X0) as above but the state is initialized to X0.
            if nargin > 1
                veh.x = x0(:);
            else
                veh.x = veh.x0;
            end
            veh.x_hist = [];
            if ~isempty(veh.driver)
                veh.driver.init()
            end
        end

        function add_driver(veh, driver)
            %Vehicle.add_driver Add a driver for the vehicle
            %
            % V.add_driver(D) connects a driver object D to the vehicle.  The driver
            % object has one public method:
            %        [speed, steer] = D.demand();
            % that returns a speed and steer angle.
            %
            % Notes::
            % - The Vehicle.step() method invokes the driver if one is attached.
            %
            % See also Vehicle.step, RandomPath.
            veh.driver = driver;
            driver.veh = veh;
        end

        function xnext = f(veh, x, odo, w)
            %Vehicle.f Predict next state based on odometry
            %
            % XN = V.f(X, ODO) is the predicted next state XN (1x3) based on current
            % state X (1x3) and odometry ODO (1x2) = [distance, heading_change].
            %
            % XN = V.f(X, ODO, W) as above but with odometry noise W.
            %
            % Notes::
            % - Supports vectorized operation where X and XN (Nx3).
            if nargin < 4
                w = [0 0];
            end

            dd = odo(1) + w(1); dth = odo(2);

            % straightforward code:
            % thp = x(3) + dth;
            % xnext = zeros(1,3);
            % xnext(1) = x(1) + (dd + w(1))*cos(thp);
            % xnext(2) = x(2) + (dd + w(1))*sin(thp);
            % xnext(3) = x(3) + dth + w(2);
            %
            % vectorized code:

            thp = x(:,3) + dth;
            xnext = x + [(dd+w(1))*cos(thp)  (dd+w(1))*sin(thp) ones(size(x,1),1)*dth+w(2)];
        end

        function odo = update(veh, u)
            %Vehicle.update Update the vehicle state
            %
            % ODO = V.update(U) is the true odometry value for
            % motion with U=[speed,steer].
            %
            % Notes::
            % - Appends new state to state history property x_hist.
            % - Odometry is also saved as property odometry.

            xp = veh.x; % previous state
            veh.x(1) = veh.x(1) + u(1)*veh.dt*cos(veh.x(3));
            veh.x(2) = veh.x(2) + u(1)*veh.dt*sin(veh.x(3));
            veh.x(3) = veh.x(3) + u(1)*veh.dt/veh.L * u(2);
            odo = [colnorm(veh.x(1:2)-xp(1:2)) veh.x(3)-xp(3)];
            veh.odometry = odo;

            veh.x_hist = [veh.x_hist; veh.x'];   % maintain history
        end

        function J = Fx(veh, x, odo)
        %Vehicle.Fx  Jacobian df/dx
        %
        % J = V.Fx(X, ODO) is the Jacobian df/dx (3x3) at the state X, for
        % odometry input ODO (1x2) = [distance, heading_change].
        %
        % See also Vehicle.f, Vehicle.Fv.
            dd = odo(1); dth = odo(2);
            thp = x(3) + dth;

            J = [
                1   0   -dd*sin(thp)
                0   1   dd*cos(thp)
                0   0   1
                ];
        end

        function J = Fv(veh, x, odo)
            %Vehicle.Fv  Jacobian df/dv
            %
            % J = V.Fv(X, ODO) is the Jacobian df/dv (3x2) at the state X, for
            % odometry input ODO (1x2) = [distance, heading_change].
            %
            % See also Vehicle.F, Vehicle.Fx.
            dd = odo(1); dth = odo(2);
            thp = x(3) + dth;

            J = [
                cos(thp)    -dd*sin(thp)
                sin(thp)    dd*cos(thp)
                0           1
                ];
        end

        function odo = step(veh, varargin)
            %Vehicle.step Advance one timestep
            %
            % ODO = V.step(SPEED, STEER) updates the vehicle state for one timestep
            % of motion at specified SPEED and STEER angle, and returns noisy odometry.
            %
            % ODO = V.step() updates the vehicle state for one timestep of motion and
            % returns noisy odometry.  If a "driver" is attached then its DEMAND() method
            % is invoked to compute speed and steer angle.  If no driver is attached
            % then speed and steer angle are assumed to be zero.
            %
            % Notes::
            % - Noise covariance is the property V.
            %
            % See also Vehicle.control, Vehicle.update, Vehicle.add_driver.

            % get the control input to the vehicle from either passed demand or driver
            u = veh.control(varargin{:});

            % compute the true odometry and update the state
            odo = veh.update(u);

            % add noise to the odometry
            if veh.V
                odo = veh.odometry + randn(1,2)*veh.V;
            end
        end

        function u = control(veh, speed, steer)
            %Vehicle.control Compute the control input to vehicle
            %
            % U = V.control(SPEED, STEER) is a control input (1x2) = [speed,steer]
            % based on provided controls SPEED,STEER to which speed and steering angle
            % limits have been applied.
            %
            % U = V.control() as above but demand originates with a "driver" object if
            % one is attached, the driver's DEMAND() method is invoked. If no driver is
            % attached then speed and steer angle are assumed to be zero.
            %
            % See also Vehicle.step, RandomPath.
            if nargin < 2
                % if no explicit demand, and a driver is attached, use
                % it to provide demand
                if ~isempty(veh.driver)
                    [speed, steer] = veh.driver.demand();
                else
                    % no demand, do something safe
                    speed = 0;
                    steer = 0;
                end
            end

            % clip the speed
            u(1) = min(veh.maxspeed, max(-veh.maxspeed, speed));

            % clip the steering angle
            u(2) = max(-veh.alphalim, min(veh.alphalim, steer));
        end

        function p = run(veh, nsteps)
            %Vehicle.run Run the vehicle simulation
            %
            % V.run(N) runs the vehicle model for N timesteps and plots
            % the vehicle pose at each step.
            %
            % P = V.run(N) runs the vehicle simulation for N timesteps and
            % return the state history (Nx3) without plotting.  Each row
            % is (x,y,theta).
            %
            % See also Vehicle.step.

            if nargin < 2
                nsteps = 1000;
            end

            %veh.clear();
            if ~isempty(veh.driver)
                veh.driver.visualize();
            end

            veh.visualize();
            for i=1:nsteps
                veh.step();
                if nargout == 0
                    % if no output arguments then plot each step
                    veh.plot();
                    drawnow
                end
            end
            p = veh.x_hist;
        end

        % TODO run and run2 should become superclass methods...

        function p = run2(veh, T, x0, speed, steer)
            %Vehicle.run2 Run the vehicle simulation with control inputs
            %
            % P = V.run2(T, X0, SPEED, STEER) runs the vehicle model for a time T with
            % speed SPEED and steering angle STEER.  P (Nx3) is the path followed and
            % each row is (x,y,theta).
            %
            % Notes::
            % - Faster and more specific version of run() method.
            % - Used by the RRT planner.
            %
            % See also Vehicle.run, Vehicle.step, RRT.
            veh.init(x0);

            for i=1:(T/veh.dt)
                veh.update([speed steer]);
            end
            p = veh.x_hist;
        end

        function h = plot(veh, varargin)
        %Vehicle.plot Plot vehicle
        %
        % V.plot(OPTIONS) plots the vehicle on the current axes at a pose given by
        % the current state.  If the vehicle has been previously plotted its
        % pose is updated.  The vehicle is depicted as a narrow triangle that
        % travels "point first" and has a length V.rdim.
        %
        % V.plot(X, OPTIONS) plots the vehicle on the current axes at the pose X.
        %
        % H = V.plotv(X, OPTIONS) draws a representation of a ground robot as an
        % oriented triangle with pose X (1x3) [x,y,theta].  H is a graphics handle.
        %
        % V.plotv(H, X) as above but updates the pose of the graphic represented
        % by the handle H to pose X.
        %
        % Options::
        % 'scale',S    Draw vehicle with length S x maximum axis dimension
        % 'size',S     Draw vehicle with length S
        % 'color',C    Color of vehicle.
        % 'fill'       Filled
        %
        % See also Vehicle.plotv.

            h = findobj(gcf, 'Tag', 'Vehicle.plot');
            if isempty(h)
                % no instance of vehicle graphical object found
                h = Vehicle.plotv(veh.x, varargin{:});
                set(h, 'Tag', 'Vehicle.plot');  % tag it
            end

            if ~isempty(varargin) && isnumeric(varargin{1})
                % V.plot(X)
                Vehicle.plotv(h, varargin{1}); % use passed value
            else
                % V.plot()
                Vehicle.plotv(h, veh.x);    % use current state
            end
        end

        function out = plot_xy(veh, varargin)
            %Vehicle.plot_xy Plots true path followed by vehicle
            %
            % V.plot_xy() plots the true xy-plane path followed by the vehicle.
            %
            % V.plot_xy(LS) as above but the line style arguments LS are passed
            % to plot.
            %
            % Notes::
            % - The path is extracted from the x_hist property.

            xyt = veh.x_hist;
            if nargout == 0
                plot(xyt(:,1), xyt(:,2), varargin{:});
            else
                out = xyt;
            end
        end

        function visualize(veh)
            grid on
        end

        function verbosity(veh, v)
        %Vehicle.verbosity Set verbosity
        %
        % V.verbosity(A) set verbosity to A.  A=0 means silent.
            veh.verbose = v;
        end

        function display(nav)
        %Vehicle.display Display vehicle parameters and state
        %
        % V.display() displays vehicle parameters and state in compact
        % human readable form.
        %
        % Notes::
        % - This method is invoked implicitly at the command line when the result
        %   of an expression is a Vehicle object and the command has no trailing
        %   semicolon.
        %
        % See also Vehicle.char.

            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            disp( char(nav) );
        end % display()

        function s = char(veh)
        %Vehicle.char Convert to a string
        %
        % s = V.char() is a string showing vehicle parameters and state in
        % a compact human readable format.
        %
        % See also Vehicle.display.

            s = 'Vehicle object';
            s = char(s, sprintf(...
            '  L=%g, maxspeed=%g, alphalim=%g, T=%f, nhist=%d', ...
                veh.L, veh.maxspeed, veh.alphalim, veh.dt, ...
                numrows(veh.x_hist)));
            if ~isempty(veh.V)
                s = char(s, sprintf(...
                '  V=(%g,%g)', ...
                    veh.V(1,1), veh.V(2,2)));
            end
            s = char(s, sprintf('  x=%g, y=%g, theta=%g', veh.x));
            if ~isempty(veh.driver)
                s = char(s, '  driven by::');
                s = char(s, [['    '; '    '] char(veh.driver)]);
            end
        end

    end % method

    methods(Static)

        function h_ = plotv(x, varargin)
        %Vehicle.plotv Plot ground vehicle pose
        %
        % H = Vehicle.plotv(X, OPTIONS) draws a representation of a ground robot as an
        % oriented triangle with pose X (1x3) [x,y,theta].  H is a graphics handle.
        % If X (Nx3) is a matrix it is considered to represent a trajectory in which case
        % the vehicle graphic is animated.
        %
        % Vehicle.plotv(H, X) as above but updates the pose of the graphic represented
        % by the handle H to pose X.
        %
        % Options::
        % 'scale',S    Draw vehicle with length S x maximum axis dimension
        % 'size',S     Draw vehicle with length S
        % 'color',C    Color of vehicle.
        % 'fill'       Filled with solid color as per 'color' option
        % 'fps',F      Frames per second in animation mode (default 10)
        %
        % Example::
        %
        % Generate some path 3xN
        %         p = PRM.plan(start, goal);
        % Set the axis dimensions to stop them rescaling for every point on the path
        %         axis([-5 5 -5 5]);
        %
        % Now invoke the static method
        %         Vehicle.plotv(p);
        %
        % Notes::
        % - This is a class method.
        %
        % See also Vehicle.plot.

            if isscalar(x) && ishandle(x)
                % plotv(h, x)
                h = x;
                x = varargin{1};
                x = x(:)';
                T = transl([x(1:2) 0]) * trotz( x(3) );
                set(h, 'Matrix', T);
                return
            end

            opt.scale = 1/60;
            opt.size = [];
            opt.fill = false;
            opt.color = 'r';
            opt.fps = 10;

            [opt,args] = tb_optparse(opt, varargin);

            lineprops = { 'Color', opt.color' };
            if opt.fill
                lineprops = [lineprops 'fill' opt.color ];
            end

            % compute the dimensions of the robot
            if ~isempty(opt.size)
                d = opt.size;
            else
                % get the current axes dimensions
                a = axis;
                d = (a(2)+a(4) - a(1)-a(3)) * opt.scale;
            end

            % draw it
            points = [
                d 0
                -d -0.6*d
                -d 0.6*d
            ]';

            h = hgtransform();
            hp = plot_poly(points, lineprops{:});
            for hh=hp
                set(hh, 'Parent', h);
            end

            if (numel(x) > 3) && (numcols(x) == 3)
                % animation mode
                for i=1:numrows(x)
                    T = transl([x(i,1:2) 0]) * trotz( x(i,3) );
                    set(h, 'Matrix', T);
                    pause(1/opt.fps);
                end
            elseif (numel(x) == 3)
                % compute the pose
                % convert vector form of pose to SE(3)

                x = x(:)';
                T = transl([x(1:2) 0]) * trotz( x(3) );
                set(h, 'Matrix', T);
            else
                error('bad pose');
            end

            if nargout > 0
                h_ = h;
            end
        end

    end % static methods

end % classdef

author: Allopart

repo: rbpf-gmapping


text:
function Rth=th_rot(phi,theta)
%TH_ROT Calculates the rotation matrix for the RF pulse flip angle, taking into
% consideration of the RF pulses phase.
%
%   phi: RF pulse flip-angle (radian)
%   theta: RF pulse phase (radian).

Rz = z_rot(-theta);
Rx = x_rot(phi);
Rth = inv(Rz)*Rx*Rz;

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [LR, HR, kernel, sigma] = degradation_model(HR, scale)

% kernel
ksize  = 15;
theta  = pi*rand(1);
l1     = 0.1+9.9*rand(1);
l2     = 0.1+(l1-0.1)*rand(1);
kernel = anisotropic_Gaussian(ksize,theta,l1,l2); % double
kernel = single(kernel);

% noise level
sigma_max = 75;
sigma = single(sigma_max*rand(1)/255); % single

% HR image
HR  = modcrop(HR, scale); % double

% xxxxxxxxxxxxxx degradation model xxxxxxxxxxxxxxxxxxxxxx
% you can change to your own degradation model

blurry_HR  = imfilter(HR,double(kernel),'replicate');
downsampled_HR = imresize(blurry_HR, 1/scale,'bicubic');
downsampled_HR = im2single(im2uint8(downsampled_HR));

noise  = single(randn(size(downsampled_HR),'single')*sigma);
LR = downsampled_HR + noise;

kernel = single(kernel);
% xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

end

author: cszn

repo: SRMD

save_path: github-repos/MATLAB/cszn-SRMD


text:
function [wave] = PL2Waves(filename, channel, unit)
% PL2Waves( filename, channel, unit ): read spike timestamps and waveforms from a .pl2 file
%                returns waveform values in millivolts
%
% wave = PL2Waves( filename, channel, unit )
% wave = PL2Waves( filename, 'SPK01', 1 );
% wave = PL2Waves( filename, 2, 3 );
%
% INPUT:
%   filename - if empty string, will use File Open dialog
%   channel - channel name or 1-based channel number
%           1-based channel number is the index of the spike channel in pl2.SpikeChannels cell array
%           shown in the Name column of the output produced by PL2Print(pl2.SpikeChannels)
%           To print the list of spike channels, use these commands:
%
%           pl2 = PL2GetFileIndex(filename);
%           PL2Print(pl2.SpikeChannels);
%
%   unit  - unit number (0 - unsorted, 1-254 sorted units, 255 - invalidated spikes)
%
% OUTPUT:
%   wave.NumPointsWave - number of points in each waveform
%   wave.Ts - array of timestamps (in seconds)
%   wave.Waves - matrix of waveforms, each row is a waveform in millivolts

wave = internalPL2Waves(filename, channel, unit);

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [trx,alltimestamps] = MergeMWTData(varargin)

MINDT = .0001;

trxall = varargin(1:3:end-2);
fnsall = varargin(2:3:end-1);
fnscopyall = varargin(3:3:end);
nsets = numel(trxall);
if nsets < 1,
  error('Usage: MergeMWTData(trx1,fns1,trx2,fns2,...)');
end

ids1 = unique([trxall{1}.id]);
ntargets = numel(ids1);
for seti = 2:nsets,
  ids2 = [trxall{seti}.id];
  if ~isempty(setxor(ids1,ids2)),
    error('ids do not match for sets 1 and %d',seti);
  end
end
idsall = cell(1,nsets);
for seti = 1:nsets,
  idsall{seti} = [trxall{seti}.id];
end

idx = cell(1,nsets);
for seti = 1:nsets,
  [~,idx{seti}] = ismember(ids1,idsall{seti});
end

alltimestamps = [];
for seti = 1:nsets,
  alltimestamps = union(alltimestamps,[trxall{seti}.timestamps]);
end
dts = diff(alltimestamps);
if any(dts < MINDT),
  error('Rounding error check failed: there are timestamps < %f apart',MINDT);
end

% frames per target
framesall = cell(ntargets,nsets);
firstframes = nan(ntargets,nsets);
endframes = nan(ntargets,nsets);
for seti = 1:nsets,
  for i = 1:ntargets,
    j = idx{seti}(i);
    [~,framesall{i,seti}] = ismember(trxall{seti}(j).timestamps,alltimestamps);
    firstframes(i,seti) = min(framesall{i,seti});
    endframes(i,seti) = max(framesall{i,seti});
  end
end

trx = [];
for i = 1:ntargets,
  trk = struct;
  trk.id = ids1(i);
  % union all frames for which there is data
  trk.firstframe = min(firstframes(i,:));
  trk.endframe = max(endframes(i,:));
  trk.nframes = trk.endframe - trk.firstframe + 1;
  trk.off = 1 - trk.firstframe;
  trk.timestamps = alltimestamps(trk.firstframe:trk.endframe);
  trk.dt = diff(trk.timestamps);
  for seti = 1:nsets,
    j = idx{seti}(i);
    frames = framesall{i,seti};
    for k = 1:numel(fnsall{seti}),
      fn = fnsall{seti}{k};
      sz = size(trxall{seti}(j).(fn));
      if numel(frames) ~= sz(2),
        error('Number of frames does not match size of %s for input %d',fn,seti);
      end
      if sz(2) ~= trk.nframes,
        fprintf('Adding %d empty frames to id %d, input %d, field %s\n',...
          trk.nframes-sz(2),ids1(i),fn);
      end
      sz(2) = trk.nframes;
      if iscell(trxall{seti}(j).(fn)),
        tmp = cell(sz);
      else
        tmp = nan(sz);
      end
      tmp(:,frames-trk.firstframe+1,:) = trxall{seti}(j).(fn);
      if isfield(trk,fn),
        if iscell(trk.(fn)) || iscell(tmp) || ...
            ndim(trk.(fn)) ~= numel(sz) || ...
            ~all(size(trk.(fn))==sz),
          warning('Field %s in multiple inputs and cannot be merged, skipping input %d',fn,seti);
        else
          idxdouble = ~isnan(trk.(fn)) & ~isnan(tmp);
          if any(trk.(fn)(idxdouble) ~= tmp(idxdouble)),
            warning('Field %s is in multiple inputs and does not match, skipping mismatches from input %d',fn,seti);
          end
          idxadd = isnan(trk.(fn)) & ~isnan(tmp);
          trk.(fn)(idxadd) = tmp(idxadd);
        end
      else
        trk.(fn) = tmp;
      end
    end
    for k = 1:numel(fnscopyall{seti}),
      fn = fnscopyall{seti}{k};
      trk.(fn) = trxall{seti}.(fn);
    end
  end
  trx = structappend(trx,trk);
end

% check for fields that got lost
fnsmerge = fieldnames(trx);
for seti = 1:nsets,
  fns = fieldnames(trxall{seti});
  missingfns = setdiff(fns,fnsmerge);
  if ~isempty(missingfns),
    fprintf('The following fields from input %d were not merged:',seti);
    fprintf(' %s',missingfns{:});
    fprintf('\n');
  end
end

author: kristinbranson

repo: JAABA


text:
function varargout = cfg_util_persistent(varargin)
%CFG_UTIL_PERSISTENT - store persistent variables for cfg_util
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_util_persistent.m 5750 2013-11-15 15:02:24Z volkmar $

rev = '$Rev: 5750 $';

persistent c0;
persistent jobs;
if nargin == 2 && nargout == 0
    c0   = varargin{1};
    jobs = varargin{2};
elseif nargin == 0 && nargout == 2
    varargout{1} = c0;
    varargout{2} = jobs;
else
    cfg_message('matlabbatch:usage', '%s: Usage error.', mfilename);
end

author: spm


text:
function [modelMin,AddedExchange] = findMinCardModel(model,Ex_Rxns)
% Function to find the minimal cardinality model.
%
% USAGE:
%    [modelMin, AddedExchange] = findMinCardModel(model, Ex_Rxns)
%
% INPUTS:
%     model:         Metabolic model
%     Ex_Rxns:       Vector of exchange reactions for FVA
%
% OUTPUTS:
%    modelMin:       Updated model with new reaction bounds
%    AddedExchange:  Vector of exchanged reactions
%
% .. Authors:
%       - Ines Thiele 2014
%       - Maike K. Aurich 27/05/15 change tolerance

tol = -1e-6;  % Default tolerance (limit for calling flux 0)

AddedExchange = '';
% convert to irrev model
[modelIrrev,matchRev,rev2irrev,irrev2rev] = convertToIrreversible(model);
%translate Ex_Rxns from rev model to irrev model
Ex_RxnsIrrev=[];
for w=1:length(Ex_Rxns)
    tmp = modelIrrev.rxns(strmatch(Ex_Rxns(w),modelIrrev.rxns));
    Ex_RxnsIrrev = [Ex_RxnsIrrev; tmp];
end

minNorm = zeros(length(modelIrrev.rxns),1);
minNorm(ismember(modelIrrev.rxns,Ex_RxnsIrrev))= 1;

modelIrrev.osense=1; %minimize linear objective

[solutionCard2]=solveCobraLPCPLEXcard(modelIrrev,0,0,[],[],minNorm,'zero');

% needs backward translation from irrev to rev
% update constraints in original model according to solution
RxnID = find(ismember(modelIrrev.rxns,Ex_RxnsIrrev));
unUsedRxns = find(solutionCard2.full<=abs(tol));
usedRxns = find(solutionCard2.full>abs(tol));
UnusedRxnID = intersect(unUsedRxns,RxnID);
UsedRxnID = intersect(usedRxns,RxnID);
modelMin = model;
for i = 1 :length(UnusedRxnID)
    if strcmp(modelIrrev.rxns{UnusedRxnID(i)}(end-1:end),'_r')
        OriModelName = modelIrrev.rxns{UnusedRxnID(i)}(1:end-2);
        modelMin = changeRxnBounds(modelMin,OriModelName,0,'b');
    elseif strcmp(modelIrrev.rxns{UnusedRxnID(i)}(end-1:end),'_b') % no uptake
        OriModelName = modelIrrev.rxns{UnusedRxnID(i)}(1:end-2);
        modelMin = changeRxnBounds(modelMin,OriModelName,0,'l');
    elseif strcmp(modelIrrev.rxns{UnusedRxnID(i)}(end-1:end),'_f') % no secretion
        OriModelName = modelIrrev.rxns{UnusedRxnID(i)}(1:end-2);
        modelMin = changeRxnBounds(modelMin,OriModelName,0,'u');
    else
        OriModelName = modelIrrev.rxns{UnusedRxnID(i)};
        modelMin = changeRxnBounds(modelMin,OriModelName,0,'b');
    end
end
cntR=1;
for i = 1 :length(UsedRxnID)
    if strcmp(modelIrrev.rxns{UsedRxnID(i)}(end-1:end),'_r')
        AddedExchange{cntR,1} = strcat(modelIrrev.rxns{UsedRxnID(i)}(1:end-2));
        cntR = cntR+1;
    elseif strcmp(modelIrrev.rxns{UsedRxnID(i)}(end-1:end),'_b') %  uptake
        AddedExchange{cntR,1}= strcat(modelIrrev.rxns{UsedRxnID(i)}(1:end-2));

        cntR = cntR+1;
    elseif strcmp(modelIrrev.rxns{UsedRxnID(i)}(end-1:end),'_f') %  uptake
        AddedExchange{cntR,1}= strcat(modelIrrev.rxns{UsedRxnID(i)}(1:end-2));

        cntR = cntR+1;

    else
        AddedExchange{cntR,1} = strcat(modelIrrev.rxns{UsedRxnID(i)});
        cntR = cntR+1;
    end
end

author: opencobra

repo: cobratoolbox


text:
function h = hmatrix ( p, xx, yy, dd, hh, varargin )

%*****************************************************************************80
%
%% HMATRIX computes the mesh size function from values specified on a Cartesian grid.
%
%  Licensing:
%
%    (C) 2004 Per-Olof Persson.
%    See COPYRIGHT.TXT for details.
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real P(NP,2), the point coordinates.
%
%    Input, real XX(NDATA), YY(NDATA), the coordinates of points at which the
%    mesh size function has been specified.
%
%    Input, real DD, a dummy parameter included for consistency with the DMATRIX routine.
%
%    Input, real HH(NDATA), the value of the mesh size function at the data points.
%
%    Input, VARARGIN, room for extra arguments.
%
%    Output, real H(NP), the interpolated value of the mesh size function at each
%    of the input points.
%
  h = interp2 ( xx, yy, hh, p(:,1), p(:,2), '*linear' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [indices, bestSqDist] = nearpoints2d(src, dst)
%
% [indices, bestSqDist] = nearpoints2d(src, dest)
%
% For each point in one set, find the nearest point in another.
%
%- src is a 2xM array of points
%- dest is a 2xN array of points
%- indices is a 1xM vector, in which each element tells which point
%  in dest is closest to the corresponding point in src.  For
%  example, dest[indices[i]] is near src[i].
%- bestSqDist is a 1xM array of the squared distance between
%  dest[indices[i]] and src[i].
%
% SEE ALSO: nearpoints
%
% HISTORY:
%  2006.05.05 RFD wrote it.

if(size(src,1)~=2) src = src'; end
if(size(dst,1)~=2) dst = dst'; end

if(size(src,1)~=2 || size(dst,1)~=2)
    error('arrays must be 2xN and 2xM!');
end

src = [src; zeros(size(src(1,:)))];
dst = [dst; zeros(size(dst(1,:)))];
[indices, bestSqDist] = nearpoints(src, dst);
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [X, info] = dominant_invariant_subspace(A, p)
% Returns an orthonormal basis of the dominant invariant p-subspace of A.
%
% function X = dominant_invariant_subspace(A, p)
%
% Input: A real, symmetric matrix A of size nxn and an integer p < n.
% Output: A real, orthonormal matrix X of size nxp such that trace(X'*A*X)
%         is maximized. That is, the columns of X form an orthonormal basis
%         of a dominant subspace of dimension p of A. These are thus
%         eigenvectors associated with the largest eigenvalues of A (in no
%         particular order). Sign is important: 2 is deemed a larger
%         eigenvalue than -5.
%
% The optimization is performed on the Grassmann manifold, since only the
% space spanned by the columns of X matters. The implementation is short to
% show how Manopt can be used to quickly obtain a prototype. To make the
% implementation more efficient, one might first try to use the caching
% system, that is, use the optional 'store' arguments in the cost, grad and
% hess functions. Furthermore, using egrad2rgrad and ehess2rhess is quick
% and easy, but not always efficient. Having a look at the formulas
% implemented in these functions can help rewrite the code without them,
% possibly more efficiently.
%
% See also: dominant_invariant_subspace_complex

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Nicolas Boumal, July 5, 2013
% Contributors:
%
% Change log:
%
%   NB Dec. 6, 2013:
%       We specify a max and initial trust region radius in the options.

    % Generate some random data to test the function
    if ~exist('A', 'var') || isempty(A)
        A = randn(128);
        A = (A+A')/2;
    end
    if ~exist('p', 'var') || isempty(p)
        p = 3;
    end

    % Make sure the input matrix is square and symmetric
    n = size(A, 1);
	assert(isreal(A), 'A must be real.')
    assert(size(A, 2) == n, 'A must be square.');
    assert(norm(A-A', 'fro') < n*eps, 'A must be symmetric.');
	assert(p<=n, 'p must be smaller than n.');

    % Define the cost and its derivatives on the Grassmann manifold
    Gr = grassmannfactory(n, p);
    problem.M = Gr;
    problem.cost = @(X)    -trace(X'*A*X);
    problem.grad = @(X)    -2*Gr.egrad2rgrad(X, A*X);
    problem.hess = @(X, H) -2*Gr.ehess2rhess(X, A*X, A*H, H);

    % Execute some checks on the derivatives for early debugging.
    % These can be commented out.
    % checkgradient(problem);
    % pause;
    % checkhessian(problem);
    % pause;

    % Issue a call to a solver. A random initial guess will be chosen and
    % default options are selected except for the ones we specify here.
    options.Delta_bar = 8*sqrt(p);
    [X, costX, info, options] = trustregions(problem, [], options); %#ok<ASGLU>

    fprintf('Options used:\n');
    disp(options);

    % For our information, Manopt can also compute the spectrum of the
    % Riemannian Hessian on the tangent space at (any) X. Computing the
    % spectrum at the solution gives us some idea of the conditioning of
    % the problem. If we were to implement a preconditioner for the
    % Hessian, this would also inform us on its performance.
    %
    % Notice that (typically) all eigenvalues of the Hessian at the
    % solution are positive, i.e., we find an isolated minimizer. If we
    % replace the Grassmann manifold by the Stiefel manifold, hence still
    % optimizing over orthonormal matrices but ignoring the invariance
    % cost(XQ) = cost(X) for all Q orthogonal, then we see
    % dim O(p) = p(p-1)/2 zero eigenvalues in the Hessian spectrum, making
    % the optimizer not isolated anymore.
    if Gr.dim() < 512
        evs = hessianspectrum(problem, X);
        stairs(sort(evs));
        title(['Eigenvalues of the Hessian of the cost function ' ...
               'at the solution']);
        xlabel('Eigenvalue number (sorted)');
        ylabel('Value of the eigenvalue');
    end

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
classdef TestComputeRecallPrecisionCurve
    %TestComputeRecallPrecisionCurve

    methods (Static)
        function test_1
            H = [0.8 -0.04 50; -0.05 0.9 50; 1e-5 1e-4 1];
            img1 = imread(fullfile(mexopencv.root(),'test','fruits.jpg'));
            img2 = cv.warpPerspective(img1, H);

            obj = cv.ORB();
            [~, feat1] = obj.detectAndCompute(img1);
            [~, feat2] = obj.detectAndCompute(img2);
            matcher = cv.DescriptorMatcher('BruteForce-Hamming');
            matches1to2 = matcher.radiusMatch(feat1, feat2, 100);

            %HACK: fake correct matches (see cv::evaluateGenericDescriptorMatcher)
            % Ideally we would look into the thresold overlap mask matrix
            correctMatches1to2Mask = cell(size(matches1to2));
            for i=1:numel(matches1to2)
                [~,ord] = sort([matches1to2{i}.distance]);
                mask = false(size(ord));
                mask(ord(1:fix(end/2))) = true;
                correctMatches1to2Mask{i} = mask;
            end

            recallPrecisionCurve = cv.computeRecallPrecisionCurve(...
                matches1to2, correctMatches1to2Mask);
            validateattributes(recallPrecisionCurve, {'numeric'}, ...
                {'2d', 'size',[NaN 2]});
            %plot(recallPrecisionCurve(:,1), recallPrecisionCurve(:,2))
        end

        function test_error_argnum
            try
                cv.computeRecallPrecisionCurve();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

%{
function nearestPointIndex = getNearestPoint(recallPrecisionCurve, l_precision)
    if (0 <= l_precision && l_precision <= 1)
        [~,nearestPointIndex] = min(abs(recallPrecisionCurve(:,1) - l_precision));
    else
        nearestPointIndex = -1;
end

function recall = getRecall(recallPrecisionCurve, l_precision)
    nearestPointIndex = getNearestPoint(recallPrecisionCurve, l_precision);
    if nearestPointIndex > 0
        recall = recallPrecisionCurve(nearestPointIndex,2);
    else
        recall = -1.0;
    end
end
%}

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [labels, conf_map] = ...
    APPtestDirectory(segDensity, vClassifier, hClassifier, ...
        imdir, imsegs, varargin)
% [labels, conf_map] = APPtestDirectory(segDensity, vClassifier,
%                                   hClassifier, imdir, imsegs, varargin)
%
% Gets the geometry for each image with superpixels given by imsegs.
%
% Input:
%   segDensity: structure giving probability of 2 sp having same label
%   vClassifier: segment classifier for ground/vert/sky
%   hClassifier: segment classifier for subclassses of vert
%   imdir: the directory of the images referred to by imsegs
%   imsegs: the superpixel structure for each image
%   varargin{1} (optional): the output directory for displaying results
% Output:
%   labels: structure containing labeling results for each image
%   conf_map: the likelihoods for each sp for each class for each image
%
% Copyright(C) Derek Hoiem, Carnegie Mellon University, 2005
% Permission granted to non-commercial enterprises for
% modification/redistribution under GNU GPL.
% Current Version: 1.0  09/30/2005

DO_PARALLEL = 0; % for running multiple parallel processes on directory

if length(varargin) > 0
    outdir = varargin{1};
end

for f = 1:length(imsegs)

    fn = imsegs(f).imname;
    bn = strtok(fn, '.');

    if ~DO_PARALLEL || ~exist('outdir') || ~exist([outdir '/' bn '.c.mat'])

        if exist('outdir') && DO_PARALLEL % to mark as being processed
            save([outdir '/' bn '.c.mat'], 'fn');
        end

        image = im2double(imread([imdir, '/', fn]));

        if size(image, 3) == 3

            disp(['processing image ' fn]);

            [labels(f), conf_map(f), maps{f}, pmaps(f)] = ...
                    APPtestImage(image, imsegs(f), vClassifier, hClassifier, segDensity);
            % for generating context
            [cimages, cnames] = APPclassifierOutput2confidenceImages(imsegs(f), conf_map(f));

            if length(varargin) > 0
                outdir = varargin{1};

                limage = APPgetLabeledImage(image, imsegs(f), labels(f).vert_labels, labels(f).vert_conf, ...
                    labels(f).horz_labels, labels(f).horz_conf);
                imwrite(limage, [outdir, '/', bn, '.l.jpg']);
                imwrite(image, [outdir, '/', fn]);

                % for generating context
                glabels = labels(f);
                gconf_map = conf_map(f);
                gmaps = maps{f};
                gpmaps = pmaps(f);
                save([outdir '/' bn '.c.mat'], 'glabels', 'gconf_map', 'cimages', 'cnames', 'gmaps', 'gpmaps');

                % make a vrml file
                if 0
                    APPwriteVrmlModel(imdir, imsegs(f), labels(f), outdir);
                end
            %else
            %    disp('warning: set up for ICCV, no output ==> no results')
            end

        end

        pause(0.05);

    end

end

for f = 1:length(imsegs)
    labels(f).imname = imsegs(f).imname;
end

author: zouchuhang

repo: LayoutNet


text:
function [V,E,H,h] = poly2VEH(poly)
  % POLY2VEH convert poly struct array to vertex, edge and hole lists
  %
  % [V,E,H,h] = poly2VEH(poly)
  % Input:
  %   poly struct array for each component of polygon boundary. Each struct
  %     contains a loop of vertices (x,y) and a flag for whether this is an
  %     inner hole
  % Output:
  %   V  #V by 2 list of polygon vertices
  %   E  #E by 2 list of polygon edge indices
  %   H  #H by 2 list of hole positions
  %
  % This has the problem that it is throwing away the information about which
  % sets of edges are inner boundaries
  %
  % Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)
  %
  % See also: png2poly
  %

  V = [];
  E = [];
  H = [];
  H
  components = size(poly,2);
  % loop over components collecting vertices, edges, and holes
  for component_index = 1:components
    % component should constitute at least one triangle
    if(size(poly(component_index).x,1) >=3)
      component_vertices = ...
        [poly(component_index).x, poly(component_index).y];
      V = [V; component_vertices];
      component_E = ...
        [ 1:size(poly(component_index).x,1) ;...
        [size(poly(component_index).x,1), ...
        1:(size(poly(component_index).x,1)-1)]]';
      E = [E ; ...
        size(E,1) + component_E(:,1) , ...
        size(E,1) + component_E(:,2) ];
      if poly(component_index).hole == 1
        % find point inside polygon to be hole marker
        poly(component_index)
        component_hole = ...
          point_inside_polygon(component_vertices);
        H = [H; component_hole];
      end
    end
  end

  if(size(H,1) >= size(poly,2))
    warning('Number of holes >= number of components. Something is wrong...');
  end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%% housekeeping
clear
close all
clc

%%

noclearall=true;

txt=which('driver_all');

lastsep=find(txt==filesep,1,'last');

folder=txt(1:lastsep-1);

w=what(folder);

w=regexp(w.m,'driver\d{2}','match');

w=[w{:}];

bigm=rise.empty();

for ii=1:numel(w)

    run(w{ii})

    bigm(ii)=mest;
end

author: jmaih


text:
function z = lambda_sum_smallest( Y, k )

% LAMBDA_SUM_SMALLEST    Sum of the k smallest eigenvalues of a symmetric matrix.
%     For square matrix X, LAMBDA_SUM_SMALLEST(X,K) is SUM_SMALLEST(EIG(X),k)
%     if X is Hermitian or symmetric and real; and +Inf otherwise.
%
%     An error results if X is not a square matrix.
%
%     Disciplined convex programming information:
%         LAMBDA_SUM_SMALLEST is concave and nonmonotonic (at least with
%         respect to elementwise comparison), so its argument must be affine.

z = - lambda_sum_largest( - Y, k );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
hdistance1 = 1393.5;
hdistance2 = 1425.2;
sourceheight=0:500;
soundspeed=350.4;
pathlength1 = sqrt(hdistance1^2 + sourceheight.^2);
traveltime1 = pathlength1 ./ soundspeed;
pathlength2 = sqrt(hdistance2^2 + sourceheight.^2);
traveltime2 = pathlength2 ./ soundspeed;
apparentC = (hdistance2 - hdistance1)./(traveltime2-traveltime1);
plot(sourceheight, apparentC)
xlabel('Source Height (m)');
ylabel('Apparent sound speed (m/sec)')

figure
theta = 180 * atan(sourceheight./hdistance1) / pi;
plot(theta, apparentC)
xlabel('Incidence angle (degrees)');
ylabel('Apparent sound speed (m/sec)')

author: geoscience-community-codes

repo: GISMO


text:
function qualityMeasure = imageQuality_sseq(sr, gt)

% Making sure that only functions in the SSEQ folder can be called
addpath(genpath('../algorithms/qualityMeasures/SSEQ/'));

path = pwd;
cd(['..', filesep, 'algorithms', filesep, 'qualityMeasures', filesep, 'SSEQ']);
qualityMeasure = - SSEQ(repmat(uint8(255*sr), 1, 1, 3));
cd(path);

% Making sure that only functions in the SSEQ folder can be called
rmpath(genpath('../algorithms/qualityMeasures/SSEQ/'));

author: thomas-koehler

repo: SupER


text:
function y = modulate2(x, type, center)
% MODULATE2	2D modulation
%
%	y = modulate2(x, type, [center])
%
% With TYPE = {'r', 'c' or 'b'} for modulate along the row, or column or
% both directions.
%
% CENTER secify the origin of modulation as floor(size(x)/2)+1+center
% (default is [0, 0])

if ~exist('center', 'var')
    center = [0, 0];
end

% Size and origin
s = size(x);
o = floor(s / 2) + 1 + center;

n1 = [1:s(1)] - o(1);
n2 = [1:s(2)] - o(2);

switch lower(type(1))
    case 'r'
	m1 = (-1) .^ n1;
	y = x .* repmat(m1', [1, s(2)]);

    case 'c'
	m2 = (-1) .^ n2;
	y = x .* repmat(m2, [s(1), 1]);

    case 'b'
	m1 = (-1) .^ n1;
	m2 = (-1) .^ n2;
	m = m1' * m2;
	y = x .* m;

    otherwise
	error('Invalid input type');
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
%compute the xcoordinate of the center in the tails coordinate system

function [data,units]=compute_xtailcentral_mm(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
xtailcentral_mm=cell(1,numlarvae);

for i=1:numlarvae
    larva=larvae(i);
    xtailcentral_mm{1,i}=trx(larva).xcentral_mm-trx(larva).xtail_mm;
end

units=parseunits('mm');
data=xtailcentral_mm;


text:
start_vy_series = {[0, 0.2, 0.4, 0.6, 0.8]
    [0, -0.2, -0.4, -0.6, -0.8]
    0
    0};
% start_vy = [0, -0.2, -0.4, -0.6, -0.8];
target_vy = [0.0];
start_vx_series = {0
    0
    [0,0.1,0.2,0.3,0.4]
    [0,-0.1,-0.2,-0.3,-0.4]};
target_vx = [0.0];
T = 0.4;
subfolder_name = 'library5';
if ~exist(fullfile('local', subfolder_name, 'transition'), 'dir')
    mkdir(fullfile('local', subfolder_name, 'transition'));
end

fit_data = load(fullfile('local', subfolder_name, 'midstep_fit.mat'));
P = fit_data.P;
dP = fit_data.dP;

fit_data_right = load(fullfile('local', subfolder_name, 'midstep_fit_right.mat'));
P_right = fit_data_right.P;
dP_right = fit_data_right.dP;
% nlp.Phase(1).removeCost('stateDeviation_RightStance');
% nlp.Phase(3).removeConstraint('periodicState_LeftStance');
start_vy = [-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8];
start_vx = [-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4];
counter = 1;
for i=1:4
    start_vy = start_vy_series{i};
    start_vx = start_vx_series{i};
    target_gait_file = fullfile('local', subfolder_name, sprintf('gait_X%0.1f_Y%.1f.mat', target_vx, target_vy));
    target_gait = load(target_gait_file);
    guess = [target_gait.gait; target_gait.gait(1:3)];
    guess(3).tspan = guess(3).tspan - guess(3).tspan(1);
    guess(7).tspan = guess(7).tspan - guess(7).tspan(1);

    for vx = start_vx
        for vy = start_vy

            data_name = fullfile('local', subfolder_name, 'transition', ...
                sprintf('gait_X%0.1f_Y%.1f_TO_X%0.1f_Y%.1f.mat', vx, vy, target_vx, target_vy));
            if exist(data_name, 'file')
                continue;
            end

            vel_lb = [min(vx,target_vx), min(vy,target_vy)];
            vel_ub = [max(vx,target_vx), max(vy,target_vy)];

            start_gait_file = fullfile('local', subfolder_name, sprintf('gait_X%0.1f_Y%.1f.mat', vx, vy));
            start_gait = load(start_gait_file);
            x0 = [start_gait.gait(1).states.x(:,11);start_gait.gait(1).states.dx(:,11)];

            target_gait_file = fullfile('local', subfolder_name, sprintf('gait_X%0.1f_Y%.1f.mat', target_vx, target_vy));
            target_gait = load(target_gait_file);
            xf = [target_gait.gait(3).states.x(:,11);target_gait.gait(3).states.dx(:,11)];

            bounds = trans_opt.GetBounds(robot, vel_lb, vel_ub, T, x0, xf);
            bounds.LeftStance1.midState_QFit = P;
            bounds.LeftStance1.midState_dQFit = dP;

            bounds.RightStance2.midState_QFit = P_right;
            bounds.RightStance2.midState_dQFit = dP_right;

            trans_opt.updateVariableBounds(nlp, bounds);
            %             nlp.Phase(3).removeConstraint('periodicState_LeftStance');
            % update desired gaits
            gait_cost = [target_gait.gait(1), target_gait.gait(3), target_gait.gait(1), target_gait.gait(3)];
            trans_opt.updateDesiredGait(nlp, system, gait_cost);

            % update initial condition

            trans_opt.updateInitCondition(nlp,guess);

            diary_name = fullfile('local', subfolder_name, 'transition', ...
                sprintf('gait_X%0.1f_Y%.1f_TO_X%0.1f_Y%.1f.txt', vx, vy, target_vx, target_vy));
            diary(diary_name);

            [gait, sol, info, total_time] = trans_opt.solve(nlp);
            pause(1);
            if info.status == 0 || info.status == 1
                data_name = fullfile('local', subfolder_name, 'transition', ...
                    sprintf('gait_X%0.1f_Y%.1f_TO_X%0.1f_Y%.1f.mat', vx, vy, target_vx, target_vy));
                fprintf('Saving gait %s\n', data_name);
            else
                data_name = fullfile('local', subfolder_name, 'transition', ...
                    sprintf('gait_X%0.1f_Y%.1f_TO_X%0.1f_Y%.1f_Failed.mat', vx, vy, target_vx, target_vy));
                fprintf('Saving (failed) gait %s\n', data_name);
            end
            pause(0.2);
            diary off;
            pause(1);
            save(data_name, 'gait', 'sol', 'info', 'bounds', 'total_time');

            counter = counter + 1;
            guess = gait;

        end

    end
end;

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function a2 = r8cmat_to_r8mat ( lda, m, n, a1 )

%*****************************************************************************80
%
%% R8CMAT_TO_R8MAT transfers data from an R8CMAT to an R8MAT.
%
%  Discussion:
%
%    An R8CMAT is an MxN array of R8's, stored with a leading dimension LD,
%    accessible as a vector:
%      (I,J) -> (I+J*LD).
%    or as a doubly-dimensioned array, if declared A(LD,N):
%      (I,J) -> A(I,J)
%
%    An R8MAT is an MxN array of R8's,
%    accessible as a vector:
%      (I,J) -> (I+J*M).
%    or as a doubly-dimensioned array, if declared A(M,N):
%      (I,J) -> A(I,J)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer LDA, the leading dimension of A1.
%    M <= LDA.
%
%    Input, integer M, the number of rows of data.
%    M <= LDA.
%
%    Input, integer N, the number of columns of data.
%
%    Input, real A1(LDA,N), the M by N matrix to be copied.
%
%    Output, real A2(M,N), a copy of the
%    information in the MxN submatrix of A1.
%
      a2(1:m,1:n) = a1(1:m,1:n);

      return
      end


text:
function [forecast_estimates]=olsforecast(data_endo_a,data_exo_p,Fperiods,betahat,Bhat,sigmahat,n,m,p,k,const,Fband)

% function [forecast_estimates]=olsforecast(data_endo_a,data_exo_p,Fperiods,betahat,Bhat,sigmahat,n,m,p,k,const,Fband)
% computes unconditional forecast values (point estimates and confidence bands) for the OLS VAR model
% inputs:  - matrix 'data_endo_a': matrix of pre-forecast endogenous data
%          - matrix 'data_exo_p': predicted values for the exogenous variables over the forecast periods
%          - integer 'Fperiods': number of forecast periods
%          - vector 'betahat': OLS VAR coefficients in vectorised form (defined in 1.1.15)
%          - matrix 'Bhat': OLS VAR coefficients, in non vectorised form (defined in 1.1.9)
%          - matrix 'sigmahat': OLS VAR variance-covariance matrix of residuals (defined in 1.1.10)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
%          - integer 'const': 0-1 value to determine if a constant is included in the model
%          - scalar 'Fband': confidence level for forecasts
% outputs: - cell 'forecast_estimates': lower bound, point estimates, and upper bound for the unconditional forecasts

% this function implements the chain rule of forecast, desibed p38 of the technical guide

% point estimates for forecasts are obtained using the chain rule for forecast in Lutkepohl (1991), equation (2.2.3) p 29
% approximate confidence interval are obtained from Lutkepohl (1991), formula (3.5.15) p 89, based on the sample mean squared error matrix (2.2.10)

% generate the matrix of predicted exogenous variables
% if the constant has been retained, augment the matrices of exogenous with a column of ones:
if const==1
data_exo_p=[ones(Fperiods,1) data_exo_p];
% if no constant was included, do nothing
else
end

% then generate the point estimates
% recover the lagged endogenous required to produce the forecasts
temp=data_endo_a(end-p+1:end,:);
% repeat the process for periods T+1 to T+h
for ii=1:Fperiods
% Define the matrix of regressors X by using lagX on temp; retain only the last row of the matrix
   % if no exogenous variable is present at all in the model (neither constant nor other exogenous), define X from the endogenous variables only
   if isempty(data_exo_p)
   X=bear.lagx(temp,p-1);
   X=X(end,:);
   % if there are exogenous vaiables, concatenate them next to the endogenous
   else
   X=bear.lagx(temp,p-1);
   X=[X(end,:) data_exo_p(ii,:)];
   end
% obtain predicted value for T+ii
yp=X*Bhat;
% concatenate the transpose of yp to the top of temp
temp=[temp;yp];
% repeat until values are obtained for T+h
end

% finally, generate the confidence bands
% this requires to estimate the forecast error matrix sigmaf for each forecast period
% to do so, it is first necessary to obtain irfs
[irfmatrix,~]=bear.irfsim(betahat,eye(n),n,m,p,k,Fperiods);
% then initiate sigmaf for period 1
sigmaf(:,:,1)=irfmatrix(:,:,1)*sigmahat*irfmatrix(:,:,1)';
% and increment for each forecast period
for ii=2:Fperiods
sigmaf(:,:,ii)=sigmaf(:,:,ii-1)+irfmatrix(:,:,ii)*sigmahat*irfmatrix(:,:,ii)';
end
% with the sigmaf series, it is possible to compute the confidence intervals
% first compute the percentile of the normal distribution corresponding to size of the confidence interval
c_low=norminv((1-Fband)/2,0,1);
c_high=norminv(Fband+(1-Fband)/2,0,1);

% finally, create and fill the forecast_estimates cell
forecast_estimates=cell(n,1);
for ii=1:n
% record forecast, point estimate
forecast_estimates{ii,1}(2,:)=temp(p+1:end,ii)';
   % then loop over forecast periods
   for jj=1:Fperiods
   % record forecast, lower bound
   forecast_estimates{ii,1}(1,jj)=forecast_estimates{ii,1}(2,jj)+c_low*sigmaf(ii,ii,jj)^0.5;
   % record forecast, upper bound
   forecast_estimates{ii,1}(3,jj)=forecast_estimates{ii,1}(2,jj)+c_high*sigmaf(ii,ii,jj)^0.5;
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function DEM_demo_factor_analysis
% Demo for Probabilistic Factor Analysis; This uses a hierarchical model
% under the constraint that the causes have a deterministic and stochastic
% components.  The aim is to recover the true subspace of the real causes.

rng('default')

% non-hierarchical linear generative model (static)
%==========================================================================
n     = 8;
m     = 2;
M     = spm_DEM_M('FA',[n m]);

% create data
%==========================================================================
N     = 8;                                        % length of data sequence
X     = randn(size(M(1).pE));
DEM   = spm_DEM_generate(M,N,{X},{4});

% Initialise parameters
%--------------------------------------------------------------------------
DEM.class = 'FA';
DEM   = spm_dem_initialise(DEM);

% DEM estimation
%==========================================================================
DEM.M(1).E.nE = 16;
DEM   = spm_DEM(DEM);

% compare real and estimated factor and causes
%==========================================================================

% plot
%--------------------------------------------------------------------------
subplot(2,2,2)
v     = DEM.qU.v{2};
u     = DEM.pU.v{2};
plot(v'*pinv(full(v'))*u')
hold on
plot(u',':')
title({'real and rotated causes','Factor analysis'},'FontSize',16)
axis square
grid on


text:
function obj = rejection(obj, mode, varargin)

% Input parsing ----------------------------------------------------------------

validMode = {'Attribute' 'dAlphaNaN'};

p = inputParser;
p.addRequired( 'mode'                , @(x) any(strcmpi(x, validMode)));
p.addParameter('AttributeName'  , '' , @ischar);          % for mode 'Attribute'
p.addParameter('AttributeMinMax', [] , @(x) numel(x)==2); % for mode 'Attribute'
p.parse(mode, varargin{:});
p = p.Results;
% Clear required inputs to avoid confusion
clear mode

% Start ------------------------------------------------------------------------

procHierarchy = {'CORRPOINTS' 'REJECTION'};
msg('S', procHierarchy);
msg('I', procHierarchy, sprintf('Corr. points label = ''%s''', obj.label));
msg('I', procHierarchy, sprintf('IN: Mode = ''%s''', p.mode));
if strcmpi(p.mode, 'Attribute')
    msg('I', procHierarchy, sprintf('IN: AttributeName = ''%s''', p.AttributeName));
    msg('V', p.AttributeMinMax(1), 'IN: AttributeMinMax(1)', 'Prec', 3);
    msg('V', p.AttributeMinMax(2), 'IN: AttributeMinMax(2)', 'Prec', 3);
end

% Correspondences present? -----------------------------------------------------

if obj.noCP == 0
    msg('I', procHierarchy, 'no correspondences present!');
    msg('E', procHierarchy);
    return;
end

% Rejection of all correspondences where isnan(dAlpha) -------------------------

if strcmpi(p.mode, 'dAlphaNaN')

    idxRejection = isnan(obj.dAlpha);

end

% Rejection based on an attribute ----------------------------------------------

if strcmpi(p.mode, 'Attribute')

    att = obj.getAttribute(p.AttributeName);

    idxRejection = att >= min(p.AttributeMinMax) & att <= max(p.AttributeMinMax);

end

% Rejection --------------------------------------------------------------------

msg('V', sum(idxRejection), 'no. of rejected correspondences', 'Prec', 0);

obj.X1  = obj.X1(~idxRejection,:);
obj.X2  = obj.X2(~idxRejection,:);

% Rejection of attributes
if ~isempty(obj.A) , att  = fields(obj.A) ; else att  = []; end
if ~isempty(obj.A1), att1 = fields(obj.A1); else att1 = []; end
if ~isempty(obj.A2), att2 = fields(obj.A2); else att2 = []; end
for a = 1:numel(att) , obj.A.(att{a})   = obj.A.(att{a})(~idxRejection);   end
for a = 1:numel(att1), obj.A1.(att1{a}) = obj.A1.(att1{a})(~idxRejection); end
for a = 1:numel(att2), obj.A2.(att2{a}) = obj.A2.(att2{a})(~idxRejection); end

if ~isempty(obj.idxPC1), obj.idxPC1 = obj.idxPC1(~idxRejection,:); end
if ~isempty(obj.idxPC2), obj.idxPC2 = obj.idxPC2(~idxRejection,:); end

% End --------------------------------------------------------------------------

msg('E', procHierarchy);
obj.info

end

author: pglira


text:
function [b b0] = circ_skewness(alpha, w, dim)

% [b b0] = circ_skewness(alpha,w,dim)
%   Calculates a measure of angular skewness.
%
%   Input:
%     alpha     sample of angles
%     [w        weightings in case of binned angle data]
%     [dim      statistic computed along this dimension, 1]
%
%     If dim argument is specified, all other optional arguments can be
%     left empty: circ_skewness(alpha, [], dim)
%
%   Output:
%     b         skewness (from Pewsey)
%     b0        alternative skewness measure (from Fisher)
%
%   References:
%     Pewsey, Metrika, 2004
%     Statistical analysis of circular data, Fisher, p. 34
%
% Circular Statistics Toolbox for Matlab

% By Philipp Berens, 2009
% berens@tuebingen.mpg.de

if nargin < 3
  dim = 1;
end

if nargin < 2 || isempty(w)
  % if no specific weighting has been specified
  % assume no binning has taken place
	w = ones(size(alpha));
else
  if size(w,2) ~= size(alpha,2) || size(w,1) ~= size(alpha,1)
    error('Input dimensions do not match');
  end
end

% compute neccessary values
R = circ_r(alpha,w,[],dim);
theta = circ_mean(alpha,w,dim);
[~, rho2 mu2] = circ_moment(alpha,w,2,true,dim);

% compute skewness
theta2 = repmat(theta, size(alpha)./size(theta));
b = sum(w.*(sin(2*(circ_dist(alpha,theta2)))),dim)./sum(w,dim);
b0 = rho2.*sin(circ_dist(mu2,2*theta))./(1-R).^(3/2);    % (formula 2.29)

author: buzsakilab


text:
function [MDP] = spm_MDP_game_KL(MDP,varargin)
% action selection using active inference (KL formulation)
% FORMAT [MDP] = spm_MDP_game_KL(MDP,[EU])
%
% EU              - optional flag to invoke expected utility only
%
% MDP.T           - process depth (the horizon)
% MDP.N           - number of variational iterations (default 4)
% MDP.S(N,1)      - true initial state
%
% MDP.A(O,N)      - Likelihood of O outcomes given N hidden states
% MDP.B{M}(N,N)   - transition probabilities among hidden states (priors)
% MDP.C(N,1)      - terminal cost probabilities (prior over hidden states)
% MDP.D(N,1)      - initial prior probabilities (prior over hidden states)
%
% MDP.V(T,P)      - P allowable policies (control sequences over T times)
%
% optional:
% MDP.s(1 x T)    - vector of true states  - for deterministic solutions
% MDP.o(1 x T)    - vector of observations - for deterministic solutions
% MDP.a(1 x T)    - vector of action       - for deterministic solutions
% MDP.w(1 x T)    - vector of precisions   - for deterministic solutions
%
% MDP.B{T,M}(N,N) - model transition probabilities for each time point
% MDP.G{T,M}(N,N) - true  transition probabilities for each time point
%                   (default: MDP.G{T,M} = MDP.G{M} = MDP.B{M})
%
% MDP.plot        - switch to suppress graphics: (default: [0])
% MDP.alpha       - upper bound on precision (Gamma hyperprior - shape [8])
% MDP.beta        - precision over precision (Gamma hyperprior - rate  [1])
%
% produces:
%
% MDP.P(M,T)   - probability of emitting an action 1,...,M at time 1,...,T
% MDP.Q(N,T)   - an array of conditional (posterior) expectations over
%                N hidden states and time 1,...,T
% MDP.O(O,T)   - a sparse matrix of ones encoding outcomes at time 1,...,T
% MDP.S(N,T)   - a sparse matrix of ones encoding states at time 1,...,T
% MDP.U(M,T)   - a sparse matrix of ones encoding action at time 1,...,T
% MDP.W(1,T)   - posterior expectations of precision
% MDP.d        - simulated dopamine responses
%
% This routine provides solutions of active inference (minimisation of
% variational free energy) using a generative model based upon a Markov
% decision process. This model and inference scheme is formulated
% in discrete space and time. This means that the generative model (and
% process) are  finite state  machines or hidden Markov models whose
% dynamics are given by transition probabilities among states and the
% likelihood corresponds to a particular outcome conditioned upon
% hidden states. For simplicity, this routine assumes that action
% and hidden controls are isomorphic. If the dynamics of transition
% probabilities of the true process are not provided, this routine will use
% the equivalent probabilities from the generative model.
%
% This particular scheme is designed for any allowable policies or control
% sequences specified in MDP.V. Constraints on allowable policies can limit
% the numerics or combinatorics considerable. For example, situations in
% which one action can be selected at one time can be reduced to T polices
% - with one (shift) control being emitted at all possible time points.
% This specification of polices simplifies the generative model, allowing a
% fairly exhaustive model of potential outcomes - eschewing a mean field
% approximation over successive control states. In brief, the agent simply
% represents the current state and states in the immediate and distant
% future.
%
% The transition probabilities are a cell array of probability transition
% matrices corresponding to each (discrete) the level of the control state.
%
% Mote that the conditional expectations are functions of time but also
% contain expectations about fictive states over time at each time point.
% To create time dependent transition probabilities, one can specify a
% function in place of the transition probabilities under different levels
% of control.
%
% Partially observed Markov decision processes can be modelled by
% specifying a likelihood (as part of a generative model) and absorbing any
% probabilistic mapping between (isomorphic) hidden states and outcomes
% into the transition probabilities G.
%
% See also: spm_MDP, which uses multiple future states and a mean field
% approximation for control states - but allows for different actions
% at all times (as in control problems).
%
% See also: spm_MDP_game, which generalises this scheme and replaces prior
% beliefs about KL control with minimisation of expected free energy.
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_MDP_game_KL.m 7679 2019-10-24 15:54:07Z spm $

% set up and preliminaries
%==========================================================================

% options and precision defaults
%--------------------------------------------------------------------------
try, PLOT  = MDP.plot;  catch, PLOT  = 0; end
try, alpha = MDP.alpha; catch, alpha = 8; end
try, beta  = MDP.beta;  catch, beta  = 1; end
try, N     = MDP.N;     catch, N     = 4; end

% set up figure if necessary
%--------------------------------------------------------------------------
if PLOT
    if ishandle(PLOT)
        figure(PLOT); clf
        PLOT = 2;
    else
        spm_figure('GetWin','MDP'); clf
    end
end

% generative model and initial states
%--------------------------------------------------------------------------
T     = MDP.T;                     % process depth (the horizon)
Ns    = size(MDP.B{1},1);          % number of hidden states
Nb    = size(MDP.B,1);             % number of time-dependent probabilities
Nu    = size(MDP.B,2);             % number of hidden controls
p0    = eps;                       % smallest probability

% likelihood model (for a partially observed MDP implicit in G)
%--------------------------------------------------------------------------
try
    A  = MDP.A + p0;
    No = size(MDP.A,1);           % number of outcomes
catch
    A  = speye(Ns,Ns) + p0;
    No = Ns;
end
A     = A*diag(1./sum(A));
lnA   = log(A);

% transition probabilities (priors)
%--------------------------------------------------------------------------
for i = 1:T
    for j = 1:Nu
        if i == 1 || Nb == T
            B{i,j}   = MDP.B{i,j} + p0;
            B{i,j}   = B{i,j}*diag(1./sum(B{i,j}));
        else
            B{i,j}   = B{1,j};
        end
    end
end

% terminal probabilities (priors)
%--------------------------------------------------------------------------
try
    C = spm_vec(MDP.C) + p0;
catch
    C = ones(Ns,1);
end
C     = C/sum(C);
lnC   = log(C);

% intital probabilities (priors)
%--------------------------------------------------------------------------
try
    D = spm_vec(MDP.D) + p0;
catch
    D = ones(Ns,1);
end
D     = D/sum(D);
lnD   = log(D);

% generative process (assume the true process is the same as the model)
%--------------------------------------------------------------------------
try
    G = MDP.G;
catch
    G = MDP.B;
end
Ng    = size(G,1);
for i = 1:T
    for j = 1:Nu
        if i == 1 || Ng == T
            G{i,j} = G{i,j} + p0;
            G{i,j} = G{i,j}*diag(1./sum(G{i,j}));
        else
            G{i,j} = G{1,j};
        end
    end
end

% policies and their expectations
%--------------------------------------------------------------------------
V      = MDP.V;
u      = zeros(size(V,2),1);
Np     = size(V,2);                % number of allowable policies
w      = 1:Np;                     % indices of allowable policies

% initial states and outcomes
%--------------------------------------------------------------------------
[p q]  = max(A*MDP.S(:,1));        % initial outcome (index)
s      = find( MDP.S(:,1));        % initial state   (index)
o      = sparse(1,1,q,1,T);        % observations    (index)
S      = sparse(s,1,1,Ns,T);       % states sampled  (1 in K vector)
O      = sparse(q,1,1,No,T);       % states observed (1 in K vector)
a      = sparse(1, T);             % action (index)
U      = sparse(Nu,T);             % action selected (1 in K vector)
P      = sparse(Nu,T);             % posterior beliefs about control
E      = sparse(T,Np);             % posterior beliefs about policies
W      = sparse(1,T);              % posterior precision

% sufficient statistics of hidden states (past, current and last)
%--------------------------------------------------------------------------
gamma  = [];                       % simulated dopamine responses
x      = zeros(Ns,T);

% solve
%==========================================================================
for t  = 1:T

    % allowable policies
    %----------------------------------------------------------------------
    if t > 1

        % record posterior expectations over policies
        %------------------------------------------------------------------
        E(t - 1,w) = u;

        % retain allowable policies (that are consistent with last action)
        %------------------------------------------------------------------
        j = ismember(V(t - 1,:),a(t - 1));
        V = V(:,j);
        u = u(j);
        w = w(j);

    end

    % conditional KL divergence (under allowable policies)
    %======================================================================
    Np    = size(V,2);                % number of allowable policies
    Q     = zeros(Np,Ns);             % value of policies x current state
    for k = 1:Np

        % compositon of future states
        %------------------------------------------------------------------
        Bj = 1;
        for j = t:T
            Bj = B{j,V(j,k)}*Bj;
        end

        % divergence or information gain
        %------------------------------------------------------------------
        if nargin > 1
            Q(k,:) = lnC'*Bj;
        else
            Q(k,:) = lnC'*Bj - sum(Bj.*log(Bj));
        end

    end

    % Variational iterations (assuming precise inference about past action)
    %======================================================================
    for i  = 1:N

        % present state (x)
        %------------------------------------------------------------------
        if t == 1
            v  = lnD;
        else
            v  = log(B{t - 1,a(t - 1)}*x(:,t - 1));
        end
        v      = v + lnA(o(t),:)' + W(t)*Q'*u;
        x(:,t) = spm_softmax(v);

        % precision (W)
        %------------------------------------------------------------------
        if isfield(MDP,'w')
            W(t) = MDP.w(t);
        else
            v    = beta - u'*Q*x(:,t);
            W(t) = alpha/v;
        end

        % policy (u)
        %------------------------------------------------------------------
        v      = W(t)*Q*x(:,t);
        u      = spm_softmax(v);
        E(t,w) = u;

        % re-compute precision for first iteration
        %------------------------------------------------------------------
        if t == 1
            v    = beta - u'*Q*x(:,t);
            W(t) = alpha/v;
        end

        % simulated dopamine responses (precision as each iteration)
        %------------------------------------------------------------------
        gamma(end + 1,1) = W(t);

    end

    % posterior expectations (control)
    %======================================================================
    for j = 1:Nu
        for k = t:T
            P(j,k) = sum(u(ismember(V(k,:),j)));
        end
    end

    % next action (the action that minimises expected free energy)
    %------------------------------------------------------------------
    try
        a(t) = MDP.a(t);
    catch
        try
            a(t) = find(rand < cumsum(P(:,t)),1);
        catch
            error('there are no more allowable policies')
        end
    end

    % save action
    %------------------------------------------------------------------
    U(a(t),t) = 1;

    % sampling of next state (outcome)
    %======================================================================
    if t < T

        % next sampled state
        %------------------------------------------------------------------
        try
            s(t + 1) = MDP.s(t + 1);
        catch
            s(t + 1) = find(rand < cumsum(G{t,a(t)}(:,s(t))),1);
        end

        % next obsverved state
        %------------------------------------------------------------------
        try
            o(t + 1) = MDP.o(t + 1);
        catch
            o(t + 1) = find(rand < cumsum(A(:,s(t + 1))),1);
        end

        % save outcome and state sampled
        %------------------------------------------------------------------
        W(1,t + 1)        = W(t);
        O(o(t + 1),t + 1) = 1;
        S(s(t + 1),t + 1) = 1;

    end

    % plot
    %======================================================================
    if PLOT > 0

        % posterior beliefs about hidden states
        %------------------------------------------------------------------
        subplot(4,2,1)
        imagesc(1 - [x C*max(max(x))/max(C)])
        if size(x,1) > 128
            hold on, spm_spy(x,16,1), hold off
        end
        title('Inferred states (and utility)','FontSize',14)
        xlabel('Time','FontSize',12)
        ylabel('Hidden state','FontSize',12)

        % posterior beliefs about control states
        %==================================================================
        subplot(4,2,2)

        % make previous plots dotted lines
        %------------------------------------------------------------------
        if T > 2
            h     = get(gca,'Children'); hold on
            for i = 1:length(h)
                set(h(i),'LineStyle',':');
            end
            plot(P')
            title('Inferred policy','FontSize',14)
            xlabel('Time','FontSize',12)
            ylabel('Control state','FontSize',12)
            spm_axis tight
        else
            bar(P)
            title('Inferred policy','FontSize',14)
            xlabel('Contol state','FontSize',12)
            ylabel('Posterior expectation','FontSize',12)
        end

        % policies
        %------------------------------------------------------------------
        subplot(4,2,3)
        imagesc(MDP.V')
        title('Allowable policies','FontSize',14)
        ylabel('Policy','FontSize',12)
        xlabel('Time','FontSize',12)

        % expectations over policies
        %------------------------------------------------------------------
        subplot(4,2,4)
        imagesc(E')
        title('Posterior probability','FontSize',14)
        ylabel('Policy','FontSize',12)
        xlabel('Time','FontSize',12)

        % true state (outcome)
        %------------------------------------------------------------------
        subplot(4,2,5)
        if size(S,1) > 128
            spm_spy(S,16)
        else
            imagesc(1 - S)
        end
        title('True states','FontSize',14)
        ylabel('State','FontSize',12)

        % sample (observation)
        %------------------------------------------------------------------
        subplot(4,2,7)
        if size(O,1) > 128
            spm_spy(O,16,1)
        else
            imagesc(1 - O)
        end
        title('Observed states','FontSize',14)
        xlabel('Time','FontSize',12)
        ylabel('State','FontSize',12)

        % action sampled (selected)
        %------------------------------------------------------------------
        subplot(4,2,6)
        if size(U,1) > 128
            spm_spy(U,16,1)
        else
            imagesc(1 - U)
        end
        title('Selected action','FontSize',14)
        ylabel('Action','FontSize',12)

        % expected action
        %------------------------------------------------------------------
        subplot(4,2,8)
        plot((1:length(gamma))/N,gamma)
        title('Expected precision (confidence)','FontSize',14)
        xlabel('Time','FontSize',12)
        ylabel('Precision','FontSize',12)
        spm_axis tight
        drawnow

    end

end

% deconvolve to simulate dopamine responses
%--------------------------------------------------------------------------
da     = pinv( tril(toeplitz(exp(-((1:length(gamma)) - 1)'/8))) )*gamma;

% assemble results and place in NDP structure
%--------------------------------------------------------------------------
MDP.P  = P;              % probability of action at time 1,...,T - 1
MDP.Q  = x;              % conditional expectations over N hidden states
MDP.O  = O;              % a sparse matrix, encoding outcomes at 1,...,T
MDP.S  = S;              % a sparse matrix, encoding the states
MDP.U  = U;              % a sparse matrix, encoding the action
MDP.W  = W;              % posterior expectations of precision
MDP.d  = gamma;          % simulated dopamine responses
MDP.da = da;             % simulated dopamine responses (deconvolved)

author: spm


text:
function [S,f,R,Serr]=mtspectrumtrigpt(data,E,win,params,fscorr)
% Multi-taper time-frequency spectrum - point process times
%
% Usage:
%
% [S,f,R,Serr]=mtspectrumtrigpt(data,E,win,params,fscorr)
% Input:
%       data        (structure array of one channel of spike times;
%                   also accepts 1d column vector of spike times) -- required
%       E           (event times) - required
%       win         (in the form [winl winr] i.e window around each event)--
%                                                 required
%       params: structure with fields tapers, pad, Fs, fpass, err, trialave
%       - optional
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                   (1) A numeric vector [TW K] where TW is the
%                       time-bandwidth product and K is the number of
%                       tapers to be used (less than or equal to
%                       2TW-1).
%                   (2) A numeric vector [W T p] where W is the
%                       bandwidth, T is the duration of the data and p
%                       is an integer such that 2TW-p tapers are used. In
%                       this form there is no default i.e. to specify
%                       the bandwidth, you have to specify T and p as
%                       well. Note that the units of W and T have to be
%                       consistent: if W is in Hz, T must be in seconds
%                       and vice versa. Note that these units must also
%                       be consistent with the units of params.Fs: W can
%                       be in Hz if and only if params.Fs is in Hz.
%                       The default is to use form 1 with TW=3 and K=5
%
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%           err  (error calculation [1 p] - Theoretical error bars; [2 p] - Jackknife error bars
%                                   [0 p] or 0 - no error bars) - optional. Default 0.
%           trialave (average over events when 1, don't average when 0) -
%           default 0
%       fscorr   (finite size corrections, 0 (don't use finite size corrections) or
%                1 (use finite size corrections) - optional
%                (available only for spikes). Defaults 0.
%
% Output:
%       S       (triggered spectrum in form frequency x events if trialave=0;
%               function of frequency if trialave=1)
%       f       (frequencies)
%       R       (spike rate)
%       Serr    (error bars) - only for err(1)>=1

if nargin < 3; error('Need data, events and window parameters'); end;
if nargin < 2; params=[]; end;
[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear tapers pad Fs fpass trialave
if nargin < 5 || isempty(fscorr); fscorr=0; end;
if nargout > 3 && err(1)==0; error('Cannot compute errors if err(1)=0'); end;
data=createdatamatpt(data,E,win);
if nargout==4;
    [S,f,R,Serr]=mtspectrumpt(data,params,fscorr);
else
    [S,f,R]=mtspectrumpt(data,params,fscorr);
end;


text:
function polygon = setLMpolygon(x,y)

polygon.x = x(:);
polygon.y = y(:);

%for i = 1:length(x)
%    polygon.pt(i).x = num2str(x(i));
%    polygon.pt(i).y = num2str(y(i));
%end

author: CSAILVision

repo: LabelMeToolbox


text:
function [pathway, notMapped]=mapPathwayRxnNames(pathway,originalLabels,newLabels)
% mapPathwayRxnNames
%   For mapping labels in the pathway object. Useful if you want to change
%   what is shown in the reaction boxes.
%
%   pathway         pathway structure representing the pathway to be drawn
%   originalLabels  cell array with the original reaction labels
%   newLabels       cell array with the new reaction labels
%
%   pathway         an updated pathway structure
%   notMapped       a cell array with labels that could not be found in the
%                   pathway object
%
%   Usage: [pathway, notMapped]=mapPathwayRxnNames(pathway,originalLabels,newLabels)

if numel(originalLabels)~=numel(newLabels)
    EM='The new label cell array must have the same length as the old label cell array';
    dispEM(EM);
end

mapped=false(numel(originalLabels),1);

for i=1:numel(pathway.listOfSpecies)
    if strcmpi(pathway.listOfSpecies(i).type,'PROTEIN')
        I=find(ismember(originalLabels,pathway.listOfSpecies(i).name));
        if any(I)
            if numel(I)==1
                pathway.listOfSpecies(i).name=newLabels{I};
                mapped(I)=true;
            else
                EM=['The label "' pathway.listOfSpecies(i).name '" was found in several positions in oldLabels'];
                dispEM(EM);
            end
        end
    end
end

notMapped=originalLabels(~mapped);
end

author: SysBioChalmers

repo: RAVEN