text:
classdef ShapeCircle < ShapeGeneral
    %ShapeCircle represents a circular geographical selection of events
    %
    % see also ShapeGeneral, ShapePolygon

    properties (SetObservable = true, AbortSet=true)
        Radius (1,1) double = 5 % active radius in units defined by the RefEllipsoid
    end

    methods
        function obj=ShapeCircle(varargin)
            % SHAPECIRCLE create a circular shape
            %
            % ShapeCircle() :
            % ShapeCircle('dlg') create via a dialg box
            %
            % CIRCLE: select using circle with a defined radius. define with 2 clicks or mouseover and press "R"

            % UNASSIGNED: clear shape

            obj@ShapeGeneral();
            if ~isempty(varargin)
               dosomething
            end

            report_this_filefun();

            %axes(findobj(gcf,'Tag','mainmap_ax')); % should be the map, with lon/lat
            obj.Type='circle';
            try
                ra=ShapeGeneral.ShapeStash.Radius;
            catch
                ra=obj.Radius;
            end
            obj.AllowVertexEditing = false;
            addlistener(obj, 'Radius', 'PostSet', @obj.notifyShapeChange);
            if numel(varargin)==0
                do_nothing;
            elseif strcmpi(varargin{1},'dlg')
                stashedshape = ShapeGeneral.ShapeStash;
                sdlg.prompt = ['Choose Radius [',obj.RefEllipsoid.LengthUnit,']:'];
                sdlg.value = ra;
                sdlg(2).prompt = 'Center X :'; sdlg(2).value=stashedshape.X0;
                sdlg(3).prompt = 'Center Y :'; sdlg(3).value=stashedshape.Y0;
                [~,cancelled,obj.Radius,obj.Points(1),obj.Points(2)]=smart_inputdlg('Define Circle',sdlg);
                if cancelled
                    beep
                    disp('Circle creation cancelled by user')
                    return
                end
            else
                oo=ShapeCircle.selectUsingMouse(gca, obj.RefEllipsoid);
                if ~isempty(oo)
                    obj=oo;
                else
                    return
                end
            end
        end

        function val=Outline(obj, col)
            if iscartesian(obj.RefEllipsoid)
                pts = exp(1i*pi*linspace(0,2*pi,3600)') .* obj.Radius;
                x = real(pts)+ obj.X0;
                y = imag(pts) + obj.Y0;
                val = [x,y];
            else
                [lat,lon]=reckon(obj.Y0,obj.X0,obj.Radius,(0:.1:360)',obj.RefEllipsoid);
                val=[lon, lat];
            end
            if exist('col','var')
                val=val(:,col);
            end
        end

        function moveTo(obj, x, y)
            if isnan(obj.Points)
                obj.Points=[0 0];
            end
            moveTo@ShapeGeneral(obj,x,y)
        end

        function s=toStruct(obj)
            s=toStruct@ShapeGeneral(obj);
            s.RadiusKm = obj.Radius;
        end

        function s = toStr(obj)
            cardinalDirs='SNWE';
            isN=obj.Y0>=0; NS=cardinalDirs(isN+1);

            isE=obj.X0>=0; EW=cardinalDirs(isE+3);
            s = sprintf('Circle with R:%s %s, centered at ( %s %s, %s %s)',...
                num2str(obj.Radius),...
                obj.RefEllipsoid.LengthUnit,...
                num2str(abs(obj.Y0)), NS,...
                num2str(abs(obj.X0)), EW);
        end

        function summary(obj)
            helpdlg(obj.toStr,'Circle');
        end

        function add_shape_specific_context(obj,c)
            uimenu(c,'label','Choose Radius','MenuSelectedFcn',@chooseRadius)
            uimenu(c,'label','Snap To N Events','MenuSelectedFcn',@snapToEvents)

            function snapToEvents(~,~)
                ZG=ZmapGlobal.Data;
                nc=inputdlg('Number of events to enclose','Edit Circle',1,{num2str(ZG.ni)});
                nc=round(str2double(nc{1}));
                if ~isempty(nc) && ~isnan(nc)
                    ZG.ni=nc;
                    zmw=get(ancestor(c,'figure'),'UserData');
                    if isa(zmw,'ZmapWindow')
                        ca=zmw.rawcatalog;
                    else
                        ca=ZG.primeCatalog;
                    end

                    [~,obj.Radius]=ca.selectClosestEvents(obj.Y0, obj.X0, [],nc,'DistanceOnly');
                    obj.Radius=obj.Radius;%+0.005;
                end
            end

            function chooseRadius(~,~)
                radiusInputText = ['Choose Radius [',obj.RefEllipsoid.LengthUnit,']'];
                nc=inputdlg(radiusInputText,'Edit Circle',1,{num2str(obj.Radius)});
                nc=str2double(nc{1});
                if ~isempty(nc) && ~isnan(nc)
                    obj.Radius=nc;
                end

            end

        end

        function [mask]=isinterior(obj,otherX, otherY, include_boundary)
            % isinterior true if value is within this circle's radius of center. Radius inclusive.
            %
            % overridden because using polygon approximation is too inaccurate for circles
            %
            % [mask]=obj.isinterior(otherX, otherY)

            if ~exist('include_boundary','var')
                include_boundary = true;
            end
            if isempty(obj.Points)||isnan(obj.Points(1))
                mask = ones(size(otherX));
            else
                otherX(ismissing(otherY))= missing;
                otherY(ismissing(otherX))= missing;
                % return a vector of size otherX that is true where item is inside polygon

                if iscartesian(obj.RefEllipsoid)
                    dists = sqrt((otherY-obj.Y0).^2 + (otherX-obj.X0).^2);
                else
                    dists = distance(obj.Y0, obj.X0, otherY, otherX, obj.RefEllipsoid);
                end

                if ~include_boundary
                    mask = dists < obj.Radius;
                else
                    mask = dists <= obj.Radius;
                end
            end
        end

        function finishedMoving(obj, movedObject, deltas)
            centerX = mean(bounds2(movedObject.XData));
            centerY = mean(bounds2(movedObject.YData));

            obj.Radius=obj.Radius.* abs(deltas(3)); % NO NEGATIVE RADII
            obj.Points=[centerX,centerY];
        end

        function save(obj, filelocation, delimiter)
            persistent savepath
            if ~exist('filelocation','var') || isempty(filelocation)
                if isempty(savepath)
                    savepath = ZmapGlobal.Data.Directories.data;
                end
                [filename,pathname,filteridx]=uiputfile(...
                    {'*.mat','MAT-files (*.mat)';...
                    '*.csv;*.txt;*.dat','ASCII files (*.csv, *.txt, *.dat)'},...
                    'Save Circle',...
                    fullfile(savepath,'zmap_shape.mat'));
                if filteridx==0
                    msg.dbdisp('user cancelled shape save');
                    return
                end
                filelocation=fullfile(pathname,filename);
            end
            [savepath,~,ext] = fileparts(filelocation);
            if ext==".mat"
                zmap_shape = obj; %#ok<NASGU>
                save(filelocation,'zmap_shape');
            else
                if ~exist('delimiter','var'), delimiter = ',';end
                radiusName = ['Radius[',shortenLengthUnit(obj.RefEllipsoid.LengthUnit),']'];
                tb=table(obj.X0, obj.Y0,obj.Radius,'VariableNames',{'Latitude','Longitude',radiusName});
                writetable(tb,filelocation,'Delimiter',delimiter);
            end

        end
    end

    methods(Static)

        function obj=selectUsingMouse(ax, ref_ellipsoid)

            [ss,ok] = selectSegmentUsingMouse(ax,'r', @circ_update);
            delete(findobj(gca,'Tag','tmp_circle_outline'));
            if ~ok
                obj=[];
                return
            end
            obj=ShapeCircle();
            obj.Points=ss.xy1;
            obj.Radius=ss.dist;

            function circ_update(stxy, ~, d)
                h=findobj(gca,'Tag','tmp_circle_outline');
                if isempty(h)
                    h=line(nan,nan,'Color','r','DisplayName','Rough Outline','LineWidth',2,'Tag','tmp_circle_outline');
                end
                if iscartesian(ref_ellipsoid)
                    pts = exp(1i*pi*linspace(0,2*pi,120)') .* d;
                    h.XData = real(pts)+ stxy(1);
                    h.YData = imag(pts) + stxy(2);
                else
                    [lat,lon]=reckon(stxy(2),stxy(1),d,(0:3:360)',ref_ellipsoid);
                    h.XData=lon;
                    h.YData=lat;
                end
            end
        end
    end

end


text:
function y=rfft(x,n,d)
%RFFT     Calculate the DFT of real data Y=(X,N,D)
% Data is truncated/padded to length N if specified.
%   N even:	(N+2)/2 points are returned with
% 			the first and last being real
%   N odd:	(N+1)/2 points are returned with the
% 			first being real
% In all cases fix(1+N/2) points are returned
% D is the dimension along which to do the DFT

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: rfft.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s=size(x);
if prod(s)==1
    y=x
else
    if nargin <3 || isempty(d)
        d=find(s>1,1);
        if nargin<2
            n=s(d);
        end
    end
    if isempty(n)
        n=s(d);
    end
    y=fft(x,n,d);
    y=reshape(y,prod(s(1:d-1)),n,prod(s(d+1:end)));
    s(d)=1+fix(n/2);
    y(:,s(d)+1:end,:)=[];
    y=reshape(y,s);
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
% goplotpso.m
% default plotting script used in PSO functions
%
% this script is not a function,
% it is a plugin for the main PSO routine (pso_Trelea_vectorized)
% so it shares all the same variables, be careful with variable names
% when making your own plugin

% Brian Birge
% Rev 2.0
% 3/1/06

% setup figure, change this for your own machine

 clf
 set(gcf,'Position',[651    31   626   474]); % this is the computer dependent part
 %set(gcf,'Position',[743    33   853   492]);
 set(gcf,'Doublebuffer','on');

% particle plot, upper right
 subplot('position',[.7,.6,.27,.32]);
 set(gcf,'color','k')

 plot3(pos(:,1),pos(:,D),out,'b.','Markersize',7)

 hold on
 plot3(pbest(:,1),pbest(:,D),pbestval,'g.','Markersize',7);
 plot3(gbest(1),gbest(D),gbestval,'r.','Markersize',25);

 % crosshairs
 offx = max(abs(min(min(pbest(:,1)),min(pos(:,1)))),...
            abs(max(max(pbest(:,1)),max(pos(:,1)))));

 offy = max(abs(min(min(pbest(:,D)),min(pos(:,D)))),...
            abs(min(max(pbest(:,D)),max(pos(:,D)))));
 plot3([gbest(1)-offx;gbest(1)+offx],...
       [gbest(D);gbest(D)],...
       [gbestval;gbestval],...
       'r-.');
 plot3([gbest(1);gbest(1)],...
       [gbest(D)-offy;gbest(D)+offy],...
       [gbestval;gbestval],...
       'r-.');

 hold off

 xlabel('Dimension 1','color','y')
 ylabel(['Dimension ',num2str(D)],'color','y')
 zlabel('Cost','color','y')

 title('Particle Dynamics','color','w','fontweight','bold')

 set(gca,'Xcolor','y')
 set(gca,'Ycolor','y')
 set(gca,'Zcolor','y')
 set(gca,'color','k')

 % camera control
 view(2)
 try
   axis([gbest(1)-offx,gbest(1)+offx,gbest(D)-offy,gbest(D)+offy]);
 catch
   axis([VR(1,1),VR(1,2),VR(D,1),VR(D,2)]);
 end

% error plot, left side
 subplot('position',[0.1,0.1,.475,.825]);
  semilogy(tr(find(~isnan(tr))),'color','m','linewidth',2)
  %plot(tr(find(~isnan(tr))),'color','m','linewidth',2)
  xlabel('epoch','color','y')
  ylabel('gbest val.','color','y')

  if D==1
     titstr1=sprintf(['%11.6g = %s( [ %9.6g ] )'],...
                gbestval,strrep(functname,'_','\_'),gbest(1));
  elseif D==2
     titstr1=sprintf(['%11.6g = %s( [ %9.6g, %9.6g ] )'],...
                gbestval,strrep(functname,'_','\_'),gbest(1),gbest(2));
  elseif D==3
     titstr1=sprintf(['%11.6g = %s( [ %9.6g, %9.6g, %9.6g ] )'],...
                gbestval,strrep(functname,'_','\_'),gbest(1),gbest(2),gbest(3));
  else
     titstr1=sprintf(['%11.6g = %s( [ %g inputs ] )'],...
                gbestval,strrep(functname,'_','\_'),D);
  end
  title(titstr1,'color','m','fontweight','bold');

  grid on
%  axis tight

  set(gca,'Xcolor','y')
  set(gca,'Ycolor','y')
  set(gca,'Zcolor','y')
  set(gca,'color','k')

  set(gca,'YMinorGrid','off')

% text box in lower right
% doing it this way so I can format each line any way I want
subplot('position',[.62,.1,.29,.4]);
  clear titstr
  if trelea==0
       PSOtype  = 'Common PSO';
       xtraname = 'Inertia Weight : ';
       xtraval  = num2str(iwt(length(iwt)));

     elseif trelea==2 | trelea==1

       PSOtype  = (['Trelea Type ',num2str(trelea)]);
       xtraname = ' ';
       xtraval  = ' ';

     elseif trelea==3
       PSOtype  = (['Clerc Type 1"']);
       xtraname = '\chi value : ';
       xtraval  = num2str(chi);

  end
  if isnan(errgoal)
    errgoalstr='Unconstrained';
  else
    errgoalstr=num2str(errgoal);
  end
  if minmax==1
     minmaxstr = ['Maximize to : '];
  elseif minmax==0
     minmaxstr = ['Minimize to : '];
  else
     minmaxstr = ['Target   to : '];
  end

  if rstflg==1
     rststat1 = 'Environment Change';
     rststat2 = ' ';
  else
     rststat1 = ' ';
     rststat2 = ' ';
  end

  titstr={'PSO Model: '      ,PSOtype;...
          'Dimensions : '    ,num2str(D);...
          '# of particles : ',num2str(ps);...
          minmaxstr          ,errgoalstr;...
          'Function : '      ,strrep(functname,'_','\_');...
          xtraname           ,xtraval;...
          rststat1           ,rststat2};

  text(.1,1,[titstr{1,1},titstr{1,2}],'color','g','fontweight','bold');
  hold on
  text(.1,.9,[titstr{2,1},titstr{2,2}],'color','m');
  text(.1,.8,[titstr{3,1},titstr{3,2}],'color','m');
  text(.1,.7,[titstr{4,1}],'color','w');
  text(.55,.7,[titstr{4,2}],'color','m');
  text(.1,.6,[titstr{5,1},titstr{5,2}],'color','m');
  text(.1,.5,[titstr{6,1},titstr{6,2}],'color','w','fontweight','bold');
  text(.1,.4,[titstr{7,1},titstr{7,2}],'color','r','fontweight','bold');

  % if we are training a neural net, show a few more parameters
  if strcmp('pso_neteval',functname)
    % net is passed from trainpso to pso_Trelea_vectorized in case you are
    % wondering where that structure comes from
    hiddlyrstr = [];
    for lyrcnt=1:length(net.layers)
       TF{lyrcnt} = net.layers{lyrcnt}.transferFcn;
       Sn(lyrcnt) = net.layers{lyrcnt}.dimensions;
       hiddlyrstr = [hiddlyrstr,', ',TF{lyrcnt}];
    end
    hiddlyrstr = hiddlyrstr(3:end);

    text(0.1,.35,['#neur/lyr = [ ',num2str(net.inputs{1}.size),'  ',...
               num2str(Sn),' ]'],'color','c','fontweight','normal',...
               'fontsize',10);
    text(0.1,.275,['Lyr Fcn: ',hiddlyrstr],...
       'color','c','fontweight','normal','fontsize',9);

  end

  legstr = {'Green = Personal Bests';...
            'Blue  = Current Positions';...
            'Red   = Global Best'};
  text(.1,0.025,legstr{1},'color','g');
  text(.1,-.05,legstr{2},'color','b');
  text(.1,-.125,legstr{3},'color','r');

  hold off

  set(gca,'color','k');
  set(gca,'visible','off');

  drawnow

author: HuangCongQing


text:
function [ccls cfun] = get_data_class()
% get class supported by ann
d = annmex();
ccls = class(d);

cfun = @(x) cast(x, ccls);

author: CHLNDDEV


text:
function [hrv, R_t, R_amp, R_index, S_t, S_amp]  = rpeakdetect(data,samp_freq,thresh,testmode)

% [hrv, R_t, R_amp, R_index, S_t, S_amp]  = rpeakdetect(data, samp_freq, thresh, testmode);
% R_t == RR points in time, R_amp == amplitude
% of R peak in bpf data & S_amp == amplitude of
% following minmum. sampling frequency (samp_freq = 256Hz
% by default) only needed if no time vector is
% specified (assumed to be 1st column or row).
% The 'triggering' threshold 'thresh' for the peaks in the 'integrated'
% waveform is 0.2 by default.  testmode = 0 (default) indicates
% no graphics diagnostics. Otherwise, you get to scan through each segment.
%
% A batch QRS detector based upon that of Pan, Hamilton and Tompkins:
% J. Pan \& W. Tompkins - A real-time QRS detection algorithm
% IEEE Transactions on Biomedical Engineering, vol. BME-32 NO. 3. 1985.
% P. Hamilton \& W. Tompkins. Quantitative Investigation of QRS
% Detection  Rules Using the MIT/BIH Arrythmia Database.
% IEEE Transactions on Biomedical Engineering, vol. BME-33, NO. 12.1986.
%
% Similar results reported by the authors above were achieved, without
% having to tune the thresholds on the MIT DB. An online version in C
% has also been written.
%
% Written by G. Clifford gari@ieee.org and made available under the
% GNU general public license. If you have not received a copy of this
% license, please download a copy from http://www.gnu.org/
%
% Please distribute (and modify) freely, commenting
% where you have added modifications.
% The author would appreciate correspondence regarding
% corrections, modifications, improvements etc.
%
% gari@ieee.org

%%%%%%%%%%% make threshold default 0.2 -> this was 0.15 on MIT data
if nargin < 4
   testmode = 0;
end
%%%%%%%%%%% make threshold default 0.2 -> this was 0.15 on MIT data
if nargin < 3
   thresh = 0.2;
end
%%%%%%%%%%% make sample frequency default 256 Hz
if nargin < 2
   samp_freq = 256;
   if(testmode==1)
       fprintf('Assuming sampling frequency of %iHz\n',samp_freq);
   end
end

%%%%%%%%%%% check format of data %%%%%%%%%%
[a b] = size(data);
if(a>b)
 len =a;
end
if(b>a)
 len =b;
end

%%%%%%%%%% if there's no time axis - make one
if (a | b == 1);
% make time axis
  tt = 1/samp_freq:1/samp_freq:ceil(len/samp_freq);
  t = tt(1:len);
  x = data;
end
%%%%%%%%%% check if data is in columns or rows
if (a == 2)
  x=data(:,1);
  t=data(:,2);
end
if (b == 2)
  t=data(:,1);
  x=data(:,2);
end

%%%%%%%%% bandpass filter data - assume 256hz data %%%%%
 % remove mean
 x = x-mean(x);

 % FIR filtering stage
 bpf=x; %Initialise
if( (samp_freq==128) & (exist('filterECG128Hz')~=0) )
        bpf = filterECG128Hz(x);
end
if( (samp_freq==256) & (exist('filterECG256Hz')~=0) )
        bpf = filterECG256Hz(x);
end

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

%%%%%%%%% differentiate data %%%%%%%%%%%%%%%%%%%%%%%%%%%
 dff = diff(bpf);  % now it's one datum shorter than before

%%%%%%%%% square data    %%%%%%%%%%%%%%%%%%%%%%%%%%%
 sqr = dff.*dff;   %
 len = len-1; % how long is the new vector now?

%%%%%%%%% integrate data over window 'd' %%%%%%%%%%%%%%%%%%%%%%%%%
 d=[1 1 1 1 1 1 1]; % window size - intialise
 if (samp_freq>=256) % adapt for higher sampling rates
   d = [ones(1,round(7*samp_freq/256))];
 end
 % integrate
 mdfint = medfilt1(filter(d,1,sqr),10);
 % remove filter delay for scanning back through ECG
 delay = ceil(length(d)/2);
 mdfint = mdfint(delay:length(mdfint));
%%%%%%%%% segment search area %%%%%%%%%%%%%%%%%%%%%%%
 %%%% first find the highest bumps in the data %%%%%%
 max_h = max (mdfint(round(len/4):round(3*len/4)));

 %%%% then build an array of segments to look in %%%%%
 %thresh = 0.2;
 poss_reg = mdfint>(thresh*max_h);

%%%%%%%%% and find peaks %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%% find indices into boudaries of each segment %%%
 left  = find(diff([0 poss_reg'])==1); % remember to zero pad at start
 right = find(diff([poss_reg' 0])==-1); % remember to zero pad at end

 %%%% loop through all possibilities
 for(i=1:length(left))
    [maxval(i) maxloc(i)] = max( bpf(left(i):right(i)) );
    [minval(i) minloc(i)] = min( bpf(left(i):right(i)) );
    maxloc(i) = maxloc(i)-1+left(i); % add offset of present location
    minloc(i) = minloc(i)-1+left(i); % add offset of present location
 end

 R_index = maxloc;
 R_t   = t(maxloc);
 R_amp = maxval;
 S_amp = minval;   %%%% Assuming the S-wave is the lowest
                   %%%% amp in the given window
 S_t   = t(minloc);

%%%%%%%%%% check for lead inversion %%%%%%%%%%%%%%%%%%%
 % i.e. do minima precede maxima?
 if (minloc(length(minloc))<maxloc(length(minloc)))
  R_t   = t(minloc);
  R_amp = minval;
  S_t   = t(maxloc);
  S_amp = maxval;
 end

%%%%%%%%%%%%
hrv  = diff(R_t);
resp = R_amp-S_amp;

%%%%%%%%%%%%%%%%%%%%
if (testmode~=0)
figure(1)
hold off
subplot(4,1,1)
plot(t,x);hold on;plot(t,bpf,'r')
title('raw ECG (blue) and zero-pahse FIR filtered ECG (red)')
ylabel('ECG')
hold off;
subplot(4,1,2)
plot(t(1:length(mdfint)),mdfint);hold on;
%plot(t(1:length(sqr)),sqr);hold on;
plot(t,max(mdfint)*bpf/(2*max(bpf)),'r')
plot(t(left),mdfint(left),'og')
plot(t(right),mdfint(right),'om')
title('Integrated data with scan boundaries over scaled ECG')
ylabel('Int ECG')
hold off;
subplot(4,1,3)
plot(t,bpf,'r');hold on;
plot(R_t,R_amp,'+k');
plot(S_t,S_amp,'+g');
title('ECG with R-peaks (black) and S-points (green) over ECG')
ylabel('ECG+R+S')
hold off;
subplot(4,1,4)
hold off
plot(R_t(1:length(hrv)),hrv,'r+')
hold on
title('RR intervals')
ylabel('RR (s)')
hold off
fprintf('Press any key for next block of data\n');
pause
end

author: peterhcharlton

repo: RRest


text:
function Sf = FreeBoundary(S,t,V,K,type)

Sf = zeros(1,length(t));
eps_star = K*1e-5;

switch type
    case 'put'
        for j = 1:length(t)
            Sf(j) = S(find(abs(V(:,j)-K+S)< eps_star, 1, 'last'));
        end
    case 'call'
        for j = 1:length(t)
            Sf(j) = S(find(abs(V(:,j)+K-S)< eps_star, 1, 'first'));
        end
end

author: Sable

repo: mcbench-benchmarks


text:
function pass = test_vectorizeOp()
%TEST_VECTORIZEOP    Test that the CHEBOP vectorization method works as expected

%% Setup
% Domain and some of CHEBFUNs
d = [1, 3];
x = chebfun(@(x) x, d);
f = exp(sin(4*x));
g = cos(x.^2).^2;
h = 1./(5+x.^2);
u = sin(5*x.^2);
v = exp(x/2);
% Some parameters
b = 4*pi;
c = -5*exp(1);

%% Simple case, already vectorized
fun = @(u) diff(u, 2) + b.*u.^3;
% Automatically vectorized
vecFun = chebop.vectorizeOp(fun);
% The answer we expect (FUN above is already vectorized, so they're identical)
exactFun = fun;
% Check that we got the correct string and the correct function
pass(1,1) = strcmp(func2str(exactFun), func2str(vecFun));
pass(1,2) = norm(exactFun(u) - vecFun(u)) == 0;

%% Simple case, already vectorized, with both x and u
fun = @(x,u) diff(u, 2) + b.*u.^3;
% Automatically vectorized
vecFun = chebop.vectorizeOp(fun);
% The answer we expect (FUN above is already vectorized, so they're identical)
exactFun = fun;
% Check that we got the correct string and the correct function
pass(2,1) = strcmp(func2str(exactFun), func2str(vecFun));
pass(2,2) = norm(exactFun(x,u) - vecFun(x,u)) == 0;

%% Vectorization needed, scalar case
fun = @(u) diff(u, 2) + b*u^3;
% Automatically vectorized
vecFun = chebop.vectorizeOp(fun);
% The answer we expect
exactFun = @(u) diff(u, 2) + b.*u.^3;
% Check that we got the correct string and the correct function
pass(3,1) = strcmp(func2str(exactFun), func2str(vecFun));
pass(3,2) = norm(exactFun(u) - vecFun(u)) == 0;

%% Vectorization needed, scalar case with both x and u
fun = @(x,u) diff(u, 2) + b*u^3;
% Automatically vectorized
vecFun = chebop.vectorizeOp(fun);
% The answer we expect
exactFun = @(x,u) diff(u, 2) + b.*u.^3;
% Check that we got the correct string and the correct function
pass(4,1) = strcmp(func2str(exactFun), func2str(vecFun));
pass(4,2) = norm(exactFun(x,u) - vecFun(x,u)) == 0;

%% Vectorization needed, scalar case with a CHEBFUN
fun = @(u) diff(u, 2) + f*u^3;
% Automatically vectorized
vecFun = chebop.vectorizeOp(fun);
% The answer we expect
exactFun = @(u) diff(u, 2) + f.*u.^3;
% Check that we got the correct string and the correct function
pass(5,1) = strcmp(func2str(exactFun), func2str(vecFun));
pass(5,2) = norm(exactFun(u) - vecFun(u)) == 0;

%% Vectorization needed, scalar case with a CHEBFUN, both x and u
fun = @(x,u) diff(u, 2) + f*u^3;
% Automatically vectorized
vecFun = chebop.vectorizeOp(fun);
% The answer we expect
exactFun = @(x,u) diff(u, 2) + f.*u.^3;
% Check that we got the correct string and the correct function
pass(6,1) = strcmp(func2str(exactFun), func2str(vecFun));
pass(6,2) = norm(exactFun(x,u) - vecFun(x,u)) == 0;

%% System case, already vectorized
fun = @(x,u,v) [diff(u, 2) + b.*u.^3+v.*u./h; diff(v,2) + v.^2./c + g.*u./v];
% Automatically vectorized
vecFun = chebop.vectorizeOp(fun);
% The answer we expect (FUN above is already vectorized, so they're identical)
exactFun = fun;
% Check that we got the correct string and the correct function
pass(7,1) = strcmp(func2str(exactFun), func2str(vecFun));
pass(7,2) = norm(exactFun(x,u,v) - vecFun(x,u,v)) == 0;

%% System case, partially vectorized
fun = @(x,u,v) [diff(u, 2) + b*u.^3+v.*u/h; diff(v,2) + v.^2./c + g*u./v];
% Automatically vectorized
vecFun = chebop.vectorizeOp(fun);
% The answer we expect
exactFun = @(x,u,v) [diff(u, 2) + b.*u.^3+v.*u./h; ...
    diff(v,2) + v.^2./c + g.*u./v];
% Check that we got the correct string
pass(8,1) = strcmp(func2str(exactFun), func2str(vecFun));
pass(8,2) = norm(exactFun(x,u,v) - vecFun(x,u,v)) == 0;

%% System case, no vectorization
fun = @(x,u,v) [diff(u, 2) + b*u^3+v*u/h; diff(v,2) + v^2/c + g*u/v];
% Automatically vectorized
vecFun = chebop.vectorizeOp(fun);
% The answer we expect
exactFun = @(x,u,v) [diff(u, 2) + b.*u.^3+v.*u./h; ...
    diff(v,2) + v.^2./c + g.*u./v];
% Check that we got the correct string and the correct function
pass(9,1) = strcmp(func2str(exactFun), func2str(vecFun));
pass(9,2) = norm(exactFun(x,u,v) - vecFun(x,u,v)) == 0;

%% For a function handle, will get the same results
s = @sin;
sVec = chebop.vectorizeOp(s);
% Check that we got the correct string and the correct function
pass(10,1) = strcmp(func2str(s), func2str(sVec));
pass(10,2) = norm(s(u)-sVec(u)) == 0;
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function obj = viewSubjectMovieFrames(obj)
	% Creates frame snippets from movies of a subject, used to get a sense of how well cross-session alignment should work.
	% Biafra Ahanonu
	% branched from controllerAnalysis
	% started: 2014.08.01 [16:09:16] (2015.09.14 14:40)
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.06.18 [21:41:07] - added modelVarsFromFilesCheck() to check and load signals if user hasn't already.
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
		% 2022.03.16 [08:45:40] - Update code standards.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% fileFilterRegexp = 'concat';
	display(repmat('#',1,21))
	display('computing signal peaks...')
	nFiles = length(obj.rawSignals);

	[fileIdxArray idNumIdxArray nFilesToAnalyze nFiles] = obj.getAnalysisSubsetsToAnalyze();
	subjectList = unique(obj.subjectStr(fileIdxArray));

	movieSettings = inputdlg({...
			'regexp for movies',...
			'group by (1) subject, (2) folder, (3) all montage, or (4) cell extraction max proj',...
			'video save directory',...
			'normalize movie',...
			'start:end frames (do not leave blank)'...
		},...
		'movie frames settings',1,...
		{...
			'concat',...
			'1',...
			obj.videoSaveDir,...
			'0',...
			'100:101'...
		}...
	);
	setNo = 1;
	fileFilterRegexp = movieSettings{setNo};setNo = setNo+1;
	usrIdxChoiceType = str2num(movieSettings{setNo});setNo = setNo+1;
	obj.videoSaveDir = movieSettings{setNo};setNo = setNo+1;
	normalizeMovieSwitch = str2num(movieSettings{setNo});setNo = setNo+1;
	frameList = str2num(movieSettings{setNo});setNo = setNo+1;

	% Miji;
	% MIJ.start;
	manageMiji('startStop','start');

	switch usrIdxChoiceType
		case 1
			for thisSubjectStr=subjectList
				validFoldersIdx = find(strcmp(thisSubjectStr,obj.subjectStr));
				% validManualIdx = find(arrayfun(@(x) isempty(x{1}),obj.validManual));
				% classifyFoldersIdx = intersect(validFoldersIdx,validManualIdx);
				movieList = getFileList({obj.inputFolders{validFoldersIdx}}, fileFilterRegexp);
				movieList
				subjectMovieFrames = loadMovieList(movieList,'convertToDouble',0,'frameList',frameList);
				% movieFrame = squeeze(movieFrame(:,:,1));
				subjectMovieFrames = subjectMovieFrames(:,:,1:2:end);
				if normalizeMovieSwitch==1
					[subjectMovieFrames] = normalizeVector(single(subjectMovieFrames),'normRange','zeroToOne');
					[subjectMovieFrames] = normalizeMovie(subjectMovieFrames,'normalizationType','meanSubtraction');
				end

				if ~isempty(obj.videoSaveDir)
					movieSavePathBase = strcat(obj.videoSaveDir,filesep,'subjectMovieFrames');
					if (~exist(movieSavePathBase,'dir')) mkdir(movieSavePathBase); end;
					movieSavePath = strcat(movieSavePathBase,filesep,thisSubjectStr{1},'.h5');
					[output] = writeHDF5Data(subjectMovieFrames,movieSavePath);
					movieSavePath = strcat(movieSavePathBase,filesep,thisSubjectStr{1},'.tiff');
					options.comp = 'no';
					options.overwrite = true;
					movieSavePath
					saveastiff(subjectMovieFrames, movieSavePath, options);
				end

				MIJ.createImage(thisSubjectStr{1}, subjectMovieFrames, true);
				% pause
			end
			uiwait(ciapkg.overloaded.msgbox('press OK to finish','Success','modal'));
			% MIJ.run('Close All Without Saving');
			manageMiji('startStop','closeAllWindows');
		case 2
			[fileIdxArray idNumIdxArray nFilesToAnalyze nFiles] = obj.getAnalysisSubsetsToAnalyze();
			for thisFileNumIdx = 1:nFilesToAnalyze
				try
					thisFileNum = fileIdxArray(thisFileNumIdx);
					obj.fileNum = thisFileNum;
					display(repmat('=',1,21))
					% display([num2str(thisFileNum) '/' num2str(length(fileIdxArray)) ': ' obj.fileIDNameArray{obj.fileNum}]);
					dispStr = [num2str(thisFileNumIdx) '/' num2str(nFilesToAnalyze) ' (' num2str(thisFileNum) '/' num2str(nFiles) '): ' obj.folderBaseSaveStr{obj.fileNum}];
					display(dispStr);
					% validManualIdx = find(arrayfun(@(x) isempty(x{1}),obj.validManual));
					% classifyFoldersIdx = intersect(validFoldersIdx,validManualIdx);
					movieList = getFileList({obj.inputFolders{obj.fileNum}}, fileFilterRegexp);
					cellfun(@display,movieList);
					subjectMovieFrames = loadMovieList(movieList,'convertToDouble',0,'frameList',frameList);
					% movieFrame = squeeze(movieFrame(:,:,1));
					subjectMovieFrames = subjectMovieFrames(:,:,1:2:end);

					if normalizeMovieSwitch==1
						[subjectMovieFrames] = normalizeVector(single(subjectMovieFrames),'normRange','zeroToOne');
						[subjectMovieFrames] = normalizeMovie(subjectMovieFrames,'normalizationType','meanSubtraction');
					end

					if ~isempty(obj.videoSaveDir)
						movieSavePathBase = strcat(obj.videoSaveDir,filesep,'subjectMovieFrames');
						if (~exist(movieSavePathBase,'dir')) mkdir(movieSavePathBase); end;
						movieSavePath = strcat(movieSavePathBase,filesep,obj.folderBaseSaveStr{obj.fileNum},'.h5');
						[output] = writeHDF5Data(subjectMovieFrames,movieSavePath);
						movieSavePath = strcat(movieSavePathBase,filesep,obj.folderBaseSaveStr{obj.fileNum},'.tiff');
						options.comp = 'no';
						display(['saving: ' movieSavePath])
						saveastiff(subjectMovieFrames, movieSavePath, options);
					end
					dispStr = [' (' num2str(thisFileNum) '/' num2str(nFiles) '): ' obj.folderBaseSaveStr{obj.fileNum}];
					MIJ.createImage(dispStr, subjectMovieFrames, true);
					for foobar=1:3; MIJ.run('In [+]'); end
				catch err
					display(repmat('@',1,7))
					disp(getReport(err,'extended','hyperlinks','on'));
					display(repmat('@',1,7))
				end
				% pause
			end
			uiwait(ciapkg.overloaded.msgbox('press OK to finish','Success','modal'));
			% MIJ.run('Close All Without Saving');
			manageMiji('startStop','closeAllWindows');
		case 3
			[fileIdxArray idNumIdxArray nFilesToAnalyze nFiles] = obj.getAnalysisSubsetsToAnalyze();
			movieListAll = {};
			primaryMovie = {};
			for thisFileNumIdx = 1:nFilesToAnalyze
				thisFileNum = fileIdxArray(thisFileNumIdx);
				obj.fileNum = thisFileNum;
				display(repmat('=',1,21))
				identifyingText{thisFileNumIdx} = obj.subjectStr{obj.fileNum};
				movieList = getFileList({obj.inputFolders{obj.fileNum}}, fileFilterRegexp);
				primaryMovie{thisFileNumIdx} = single(loadMovieList(movieList,'convertToDouble',0,'frameList',frameList));

				% movieListAll = {movieListAll{:} movieList{:}};
			end
			% subjectMovieFrames
			[primaryMovie] = createMontageMovie(primaryMovie,'identifyingText',identifyingText,'normalizeMovies', zeros([length(primaryMovie) 1]),'singleRowMontage',0);
			if ~isempty(obj.videoSaveDir)
				movieSavePathBase = strcat(obj.videoSaveDir,filesep,'subjectMovieFrames');
				if (~exist(movieSavePathBase,'dir')) mkdir(movieSavePathBase); end;
				movieSavePath = strcat(movieSavePathBase,filesep,obj.folderBaseSaveStr{obj.fileNum},'.h5');
				[output] = writeHDF5Data(primaryMovie,movieSavePath);
				movieSavePath = strcat(movieSavePathBase,filesep,obj.folderBaseSaveStr{obj.fileNum},'.tiff');
				options.comp = 'no';
				display(['saving: ' movieSavePath])
				saveastiff(primaryMovie, movieSavePath, options);
			end
			disp(identifyingText)
			MIJ.createImage('Montage', primaryMovie, true);
			% for foobar=1:3; MIJ.run('In [+]'); end
			uiwait(ciapkg.overloaded.msgbox('press OK to finish','Success','modal'));
			% MIJ.run('Close All Without Saving');
			manageMiji('startStop','closeAllWindows');
		case 4
			[fileIdxArray, idNumIdxArray, nFilesToAnalyze, nFiles] = obj.getAnalysisSubsetsToAnalyze();
			for thisSubjectStr = subjectList
				try
					disp(repmat('=',1,7))
					fprintf('Subject %s', thisSubjectStr{1});
					validFoldersIdx = find(strcmp(thisSubjectStr,obj.subjectStr));
					validFoldersIdx = intersect(fileIdxArray,validFoldersIdx);
					if isempty(validFoldersIdx)
						disp('Skipping...')
						continue;
					end
					subjectMovieFrames = single([]);
					for folderNo = 1:length(validFoldersIdx)
						disp('===')
						thisFileNum = validFoldersIdx(folderNo);
						obj.fileNum = thisFileNum;

						% Check that signal extraction information is loaded.
						obj.modelVarsFromFilesCheck(thisFileNum);

						[inputSignals, inputImages, signalPeaks, signalPeakIdx, valid] = modelGetSignalsImages(obj,'returnType','raw');
						if isempty(inputSignals);disp('no input signals');continue;end

						inputImages = thresholdImages(inputImages,'binary',0,'getBoundaryIndex',0,'threshold',0.4,'imageFilter','none');

						goodImages = inputImages(:,:,logical(valid));
						goodImages = max(goodImages,[],3,'omitnan');

						badImages = inputImages(:,:,~logical(valid));
						badImages = max(badImages,[],3,'omitnan');
						disp(num2str([sum(valid) size(goodImages) NaN size(inputImages)]))
						goodImages = goodImages + 0.5*badImages;

						% [goodImages] = viewAddTextToMovie(goodImages,obj.assay{obj.fileNum},12);
						try
							size(subjectMovieFrames)
							size(goodImages)
							subjectMovieFrames(:,:,folderNo) = goodImages;
						catch err
							subjectMovieFrames(:,:,folderNo) = NaN;
							display(repmat('@',1,7))
							disp(getReport(err,'extended','hyperlinks','on'));
							display(repmat('@',1,7))
						end

						try
							[subjectMovieFrames(:,:,folderNo)] = viewAddTextToMovie(subjectMovieFrames(:,:,folderNo),obj.assay{obj.fileNum},12);
						catch err
							display(repmat('@',1,7))
							disp(getReport(err,'extended','hyperlinks','on'));
							display(repmat('@',1,7))
						end
					end

					% subjectMovieFrames

					% validManualIdx = find(arrayfun(@(x) isempty(x{1}),obj.validManual));
					% classifyFoldersIdx = intersect(validFoldersIdx,validManualIdx);
					% movieList = getFileList({obj.inputFolders{validFoldersIdx}}, fileFilterRegexp);
					% movieList
					% subjectMovieFrames = loadMovieList(movieList,'convertToDouble',0,'frameList',frameList);
					% movieFrame = squeeze(movieFrame(:,:,1));
					% subjectMovieFrames = subjectMovieFrames(:,:,1:2:end);
					% if normalizeMovieSwitch==1
					% 	[subjectMovieFrames] = normalizeVector(single(subjectMovieFrames),'normRange','zeroToOne');
					% 	[subjectMovieFrames] = normalizeMovie(subjectMovieFrames,'normalizationType','meanSubtraction');
					% end

					if ~isempty(obj.videoSaveDir)
						movieSavePathBase = strcat(obj.videoSaveDir,filesep,'subjectMovieFrames');
						if (~exist(movieSavePathBase,'dir')) mkdir(movieSavePathBase); end;

						movieSavePath = strcat(movieSavePathBase,filesep,thisSubjectStr{1},'.h5');
						[success] = saveMatrixToFile(subjectMovieFrames,movieSavePath);

						movieSavePath = strcat(movieSavePathBase,filesep,thisSubjectStr{1},'.tiff');
						[success] = saveMatrixToFile(subjectMovieFrames,movieSavePath);
					end

					MIJ.createImage(thisSubjectStr{1}, subjectMovieFrames, true);
					% pause
				catch err
					display(repmat('@',1,7))
					disp(getReport(err,'extended','hyperlinks','on'));
					display(repmat('@',1,7))
				end
			end
			uiwait(ciapkg.overloaded.msgbox('press OK to finish','Success','modal'));
			% MIJ.run('Close All Without Saving');
			manageMiji('startStop','closeAllWindows');
		otherwise
			% body
	end
	% MIJ.exit;
	manageMiji('startStop','exit');
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [ h,mycolormap ] = PlotVorticity( STATE,varargin )

mycolormap = VorticityColormap( 64 );
        h=0;

if nargout>1
    return;
end
% cavity grid and PDE operators

N = sqrt(size(STATE,1))-1;

[Grid]=CollocationGrid_q(N);    % the computational grid
[~,Grid.wc] = clencurt(N);Grid.W = kron(Grid.wc,Grid.wc);   % integration weights
 disp('grid created')

[ Operators ] = CreateOperators_psi( Grid.D,eye(N+1));
disp('operators for psi created')

% grid for the plot
[xx2,yy2] = meshgrid(Grid.x,Grid.x);
[xxx,yyy] = meshgrid(-1:.005:1,-1:.005:1);
% vorticity
w = - Operators.del2*STATE;
ww = reshape(real(w),N+1,N+1);
www = interp2(xx2,yy2,ww,xxx,yyy,'spline');

a = 2;
if ~isempty(varargin)
    a = varargin{1};
end

www(www>a)=a;www(www<-a)=-a;

%         figure(100+randi(100))
        contourf(xxx,yyy,www,100,'LineStyle','None') ;    % vorticty
        axis square;
        box on;
        colormap(mycolormap);
        caxis([-a,a]);
        ax = gca;
        ax.XTick = [-1  1];
        ax.YTick = [-1  1];
end

function [ cmap ] = VorticityColormap( ncol )
%VORTICITYCOLORMAP Summary of this function goes here
%   Detailed explanation goes here
cmap=jet(ncol);
for i=ncol/8+1:7*ncol/16
    cmap(i,:)=[(i-ncol/8)/(5*ncol/16),(i-ncol/8)/(5*ncol/16),1];    % blue fading to white
end

    cmap(7*ncol/16+1:9*ncol/16,:)=1;  % white

for i=9*ncol/16+1:7*ncol/8
    cmap(i,:)=[1,1-(i-9*ncol/16)/(5*ncol/16),1-(i-9*ncol/16)/(5*ncol/16)];    % white turning red
end

end

author: arbabiha


text:
classdef SetPowerSinkActiveStateAction < AbstractEventAction
    %SetPowerSinkActiveStateAction Summary of this class goes here
    %   Detailed explanation goes here

    properties
        sink AbstractLaunchVehicleElectricalPowerSrcSnk
        activeStateToSet(1,1) logical = false;
    end

    properties(Constant)
        emptyVarArr = AbstractOptimizationVariable.empty(0,1);
    end

    methods
        function obj = SetPowerSinkActiveStateAction(sink, activeStateToSet)
            if(nargin > 0)
                obj.sink = sink;
                obj.activeStateToSet = activeStateToSet;
            end

            obj.id = rand();
        end

        function newStateLogEntry = executeAction(obj, stateLogEntry)
            newStateLogEntry = stateLogEntry;
            stageStates = newStateLogEntry.stageStates;

            pwrSinkStates = AbstractLaunchVehicleElectricalPowerSnkState.empty(1,0);
            for(i=1:length(stageStates)) %#ok<*NO4LP>
                pwrSinkStates = horzcat(pwrSinkStates, stageStates(i).powerSinkStates); %#ok<AGROW>
            end

            pwrSinks = AbstractLaunchVehicleElectricalPowerSrcSnk.empty(1,0);
            for(i=1:length(pwrSinkStates))
                pwrSinks = horzcat(pwrSinks, pwrSinkStates(i).getEpsSinkComponent()); %#ok<AGROW>
            end

            pwrSinkInd = find(pwrSinks == obj.sink,1,'first');
            pwrSinkState = pwrSinkStates(pwrSinkInd);

            pwrSinkState.setActiveState(obj.activeStateToSet);
        end

        function initAction(obj, initialStateLogEntry)
            %nothing
        end

        function name = getName(obj)
            if(obj.activeStateToSet)
                tf = 'Active';
            else
                tf = 'Inactive';
            end

            name = sprintf('Set Power Sink State (%s = %s)',obj.sink.getName(),tf);
        end

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

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

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

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

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

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

        function tf = usesTankToTankConn(obj, tankToTank)
            tf = false;
        end

        function tf = usesPwrSink(obj, powerSink)
            tf = [obj.sink] == powerSink;
        end

        function tf = usesPwrSrc(obj, powerSrc)
            tf = false;
        end

        function tf = usesPwrStorage(obj, powerStorage)
            tf = false;
        end

        function [tf, vars] = hasActiveOptimVar(obj)
            tf = false;
            vars = obj.emptyVarArr;
        end
    end

    methods(Static)
        function addActionTf = openEditActionUI(action, lv)
            [~, powerSinks] = lv.getPowerSinksListBoxStr();

            if(not(isempty(powerSinks)))
%                 addActionTf = lvd_EditActionSetPwrSinkStateGUI(action, lv);

                output = AppDesignerGUIOutput({false});
                lvd_EditActionSetPwrSinkStateGUI_App(action, lv, output);
                addActionTf = output.output{1};
            else
                addActionTf = false;
                warndlg('There are no power sinks on the vehicle.  Create one first.','Cannot Create Action','modal');
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function y = LSQRwrapper(x,transp,A,At)
% y = LSQRwrapper(x,trans,A,At)
% NESTA Version 1.1

% Stephen Becker, 11/12/09

if strcmpi(transp,'transp')
    y = At(x);
elseif strcmpi(transp,'notransp')
    y = A(x);
end

repo: REST


text:
function [ Y_hat, W_hat] = MCWNNM_ADMM1_Estimation( NL_mat, Sigma_arr, CurPat, Par )

Y_hat = zeros(size(CurPat));
W_hat    = zeros(size(CurPat));
for  i      =  1 : length(Par.SelfIndex) % For each keypatch group
    Y    =   CurPat(:, NL_mat(1:Par.nlsp,i)); % Non-local similar patches to the keypatch
    mY  =   repmat(mean( Y, 2 ),1,Par.nlsp);
    Y    =   Y-mY;
    X 	=   MCWNNM_ADMM1( Y, Sigma_arr(:, Par.SelfIndex(i)), Par); % WNNM Estimation
    Y_hat(:,NL_mat(1:Par.nlsp,i))  = Y_hat(:,NL_mat(1:Par.nlsp,i))+X+mY;
    W_hat(:,NL_mat(1:Par.nlsp,i))     = W_hat(:,NL_mat(1:Par.nlsp,i))+ones(Par.ps2ch, Par.nlsp);
end
end

author: csjunxu


text:
function testModel(downsample,model,trainNum,warpsize,cdim,ndim,seed,RoIRoD)

switch model
    case 'vgg16'
        Epoches = 40;
    case 'resnet50'
        Epoches = 20;
end

modelDir = fullfile('data','GAIC',model,RoIRoD,[num2str(trainNum) '_down' num2str(downsample)...
           '_warp' num2str(warpsize) 'cdim' num2str(cdim) 'ndim' num2str(ndim) 'seed' num2str(seed)]);

load(['imdb_GAIC',num2str(trainNum),'.mat'],'imdb');
gt_scores = cat(2,imdb.bbox.gt_scores{imdb.images.set==1});
gt_scores_means = mean(gt_scores);
gt_scores_stds = std(gt_scores);

for epoch = 1:Epoches
    acc5 = zeros(200,4);
    acc10 = zeros(200,4);
    fprintf('processing epoch %d\n',epoch);
    netStruct = load(fullfile(modelDir,['net-epoch-',num2str(epoch) '.mat']),'net');
    net = dagnn.DagNN.loadobj(netStruct.net) ;
    net.mode = 'test' ;
    net.move('gpu') ;
    probVarI = net.getVarIndex('predcls');
    net.vars(probVarI).precious = 1;
    minScale = net.meta.minScale;

    testID = imdb.images.testSet;
    tic;
    for i = 1:numel(testID)
        img = imread(imdb.meta.img_path{testID(i)});
        [x1,x2,x3] = size(img);
        boxes = imdb.bbox.boxes{testID(i)};
        ss = (boxes(3,:)-boxes(1,:)).*(boxes(4,:)-boxes(2,:));

        imre = single(imresize(img,minScale/min([x1,x2]),'bilinear'));
        imre = bsxfun(@minus,imre,net.meta.normalization.averageImage);

        [r1,r2,r3] = size(imre);
        r1 = 32*round(r1/32);
        r2 = 32*round(r2/32);
        imre = imresize(imre,[r1,r2],'bilinear');

        scale1 = r1/x1;
        scale2 = r2/x2;
        boxes_s = [];
        boxes_s(1,:) = max(floor(boxes(1,:) * scale1),1);
        boxes_s(2,:) = max(floor(boxes(2,:) * scale2),1);
        boxes_s(3,:) = min(ceil(boxes(3,:) * scale1),r1);
        boxes_s(4,:) = min(ceil(boxes(4,:) * scale2),r2);

        inputs = {'input', gpuArray(imre), 'rois', gpuArray(single([ones(1,size(boxes_s,2));boxes_s]))} ;

        net.eval(inputs) ;
        preds{i} = squeeze(gather(net.vars(probVarI).value)) ;
        preds{i} = preds{i} * gt_scores_stds + gt_scores_means;
        gts = imdb.bbox.gt_scores{testID(i)};

        srcc(i) = corr(preds{i}, gts', 'type', 'Spearman');

        [gts_sorted,id_gts] = sort(gts,'descend');
        [~,id_preds] = sort(preds{i},'descend');
        [~,id_baseline] = sort(ss,'descend');

        for k = 1:4
            for j = 1:k
                if gts(id_preds(j)) >= gts_sorted(5)
                    acc5(i,k) = acc5(i,k) + 1;
                end
            end
            acc5(i,k) = acc5(i,k) / k;
        end

        for k = 1:4
            for j = 1:k
                if gts(id_preds(j)) >= gts_sorted(10)
                    acc10(i,k) = acc10(i,k) + 1;
                end
            end
            acc10(i,k) = acc10(i,k) / k;
        end
    end
    toc;
    Acc5(epoch,:) = sum(acc5)/numel(testID);
    Acc10(epoch,:) = sum(acc10)/numel(testID);
    SRCC(epoch) = mean(srcc);
end
AllRes = cat(2,Acc5,mean(Acc5,2),Acc10,mean(Acc10,2),SRCC');
save(fullfile(modelDir,'Result.mat'),'Acc5','Acc10','SRCC','AllRes');

end

author: HuiZeng

repo: Grid-Anchor-based-Image-Cropping


text:
%% ------------------------------------------------------------------------------
% Configuration File for Visual Odometry Algorithm
%% -------------------------------------------------------------------------------

% Path to the directories containing images
data_params.path1 = '../data/kitti/00/image_0/';
data_params.path2 = '../data/kitti/00/image_1/';

% Path to calibration text file
data_params.calib_file = '../data/kitti/00/calib.txt';

% Path to groundtruth poses. Set flag to 1 to plot groundtruth as well
data_params.gt_file = '../data/kitti/poses/00.txt';
data_params.show_gt_flag = 1;

% Use parallel threads (requires Parallel Processing Toolbox)
% !! TO-DO: fix parfor and for loops for this functionality!
data_params.use_multithreads = 1;                % 0: disabled, 1: enabled

%% Read the calibration file to find parameters of the cameras
% !! TO-DO: Read from the calib_file instead

% calibration parameters for sequence 2010_03_09_drive_0000
cam_params.fx = 7.188560000000e+02;               % focal length (u-coordinate) in pixels
cam_params.cx = 6.071928000000e+02;               % principal point (u-coordinate) in pixels
cam_params.fy = 7.188560000000e+02;               % focal length (v-coordinate) in pixels
cam_params.cy = 1.852157000000e+02;               % principal point (v-coordinate) in pixels
cam_params.base = 3.861448000000e+02;             % baseline in meters (absolute value)

%% Parameters for Feature Extraction
vo_params.feature.nms_n = 8;                      % non-max-suppression: min. distance between maxima (in pixels)
vo_params.feature.nms_tau = 50;                   % non-max-suppression: interest point peakiness threshold
vo_params.feature.margin = 21;                    % leaving margin for safety while computing features ( >= 25)

%% Parameters for Feature Matching
vo_params.matcher.match_binsize = 50;             % matching bin width/height (affects efficiency only)
vo_params.matcher.match_radius = 200;             % matching radius (du/dv in pixels)
vo_params.matcher.match_disp_tolerance = 1;       % dx tolerance for stereo matches (in pixels)
vo_params.matcher.match_ncc_window = 21;          % window size of the patch for normalized cross-correlation
vo_params.matcher.match_ncc_tolerance = 0.3;      % threshold for normalized cross-correlation
% !! TO-DO: add subpixel-refinement using parabolic fitting
vo_params.matcher.refinement = 2;                 % refinement (0=none,1=pixel,2=subpixel)

%% Paramters for Feature Selection using bucketing
vo_params.bucketing.max_features = 1;             % maximal number of features per bucket
vo_params.bucketing.bucket_width = 50;            % width of bucket
vo_params.bucketing.bucket_height = 50;           % height of bucket
% !! TO-DO: add feature selection based on feature tracking
vo_params.bucketing.age_threshold = 10;           % age threshold while feature selection

%% Paramters for motion estimation
% !! TO-DO: use Nister's algorithm for Rotation estimation (along with SLERP) and
% estimate translation using weighted optimization equation
vo_params.estim.ransac_iters = 200;              % number of RANSAC iterations
vo_params.estim.inlier_threshold = 2.0;          % fundamental matrix inlier threshold
vo_params.estim.reweighing = 1;                  % lower border weights (more robust to calibration errors)

repo: Stereo-Odometry-SOFT


text:
function [g1, g2] = lfmaXlfmaKernGradient(lfmKern1, lfmKern2, t1, t2, covGrad, meanVector)

% LFMAXLFMAKERNGRADIENT Compute a cross gradient between a LFMA and a LFMA.
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel
% between two lfm kernels, both lfm kernels correspond to the accelerations.
% with the multiple output kernel.
% ARG lfmKern1 : the kernel structure associated with the first LFM
% kernel (acceleration).
% ARG lfmKern2 : the kernel structure associated with the second LFM
% kernel (acceleration).
% ARG t1 : inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see lfmKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see lfmKernExtractParam.
%
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel
% between two lfm kernels, both lfm kernels correspond to the accelerations.
% with the multiple output kernel.
% ARG lfmKern1 : the kernel structure associated with the first LFM
% kernel (acceleration).
% ARG lfmKern2 : the kernel structure associated with the second LFM
% kernel (acceleration).
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see lfmKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see lfmKernExtractParam.
%
% FORMAT
% DESC computes cross gradient of parameters of a cross kernel
% between two lfm kernels. Both kernels correspond to accelerations.
% It is supposed to be used together with the multiple output kernel.
% ARG lfmKern1 : the kernel structure associated with the first LFM
% kernel (acceleration).
% ARG lfmKern2 : the kernel structure associated with the second LFM
% kernel (acceleration).
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% ARG meanVec : precomputed factor that is used for the switching dynamical
% latent force model.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see lfmKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see lfmKernExtractParam.
%
% SEEALSO : multiKernParamInit, multiKernCompute, lfmvKernParamInit
%
% COPYRIGHT : Mauricio Alvarez, 2010

% KERN

subComponent = false; % This is just a flag that indicates if this kernel is part of a bigger kernel (SDLFM)

if nargin == 4
    covGrad = t2;
    t2 = t1;
elseif nargin == 6
    subComponent = true;
    if numel(meanVector)>1
        if size(meanVector,1) == 1,
            if size(meanVector, 2)~=size(covGrad, 2)
                error('The dimensions of meanVector don''t correspond to the dimensions of covGrad')
            end
        else
            if size((meanVector'), 2)~=size(covGrad,2)
                error('The dimensions of meanVector don''t correspond to the dimensions of covGrad')
            end
        end
    else
        if numel(t1)==1 && numel(t2)>1
            % matGrad will be row vector and so should be covGrad
            dimcovGrad = length(covGrad);
            covGrad = reshape(covGrad, [1 dimcovGrad]);
        elseif numel(t1)>1 && numel(t2)==1
            % matGrad will be column vector and sp should be covGrad
            dimcovGrad = length(covGrad);
            covGrad = reshape(covGrad, [dimcovGrad 1]);
        end
    end
end

if size(t1, 2) > 1 || size(t2, 2) > 1
    error('Input can only have one column');
end
if lfmKern1.inverseWidth ~= lfmKern2.inverseWidth
    error('Kernels cannot be cross combined if they have different inverse widths.')
end

% Parameters of the simulation (in the order providen by kernExtractParam)

m = [lfmKern1.mass lfmKern2.mass];                  % Par. 1
D = [lfmKern1.spring lfmKern2.spring];              % Par. 2
C = [lfmKern1.damper lfmKern2.damper];              % Par. 3
sigma2 = 2/lfmKern1.inverseWidth;                   % Par. 4
sigma = sqrt(sigma2);
S = [lfmKern1.sensitivity lfmKern2.sensitivity];    % Par. 5

alpha = C./(2*m);
omega = sqrt(D./m-alpha.^2);

% Initialization of vectors and matrices

g1 = zeros(1,5);
g2 = zeros(1,5);

% Precomputations
computeH  = cell(4,1);
computeUpsilonMatrix = cell(2,1);
computeUpsilonVector = cell(2,1);
gradientUpsilonMatrix = cell(4,1);
gradientUpsilonVector = cell(4,1);
gamma1_p = alpha(1) + j*omega(1);
gamma1_m = alpha(1) - j*omega(1);
gamma2_p = alpha(2) + j*omega(2);
gamma2_m = alpha(2) - j*omega(2);
preExp1 = zeros(length(t1),2);
preExp2 = zeros(length(t2),2);
preExpg1 = zeros(length(t1),2);
preExpgg1 = zeros(length(t1),2);
preExpg2 = zeros(length(t2),2);
preExpgg2 = zeros(length(t2),2);
preExpt1 = zeros(length(t1),2);
preExpt2 = zeros(length(t2),2);
preGamma(1) = gamma1_p + gamma2_p;
preGamma(2) = gamma1_p + gamma2_m;
preGamma(3) = gamma1_m + gamma2_p;
preGamma(4) = gamma1_m + gamma2_m;
preGamma2 = preGamma.^2;
preConst = 1./preGamma;
preConst2 = 1./(preGamma2);
preFactors(1) = preConst(2) - preConst(1);
preFactors(2) = preConst(3) - preConst(4);
preFactors(3) = preConst(3) - preConst(1);
preFactors(4) = preConst(2) - preConst(4);
preFactors2(1) = -preConst2(2) + preConst2(1);
preFactors2(2) = -preConst2(3) + preConst2(4);
preFactors2(3) = -preConst2(3) + preConst2(1);
preFactors2(4) = -preConst2(2) + preConst2(4);
preExp1(:,1) = exp(-gamma1_p*t1);
preExp1(:,2) = exp(-gamma1_m*t1);
preExp2(:,1) = exp(-gamma2_p*t2);
preExp2(:,2) = exp(-gamma2_m*t2);
preExpg1(:,1) = (2*gamma1_p)*preExp1(:,1);
preExpg1(:,2) = (2*gamma1_m)*preExp1(:,2);
preExpg2(:,1) = (2*gamma2_p)*preExp2(:,1);
preExpg2(:,2) = (2*gamma2_m)*preExp2(:,2);
preExpgg1(:,1) = (gamma1_p^2)*preExp1(:,1);
preExpgg1(:,2) = (gamma1_m^2)*preExp1(:,2);
preExpgg2(:,1) = (gamma2_p^2)*preExp2(:,1);
preExpgg2(:,2) = (gamma2_m^2)*preExp2(:,2);
preExpt1(:,1) = t1.*preExpgg1(:,1);
preExpt1(:,2) = t1.*preExpgg1(:,2);
preExpt2(:,1) = t2.*preExpgg2(:,1);
preExpt2(:,2) = t2.*preExpgg2(:,2);
[computeH{1}, computeUpsilonMatrix{1}, computeUpsilonMatrixLocal{1}] =  lfmComputeH3AA(gamma1_p, gamma1_m, sigma2, t1,t2,preFactors([1 2]), 0);
[computeH{2}, computeUpsilonMatrix{2}, computeUpsilonMatrixLocal{2}] =  lfmComputeH3AA(gamma2_p, gamma2_m, sigma2, t2,t1,preFactors([3 4]), 1);
[computeH{3}, computeUpsilonVector{1}, computeUpsilonVectorLocal{1}] =  lfmComputeH4AA(gamma1_p, gamma1_m, sigma2, t1, preGamma([1 2 4 3]), preExpgg2, 0 );
[computeH{4}, computeUpsilonVector{2}, computeUpsilonVectorLocal{2}] =  lfmComputeH4AA(gamma2_p, gamma2_m, sigma2, t2, preGamma([1 3 4 2]), preExpgg1, 1 );
preKernel = ( computeH{1} + computeH{2}.' + computeH{3} + computeH{4}.');

% Precompute derivatives
gradientUpsilonMatrix{1} = lfmaaGradientUpsilonMatrix(gamma1_p,sigma2,t1, t2, 0, computeUpsilonMatrixLocal{1}{1});
gradientUpsilonMatrix{2} = lfmaaGradientUpsilonMatrix(gamma1_m,sigma2,t1, t2, 0, computeUpsilonMatrixLocal{1}{2});
gradientUpsilonMatrix{3} = lfmaaGradientUpsilonMatrix(gamma2_p,sigma2,t2, t1, 1, computeUpsilonMatrixLocal{2}{1});
gradientUpsilonMatrix{4} = lfmaaGradientUpsilonMatrix(gamma2_m,sigma2,t2, t1, 1, computeUpsilonMatrixLocal{2}{2});
gradientUpsilonVector{1} = lfmapGradientUpsilonVector(gamma1_p,sigma2,t1, computeUpsilonVectorLocal{1}{1});
gradientUpsilonVector{2} = lfmapGradientUpsilonVector(gamma1_m,sigma2,t1, computeUpsilonVectorLocal{1}{2});
gradientUpsilonVector{3} = lfmapGradientUpsilonVector(gamma2_p,sigma2,t2, computeUpsilonVectorLocal{2}{1});
gradientUpsilonVector{4} = lfmapGradientUpsilonVector(gamma2_m,sigma2,t2, computeUpsilonVectorLocal{2}{2});

if lfmKern1.isNormalised
    K0 =  lfmKern1.sensitivity*lfmKern2.sensitivity/(8*sqrt(2)*lfmKern1.mass*lfmKern2.mass*prod(omega));
    K02 = 1/(8*sqrt(2)*prod(m)*prod(omega));
else
    K0 =  sigma*sqrt(pi)*lfmKern1.sensitivity*lfmKern2.sensitivity/(8*lfmKern1.mass*lfmKern2.mass*prod(omega));
    K02 = sigma*sqrt(pi)/(8*prod(m)*prod(omega));
end

% Gradient with respect to m, D and C
for ind_theta = 1:3 % Parameter (m, D or C)
    for ind_par = 0:1 % System (1 or 2)
        % Choosing the right gradients for m, omega, gamma1 and gamma2
        switch ind_theta
            case 1  % Gradient wrt m
                gradThetaM = [1-ind_par ind_par];
                gradThetaAlpha = -C./(2*(m.^2));
                gradThetaOmega = (C.^2-2*m.*D)./(2*(m.^2).*sqrt(4*m.*D-C.^2));
            case 2  % Gradient wrt D
                gradThetaM = zeros(1,2);
                gradThetaAlpha = zeros(1,2);
                gradThetaOmega = 1./sqrt(4*m.*D-C.^2);
            case 3  % Gradient wrt C
                gradThetaM = zeros(1,2);
                gradThetaAlpha = 1./(2*m);
                gradThetaOmega = -C./(2*m.*sqrt(4*m.*D-C.^2));
        end
        gradThetaGamma1 = gradThetaAlpha + j*gradThetaOmega;
        gradThetaGamma2 = gradThetaAlpha - j*gradThetaOmega;
        % Gradient evaluation
        gradThetaGamma11 = [gradThetaGamma1(1) gradThetaGamma2(1)];
        gradThetaGamma2 = [gradThetaGamma1(2) gradThetaGamma2(2)];
        gradThetaGamma1 = gradThetaGamma11;

        if ~ind_par % ind_par = k or d
            matGrad = K0 ...
                * ( lfmGradientH31( preFactors([1 2]), preFactors2([1 2]), gradThetaGamma1, ...
                gradientUpsilonMatrix{1}, gradientUpsilonMatrix{2}, computeUpsilonMatrix{1}{1}, ...
                computeUpsilonMatrix{1}{2}, 1) + ...
                lfmGradientH32( preGamma2,  gradThetaGamma1, computeUpsilonMatrix{2}{1}, ...
                computeUpsilonMatrix{2}{2}, 1).' + ...
                lfmGradientH41( preGamma, preGamma2, gradThetaGamma1, preExpgg2, ...
                gradientUpsilonVector{1}, gradientUpsilonVector{2}, computeUpsilonVector{1}{1},...
                computeUpsilonVector{1}{2}, 1) + ...
                lfmGradientH42AP(preGamma, preGamma2, gradThetaGamma1, preExpg1, preExpgg1, preExpt1, ...
                computeUpsilonVector{2}{1}, computeUpsilonVector{2}{2}).'...
                - (gradThetaM(1+ind_par)/m(1+ind_par) ...
                + gradThetaOmega(1+ind_par)/omega(1+ind_par)) ...
                *preKernel);
        else % ind_par = r or d'
            matGrad = K0 ...
                * ( lfmGradientH31( preFactors([3 4]), preFactors2([3 4]), gradThetaGamma2, ...
                gradientUpsilonMatrix{3}, gradientUpsilonMatrix{4}, computeUpsilonMatrix{2}{1}, ...
                computeUpsilonMatrix{2}{2}, 1).' + ...
                lfmGradientH32( preGamma2([1 3 2 4]),  gradThetaGamma2, computeUpsilonMatrix{1}{1}, ...
                computeUpsilonMatrix{1}{2}, 1) + ...
                lfmGradientH41( preGamma([1 3 2 4]), preGamma2([1 3 2 4]), gradThetaGamma2, preExpgg1, ...
                    gradientUpsilonVector{3}, gradientUpsilonVector{4}, computeUpsilonVector{2}{1},...
                    computeUpsilonVector{2}{2}, 1).' + ...
                lfmGradientH42AP(preGamma([1 3 2 4]), preGamma2([1 3 2 4]), gradThetaGamma2, preExpg2, preExpgg2, preExpt2, ...
                computeUpsilonVector{1}{1}, computeUpsilonVector{1}{2})...
                - (gradThetaM(1+ind_par)/m(1+ind_par) ...
                + gradThetaOmega(1+ind_par)/omega(1+ind_par)) ...
                *preKernel);
        end

        if subComponent
            if size(meanVector,1) ==1,
                matGrad = matGrad*meanVector;
            else
                matGrad = (meanVector*matGrad).';
            end
        end

        % Check the parameter to assign the derivative
        if ind_par == 0
            g1(ind_theta) = sum(sum(matGrad.*covGrad));
        else
            g2(ind_theta) = sum(sum(matGrad.*covGrad));
        end
    end
end

% Gradients with respect to sigma
if lfmKern1.isNormalised
    matGrad = K0*(lfmGradientSigmaH3AA(gamma1_p, gamma1_m, sigma2, t1, t2, preFactors([1 2]), 0)...
        +  lfmGradientSigmaH3AA(gamma2_p, gamma2_m, sigma2, t2, t1, preFactors([3 4]), 1).'...
        +  lfmGradientSigmaH4AA(gamma1_p, gamma1_m, sigma2, t1, preGamma([1 2 4 3]), preExpgg2)...
        +  lfmGradientSigmaH4AA(gamma2_p, gamma2_m, sigma2, t2, preGamma([1 3 4 2]), preExpgg1).');
else
    matGrad = (prod(S)*sqrt(pi)/(8*prod(m)*prod(omega))) ...
        * (preKernel ...
        + sigma*(lfmGradientSigmaH3AA(gamma1_p, gamma1_m, sigma2, t1, t2, preFactors([1 2]), 0)...
        +  lfmGradientSigmaH3AA(gamma2_p, gamma2_m, sigma2, t2, t1, preFactors([3 4]), 1).'...
        +  lfmGradientSigmaH4AA(gamma1_p, gamma1_m, sigma2, t1, preGamma([1 2 4 3]), preExpgg2)...
        +  lfmGradientSigmaH4AA(gamma2_p, gamma2_m, sigma2, t2, preGamma([1 3 4 2]), preExpgg1).' ));
end

if subComponent
    if size(meanVector,1) ==1,
        matGrad = matGrad*meanVector;
    else
        matGrad = (meanVector*matGrad).';
    end
end

g1(4) = sum(sum(matGrad.*covGrad))*(-(sigma^3)/4);
g2(4) = g1(4);

% Gradients with respect to S

matGrad = K02*preKernel;

if subComponent
    if size(meanVector,1) ==1,
        matGrad = matGrad*meanVector;
    else
        matGrad = (meanVector*matGrad).';
    end
end

g1(5) = sum(sum(S(2)*matGrad.*covGrad));
g2(5) = sum(sum(S(1)*matGrad.*covGrad));

g2(4) = 0; % Otherwise is counted twice, temporarly changed by Mauricio Alvarez

g1 = real(g1);
g2 = real(g2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test_bug2524

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY ft_write_mri

data = zeros(5,5,5);
trans = diag([5 5 5 1]);

tmp = tempname;
ft_write_mri([tmp,'.nii'],data,'transform',trans,'dataformat','nifti');
mri = ft_read_mri([tmp,'.nii']);

assert(mri.transform(1,1)==5);

% use a command line tool from FSL
[r,s] = system('which fslhd');
if r==0
  disp(r)
  disp(s)
else
  % use hard-coded FSL version
  s = '/opt/fsl/5.0.9/bin/fslhd';
end
[r,s] = system([deblank(s) ' ' tmp '.nii | grep pixdim']);
disp(r)
disp(s)

assert(numel(strfind(s,'5.0'))==3);

delete([tmp,'.nii']);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function mergedFG = dtiMergeFiberGroups(fg1,fg2,name)
% Merge two fiber groups
%
%   mergedFG = dtiMergeFiberGroups(fg1,fg2,[name])
%
% The parameters are all added into a long list in the params structure.
%
%
% Author: Dougherty, Wandell

if notDefined('fg1'), error('fg1 not defined'); end
if notDefined('fg2'), error('fg2 not defined'); end

% Merge just the fiber data
mergedFG = fg1;
mergedFG.fibers = vertcat(fg1.fibers, fg2.fibers);
if ~(~isfield(fg1, 'subgroup')||~isfield(fg2, 'subgroup'))
mergedFG.subgroup = horzcat(fg1.subgroup, fg2.subgroup);
end
mergedFG.visible = 1;
mergedFG.seeds = [];

% Do the merge of the parameters the parameters
n1 = length(fg1.params);
n2 = length(fg2.params);
for ii=1:n1, mergedFG.params{ii}=fg1.params{ii}; end
for ii=(n1+1):(n1+n2), mergedFG.params{ii} = fg2.params{ii - n1}; end

if ~exist('name','var'),  mergedFG.name = [fg1.name,'_',fg2.name '_MERGED'];
else
    mergedFG.name = name;
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function out = PeubSource(k, p, d, L, accuracy)
%upper bound for the error probability of the binary source
%achievability via sphere covering
%k - source blocklength
%p - source bias
%d - distortion threshold
%L = log M, where M is the number of representation points
%accuracy = accuracy of Csum: 1 - more accurate, 0 - faster

%
%   Created in 2012 by Victoria Kostina (vkostina@caltech.edu)
%

Kbig = 80;

if accuracy < 0
    %normal approximation
    Rd = h(p) - h(d);
    Vs = p*(1-p)*(log2((1-p)/p))^2;
    if p ~= 1/2 && k >= Kbig
        out = Q( (L - k*Rd - log2(k))/sqrt(k*Vs) );
    else
        out = PeubSource(L, k, p, d);
    end
    return;
end

if d == 0
    %lossless
    out = PeubSourceLossless(L, k, p);
elseif p == 1/2
    %FAIR coin flip source
    out = PeubSourceFair(L, k, d);
else
    %arbitrary p source
    out = PeubSourceBiased(L, k, p, d);
end
%        fprintf('SphereCoveringA: L = %i, k = %i, Peub = %f\n', L, k, out);

    function out = PeubSourceLossless(L, k, p)
        csum = 0;
        kstar = k;
        for i = 0:k
            csum = csum + C(k, i, 'ub');
            if log2(csum) > L
                kstar = i - 1;
                break;
            end
        end
        out =  1 - binocdf(kstar,k,p);
    end

    function out = PeubSourceFair(L, k, d)
        if k <= Kbig
            out =  2^L *log(1 - Csum(k, floor(k*d), 'lb', accuracy) /(2^k));
        else
            out = -2^(L-k)*Csum(k, floor(k*d), 'lb', accuracy);
        end
        out = exp(out);
    end

    function out = PeubSourceBiased(L, k, p, d)
        q = max ( 0, (p-d)/(1-2*d));
        if q == 0
            q = p;
        end
        out = 0;
        for t = 0:k
            out = out + binopdfbound(t,k,p, 'ub')*e(k,t);
            if isnan(out)
                break;
            end
        end

        function out = e(k,t)
            out = 0;
            if 1
                for T = 0:k
                    out = out + binopdfbound(T, k, q, 'ub')*Pwithin(k, t, T);
                end
            else
                T = round(q*k); %#ones
                out = Pwithin(k, t, T);
            end

            if k <= Kbig
                out = 2^L*log( 1 - out);
                out = exp(out);
            else
                out = exp(-2^L*out);
            end

        end

        function out = Pwithin(k,t, T)
            %y of type T is within a given x of type t - almost exact lower bound
            t0 = max ( 0, ceil( (t+T - k*d)/2) ); %# ones  in 11111 zone (T ones in codewords)

            out = logC(T, t0, 'lb') + logC(k-T, t - t0, 'lb') - logC(k, t, 'ub');
            out = 2^out;
        end

    end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function obj = Make_f19(obj,filename)
% obj = Make_f19(obj,filename)
% Input a msh class object get the values of the elevations for all times
% in the filename
%
%  Author:      William Pringle
%  Created:     March 30 2018
%               July 19 2018, Updated to do it just for one file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if isempty(obj.op)
    error('No boundary information to interpolate to')
end

%% Do some projection for obj.p
proj = 'Mercator';
m_proj(proj,'lon',[ min(obj.p(:,1)) max(obj.p(:,1)) ],...
            'lat',[ min(obj.p(:,2)) max(obj.p(:,2))])

%% Get boundary info
opedat = obj.op;
b_lon = zeros(opedat.neta,1);
b_lat = zeros(opedat.neta,1);
node_num = zeros(opedat.neta,1);
ns = 1; ne = 0;
for n = 1:opedat.nope
    ne = ne + opedat.nvdll(n);
    node_num(ns:ne) = opedat.nbdv(1:opedat.nvdll(n),n);
    b_lon(ns:ne)    = obj.p(node_num(ns:ne),1);
    b_lat(ns:ne)    = obj.p(node_num(ns:ne),2);
    ns = ne + 1;
end
% Do projection
[b_x,b_y] = m_ll2xy(b_lon,b_lat);

time = ncread(filename,'time');
DT = median(diff(time))*3600; %[s]
L = length(time);
BZ = zeros(length(b_x),L);
disp(['Interpolating ' num2str(L) ' snaps'])
for t = 1:L
    if mod(t,100) == 0
        disp(['Snap #' num2str(t)])
    end
    %% Read the grid and density data
    if t == 1
        lon = ncread(filename,'lon');
        lat = ncread(filename,'lat');
        lon(lon > 180) = lon(lon > 180) - 360;
        [lon,lat] = ndgrid(lon,lat);
        lon_x = reshape(lon,[],1);
        lat_y = reshape(lat,[],1);
        % Delete uncessecary portions
        % First delete by square
        I = find(lon_x < min(obj.p(:,1)) | lon_x > max(obj.p(:,1)) | ...
                 lat_y < min(obj.p(:,2)) | lat_y > max(obj.p(:,2)));
        lon_x(I) = []; lat_y(I) = [];

        % Delete by range search
        Kd = ourKNNsearch([lon_x,lat_y]',[b_lon, b_lat]',20);
        Kd = unique(Kd);

        % The new lon and lat vectors of data
        lon_x = lon_x(Kd); lat_y = lat_y(Kd);

        % Do the projection
        [x,y] = m_ll2xy(lon_x,lat_y);
    end
    zeta = ncread(filename,'surf_el',[1 1 t],[size(lon) 1]);
    Z = reshape(zeta,[],1);
    Z(I) = []; Z = Z(Kd);
    F = scatteredInterpolant(x(~isnan(Z)),y(~isnan(Z)),Z(~isnan(Z)),'natural');
    BZ(:,t) = F(b_x,b_y);
end

%% Make into f19 struct
ts = datenum('2000-01-01') + time(1)/24;
te = datenum('2000-01-01') + time(end)/24;
obj.f19.DataTitle = [datestr(ts) ' -> ' datestr(te)];
obj.f19.ETIMINC = DT;
obj.f19.NumOfNodes = size(BZ,1);
obj.f19.NumOfTimes = size(BZ,2);
obj.f19.Val = BZ(:);
%EOF
end


text:
function tan_test ( )

%*****************************************************************************80
%
%% TAN_TEST tests R4_TAN and R8_TAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TAN_TEST:\n' );
  fprintf ( 1, '  Test TAN_VALUES, R4_TAN, R8_TAN.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X         TAN(X)\n' );
  fprintf ( 1, '                    R4_TAN(X)         Diff\n' );
  fprintf ( 1, '                    R8_TAN(X)         Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = tan_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_tan ( single ( x ) );
    fx3 = r8_tan ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function self = text_test_result(stream, verbosity)
%text_test_result contructor.
%  The constructor creates an object of the class text_test_result. The
%  output is written to the value of stream. The parameter verbosity
%  defines, how much output is written. Possible values are 0, 1 and 2.
%
%  Class Info
%  ==========
%  The class text_test_result is inherited from test_result and prints
%  formatted test results to a stream.
%
%  Example
%  =======
%  Output all results to the Matlab Command Window:
%         result = text_test_result(1, 0)
%
%  See also TEST_RESULT.

self.stream = stream;
if (verbosity == 1)
    self.dots = 1;
    self.show_all = 0;
elseif (verbosity > 1)
    self.dots = 0;
    self.show_all = 1;
else
    self.dots = 0;
    self.show_all = 0;
end;
self.verbosity = verbosity;
result = test_result;
self = class(self, 'text_test_result', result);


text:
function op = smooth_logsumexp(sigma)
% SMOOTH_LOGSUMEXP The function log(sum(exp(x)))
%   returns a smooth function to calculate
%   log( sum( exp(x) ) )
%
% SMOOTH_LOGSUMEXP( SIGMA ) is a scaled version
%   that calclates sigma*log(sum(exp(x/sigma)), for sigma > 0.
%   As sigma --> 0, this becomes a good approximation
%   of max(x).
%   The Lipschitz constant of the gradient is 1/sigma.
%   By default, sigma = 1.
%
% For a fancier version (with offsets),
% see also smooth_logLLogistic.m

if nargin < 1 || isempty(sigma), sigma = 1; end
op = @(x)smooth_logsumexp_impl(x,sigma);

function [ v, g ] = smooth_logsumexp_impl( x, sigma )

% Even for moderate values of x/sigma, exp(x/sigma)
%   will overflow before we have a chance to take
%   its logarithm. So we subtract off the max value
%   and treat it separately:

c    = max(x);
expx = exp((x-c)/sigma);
sum_expx = sum(expx(:));
v = sigma*log(sum_expx) + c;

if nargout > 1,
    g = expx ./ sum_expx;
    % (the factor of e^{-c} cancels from both the numerator
    %  and denominator)
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [brainAreas,vals,labels] = fs_loadCtab(ctabFile)
%
% Load a FreeSurfer .ctab file.  This is generally a file containing the
% labels and indices into a FreeSurfer segmaentation a segmentation.
%
% ctab = fs_loadCtab(ctabFile)
%
% INPUTS:
%    ctabFile - the fullpath to a file ending eith the .ctab exstension.
%               Freesurfer generates two as result of the full segmentation
%               and parcellation process. These files are generally stored
%               in: $SUBJECTS_DIR/<subject>/label/
%
% OUTPUT:
%    brainAreas - is a cell array containing:
%                 brainAreas.val:  The values for each brain areas as coded
%                                  in FreeSurfer.
%                 brainAreas.name: The corresponding labels, names for each
%                                  brain areas.
%    vals       - The values for each brain areas as coded in FreeSurfer
%    labels     - The corresponding labels, names for each
%                 brain areas.
%    colors     - The RGB colors for each brain area in the segmentation.
%
% EXAMPLE:
%    fsDir = getenv('SUBJECTS_DIR');
%    subject = 'pestilli_test';
%    segmentation = 'aparc.annot.a2009s';
%    segFile = fullfile(fsDir,subject,'label',sprintf('%s.ctab',segmentation));
%    ba = fs_loadCtab(segFile,' ',1);
%
% Written by Franco Pestilli (c) Stanford University 2013, Vistasoft

% Read the .ctab file
ctab = importdata(ctabFile,' ',1);

% Extract the labels names and corresponding values, and  colors.
vals             = cellfun(@str2num,ctab.textdata(:,1));
labels           = ctab.textdata(:,2);
brainAreas.name  = labels(:);
brainAreas.val   = vals;
brainAreas.colors= zeros(size(labels,1),3);
brainAreas.colors(vals>0,:) = ctab.data(:,1:3);

end


text:
function r = uplus(p)
% MEAS/UPLUS  Implement +p for meas.

% make a meas called r
r = meas();

% give it the right entries
r.value = p.value;
r.error = p.error;


text:
function a = r8col_swap ( m, n, a, i, j )

%*****************************************************************************80
%
%% R8COL_SWAP swaps columns I and J of an R8COL.
%
%  Example:
%
%    Input:
%
%      M = 3, N = 4, I = 2, J = 4
%
%      A = (
%        1  2  3  4
%        5  6  7  8
%        9 10 11 12 )
%
%    Output:
%
%      A = (
%        1  4  3  2
%        5  8  7  6
%        9 12 11 10 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, real A(M,N), an array of N columns of length M.
%
%    Input, integer I, J, the columns to be swapped.
%
%    Output, real A(M,N), the array, with columns I and J swapped.
%
  if ( i < 1 || n < i || j < 1 || n < j )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8COL_SWAP - Fatal error!\n' );
    fprintf ( 1, '  I or J is out of bounds.\n' );
    fprintf ( 1, '  I =    %d\n', i );
    fprintf ( 1, '  J =    %d\n', j );
    fprintf ( 1, '  N =    %d\n', n );
    error ( 'R8COL_SWAP - Fatal error!' );
  end

  if ( i == j )
    return
  end

  col(1:m) = a(1:m,i)';
  a(1:m,i) = a(1:m,j);
  a(1:m,j) = col(1:m)';

  return
end


text:
function [HP, LP] = sample_patches(im, patch_size, patch_num, upscale)

if size(im, 3) == 3,
    hIm = rgb2gray(im);
else
    hIm = im;
end

% generate low resolution counter parts
lIm = imresize(hIm, 1/upscale, 'bicubic');
lIm = imresize(lIm, size(hIm), 'bicubic');
[nrow, ncol] = size(hIm);

x = randperm(nrow-2*patch_size-1) + patch_size;
y = randperm(ncol-2*patch_size-1) + patch_size;

[X,Y] = meshgrid(x,y);

xrow = X(:);
ycol = Y(:);

if patch_num < length(xrow),
    xrow = xrow(1:patch_num);
    ycol = ycol(1:patch_num);
end

patch_num = length(xrow);

hIm = double(hIm);
lIm = double(lIm);

H = zeros(patch_size^2,     length(xrow));
L = zeros(4*patch_size^2,   length(xrow));

% compute the first and second order gradients
hf1 = [-1,0,1];
vf1 = [-1,0,1]';

lImG11 = conv2(lIm, hf1,'same');
lImG12 = conv2(lIm, vf1,'same');

hf2 = [1,0,-2,0,1];
vf2 = [1,0,-2,0,1]';

lImG21 = conv2(lIm,hf2,'same');
lImG22 = conv2(lIm,vf2,'same');

for ii = 1:patch_num,
    row = xrow(ii);
    col = ycol(ii);

    Hpatch = hIm(row:row+patch_size-1,col:col+patch_size-1);

    Lpatch1 = lImG11(row:row+patch_size-1,col:col+patch_size-1);
    Lpatch2 = lImG12(row:row+patch_size-1,col:col+patch_size-1);
    Lpatch3 = lImG21(row:row+patch_size-1,col:col+patch_size-1);
    Lpatch4 = lImG22(row:row+patch_size-1,col:col+patch_size-1);

    Lpatch = [Lpatch1(:),Lpatch2(:),Lpatch3(:),Lpatch4(:)];
    Lpatch = Lpatch(:);

    HP(:,ii) = Hpatch(:)-mean(Hpatch(:));
    LP(:,ii) = Lpatch;
end

author: thomas-koehler

repo: SupER


text:
function DEMO_niche_construction
% Demo of niche construction using active inference
%__________________________________________________________________________
%
% The free-energy principle is an attempt to explain the structure of the
% agent and its brain, starting from the fact that an agent exists (Friston
% and Stephan, 2007; Friston et al., 2010). More specifically, it can be
% regarded as a systematic attempt to understand the 'fit' between an
% embodied agent and its niche, where the quantity of free-energy is a
% measure for the 'misfit' or disattunement (Bruineberg and Rietveld, 2014)
% between agent and environment. This paper offers a proof-of-principle
% simulation of niche construction under the free-energy principle. The key
% point of this paper is that the minimum of free-energy is not at a point
% in which the agent is maximally adapted to the statistics of a static
% environment, but can better be conceptualized an attracting manifold
% within the joint  agent-environment state-space as a whole, which the
% system tends toward through mutual interaction. We will provide a general
% introduction to active inference and the free-energy principle. Using
% Markov Decision Processes (MDPs), we then describe a canonical generative
% model and the ensuing update equations that minimize free-energy. We then
% apply these equations to simulations of foraging in an environment; in
% which an agent learns the most efficient path to a pre-specified
% location. In some of those simulations, unbeknownst to the agent, the
% environment changes as a function of the activity of the agent (i.e.
% unintentional niche construction occurs). We will show how, depending on
% the relative inertia of the environment and agent, the joint
% agent-environment system moves to different attracting sets of jointly
% minimized free-energy.
%
% see also: spm_MPD_VB_X.m
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

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

%% Set up and preliminaries
%--------------------------------------------------------------------------
rng('default')

label.factor     = {'where'};
label.modality   = {'what','where'};
label.outcome{1} = {'open','closed'};
MAZE =  [...
    0 0 0 0 0 0 0 0;
    1 1 0 1 1 1 1 0;
    1 1 0 1 1 1 1 0;
    1 1 0 1 1 1 1 0;
    1 1 0 1 1 1 1 0;
    0 0 0 1 1 1 1 0;
    0 1 1 1 1 1 1 0;
    0 0 0 0 0 0 0 0];

END   = sub2ind(size(MAZE),1,1);                  % goal or target location
START = sub2ind(size(MAZE),8,1);

% prior beliefs about initial states: D
%--------------------------------------------------------------------------
D{1}  = zeros(numel(MAZE),1);
Ns    = numel(D{1});

% Capture generative process in form of A: Aaux
%--------------------------------------------------------------------------
Aaux{1} = [1 - MAZE(:), MAZE(:)]';                % what
Aaux{2} = eye(Ns,Ns);                             % where
Ng    = numel(Aaux);
for g = 1:Ng
    No(g)  = size(Aaux{g},1);
end

% controlled transitions: B (up, down, left, right, stay)
%--------------------------------------------------------------------------
u    = [1 0; -1 0; 0 1; 0 -1; 0 0];               % allowable actions
nu   = size(u,1);                                 % number of actions
B{1} = zeros(Ns,Ns,nu);
[n,m] = size(MAZE);
for i = 1:n
    for j = 1:m

        % allowable transitions from state s to state ss
        %------------------------------------------------------------------
        s     = sub2ind([n,m],i,j);
        for k = 1:nu
            try
                ss = sub2ind([n,m],i + u(k,1),j + u(k,2));
                B{1}(ss,s,k) = 1;
            catch
                B{1}(s, s,k) = 1;
            end
        end
    end
end

% priors: (negative cost) C:
%--------------------------------------------------------------------------
for g = 1:Ng
    C{g} = zeros(No(g),1);
end

% allowable policies (1 move): V
%--------------------------------------------------------------------------
V = 1:nu;

% basic MDP structure
%--------------------------------------------------------------------------
mdp.V = V;                      % allowable policies
mdp.B = B;                      % transition probabilities
mdp.C = C;                      % preferred outcomes
mdp.D = D;                      % prior over initial states

mdp.label = label;
mdp.alpha = 8;
mdp.tau   = 8;

%% Simulation 1: Learning over trials
%==========================================================================
% Define concentration parameters for environment:
%--------------------------------------------------------------------------
Acp.open   = [4;0];
Acp.closed = [0;4];

% Define concentration parameters for agent:
%--------------------------------------------------------------------------
acp     = [1/8;1/8];

% Generative process: Generate probabilistic mapping from hidden states to
% outcomes: A and matrix of concentration parameters
%--------------------------------------------------------------------------
[A, As] = spm_gen_process(Aaux, Acp);

% Generative model: Generate probabilistic mapping from hidden states to
% outcomes: a matrix
a{1} = repmat(acp, 1,size(A{1},2));
a{2} = A{2}*128;

% MDP structure
%--------------------------------------------------------------------------
mdp.A  = A;
mdp.As = As;
mdp.a  = a;
mdp    = spm_MDP_check(mdp);
ntime  = 16;    % number of moves per path
ntrial = 4;     % number of paths
mdpt   = mdp;

for i = 1:ntrial
    SDP  = spm_maze_active(mdpt,ntime,START,END);
    mdpt = SDP(end);
    trial(i).MDP = SDP;
end

% show results - behavioural
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1'); clf
spm_print_agentstrials(trial, START, END, a)

%% Simulation 2: sensitivity to concentration parameters
%==========================================================================

% Show typical trial for
%--------------------------------------------------------------------------
acp    = [1/8,1/2,2;1/8,1/2,2];   % Initial concentration parameters of agent (cp1 = cp2)
ntime  = 16;                      % number of timesteps per trial
na     = length(acp);
ntrial = 4;
trial  = [];

for i = 1:na
    mdp.a{1} = repmat(acp(:,i), 1,size(A{1},2));
    mdpt = mdp;
    for j = 1:ntrial
        SDP = spm_maze_active(mdpt,ntime,START,END);
        mdpt = SDP(end);
    end
    trial(i).MDP = SDP;
end
spm_figure('GetWin','Figure 2'); clf
spm_plot_3a(trial, START, END, acp);

%% Simulation 3: sensitivity of concentration parameters agent and environment
%======================================================================
acp   = [1/8,4; 1/8,4];      % Initial concentration parameters of agent (cp1 = cp2)

Acp(1).open   = [1;0];
Acp(1).closed = [0;1];
Acp(2).open   = [16;0];
Acp(2).closed = [0;16];

ntime  = 16;                 % number of moves per path
ntrial = 32;                 % number of paths

% Inital concentration parameters
%--------------------------------------------------------------------------
na    = size(acp,2);
nenv  = size(Acp,2);
agent = [];

for i = 1:na
    mdp.a{1} = repmat(acp(:,i), 1,size(A{1},2));
    for j = 1:nenv
        [A,As] = spm_gen_process(Aaux, Acp(j));
        mdp.A  = A;
        mdp.As = As;
        mdpt   = mdp;
        for k = 1:ntrial
            SDP  = spm_maze_active(mdpt,ntime,START,END);
            mdpt = SDP(end);
            agent(i).env(j).trial(k).MDP = SDP;
        end
    end
end
spm_figure('GetWin','Figure 3'); clf
spm_plot_3b(agent, START, END);

%% Simulation 4: Singular value decomposition
%==========================================================================
spm_figure('GetWin','Figure 4'); clf
nenv  = 2;
ntt   = ntime*ntrial;
c     = linspace(10,ntt/4,ntt);
for i = 1:na
    for j = 1:nenv
        Pag  = [];
        Penv = [];
        for k = 1:ntrial
            SDP = agent(i).env(j).trial(k).MDP;
            for l = 1:ntime

                % get expectations
                %----------------------------------------------------------
                Pa   = SDP(l).a{1};
                Pa   = Pa/diag(sum(Pa));
                Pag  = [Pag, Pa(1,:)'];
                Pa   = SDP(l).As;
                Pa   = Pa/diag(sum(Pa));
                Penv = [Penv, Pa(1,:)'];
            end
        end
        Pag  = Pag  -  Pag(:,end)*ones(1,ntt);
        Penv = Penv - Penv(:,end)*ones(1,ntt);
        X    = [Pag'; Penv'];

        % eigenvariates
        %------------------------------------------------------------------
        subplot(2,2,2*(i-1) + j)
        [U,S] = svd(X);
        U   = U*S;
        V1a = U((1:ntt),1);
        V2a = U((1:ntt),2);
        V1e = U((1:ntt) + ntt,1);
        V2e = U((1:ntt) + ntt,2);

        scatter(V1a, V2a, [], c, 'filled'), hold on
        scatter(V1e, V2e, [], c, 'ro');
        axis([-1 1 -1 1]*2.5), axis square, hold off
        xlabel('First principle component')
        ylabel('Second principle component')
        title('phenotypic trajectories','FontSize',16)
    end
end

% get and show various free energies
%--------------------------------------------------------------------------
for i = 1:na
    for j = 1:nenv
        for k = 1:ntrial
            SDP = agent(i).env(j).trial(k).MDP;
            for l = 1:ntime;
                Gu     = log(spm_softmax(SDP(l).G));
                r(l,k) = mean(SDP(l).rt);
                f(l,k) = mean(SDP(l).H);
                g(l,k) = mean(SDP(l).P'*Gu);
            end
        end
        R(i,j,:) = sum(r);
        F(i,j,:) = sum(f);
        G(i,j,:) = sum(g);
    end
end

% plot
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 5'); clf
col = {'k','r';':k',':r'};

for i = 1:na
    for j = 1:nenv
        subplot(2,2,1)
        plot(spm_vec(F(i,j,:)),col{i,j}), hold on
        xlabel('exposure'), ylabel('seconds')
        title('(negative) free energy','FontSize',16)
        axis square, spm_axis tight, %set(gca,'YLim',[-3 1])

        subplot(2,2,2)
        plot(spm_vec(G(i,j,:)),col{i,j}), hold on
        xlabel('exposure'), ylabel('seconds')
        title('expected free energy','FontSize',16)
        axis square, spm_axis tight, %set(gca,'YLim',[-16 2])

        subplot(2,1,2)
        plot(spm_vec(R(i,j,:)),col{i,j}), hold on
        xlabel('exposure'), ylabel('seconds')
        title('reaction time','FontSize',16)
        axis square, spm_axis tight

    end
end

legend({'A:low-E:low','A:low-E:high','A:high-E:low','A:high-E:high'})

%%
return

function MDP = spm_maze_active(mdp,N,START,END)
% FORMAT MDP = spm_maze_search(mdp,N,START,END,alpha,beta,z)
% mdp   - MDP structure
% N     - number of moves (i.e., default 8)
% START - index of intial state (default 1)
% END   - index of target state (default 1)
% MDP   - MDP structure array

% preliminaries
%--------------------------------------------------------------------------
try, N;     catch, N     = 8;   end
try, START; catch, START = 1;   end
try, END;   catch, END   = 1;   end

% initialise concentration parameters: a (if unspecified)
%--------------------------------------------------------------------------
if ~isfield(mdp,'o')
    mdp.o = [];
end
if ~isfield(mdp,'u')
    mdp.u = [];
end
mdp.s    = START;

% Evaluate a sequence moves
%==========================================================================
for i = 1:N

    % Evaluate preferred states (subgoals) on the basis of current beliefs
    %----------------------------------------------------------------------
    mdp.C{2} = spm_maze_cost(mdp,END);

    % proceed with subsequent move
    %----------------------------------------------------------------------
    MDP(i)   = spm_MDP_VB_X(mdp);
    mdp      = MDP(i);
    mdp.s    = mdp.s(:,end);
    mdp.D{1} = MDP(i).X{1}(:,end);
    mdp.o    = [];
    mdp.u    = [];
    mdp.a{1} = MDP(i).a{1};

    % UPDATE: Change A matrix dependent on current position of agent
    %----------------------------------------------------------------------
    pos      = MDP(i).s(:,end);

    % Update environment
    %----------------------------------------------------------------------
    MDP(i).As(1,pos) = MDP(i).As(1,pos) + 1;

    % place in MDP structure
    %----------------------------------------------------------------------
    mdp.As   = MDP(i).As;
    mdp.A{1} = MDP(i).As/diag(sum(MDP(i).As));

end
return

function C = spm_maze_cost(MDP,END)
% Evaluate subgoals using
%==========================================================================
START = MDP.s(1);
if isfield(MDP,'a')
    Q = MDP.a{1};
else
    Q = MDP.A{1};
end
Q   = Q/diag(sum(Q));
Q   = Q(1,:);                                % open states
P   = diag(Q)*any(MDP.B{1},3);               % allowable transitions
ns  = length(Q);                             % number of states
X   = zeros(ns,1);X(START) = 1;              % initial state
Y   = zeros(ns,1);Y(END)   = 1;              % target state

% Preclude transitions to closed states and evaluate graph Laplacian
%--------------------------------------------------------------------------
P   = P - diag(diag(P));
P   = P - diag(sum(P));
P   = expm(P);

% evaluate (negative) cost as a path integral conjunctions
%--------------------------------------------------------------------------
for t = 1:size(MDP.V,1)
    X = P*X;
end
X     = X > exp(-3);                    % Default cut-off is -3
C     = log(X.*(P*Y) + exp(-32));

return

function spm_plot_3a (trial, START, GOAL, cp_ag)
% unpack properties
ntrial = size(trial,2);
ntime  = size(trial(1).MDP,2);
A  = spm_vec(trial(1).MDP(1).A{1}(1,:));
ns = numel(A);
ni = sqrt(ns);

for j = 1:ntrial

    A  = spm_vec(trial(j).MDP(end).A{1}(1,:));
    A  = reshape(A,ni,ni);
    subplot(2,ntrial,j), imagesc(A, [0,1]), axis image
    title(['CP = ' num2str(cp_ag(1,j))], 'fontsize',14)
    hold on, colormap gray

    % Plot trajectory
    %----------------------------------------------------------------------
    for k = 1:ntime
        step    = trial(j).MDP(k).s;
        [l1,m1] = ind2sub([ni,ni],step(1));
        [l2,m2] = ind2sub([ni,ni],step(2));
        plot([m1,m2], [l1,l2],':m');
    end

    [p,q] = ind2sub([ni,ni],START);
    plot(q,p,'.','MarkerSize',32,'Color','g');
    [p,q] = ind2sub([ni,ni],GOAL);
    plot(q,p,'.','MarkerSize',32,'Color','r');
    [p,q] = ind2sub([ni,ni],trial(j).MDP(end).s(2));
    plot(q,p,'.','MarkerSize',32,'Color','b');

    hold off

    % Plot Likelihood matrix
    %----------------------------------------------------------------------
    Q   = trial(j).MDP(end).a{1};
    Q   = Q/diag(sum(Q));
    Q   = Q(1,:);
    a   = reshape(Q(:),ni,ni);
    subplot(2,ntrial,j+3), imagesc(a, [0,1]), axis image
    title('Likelihood','fontsize',14);

end

return

function spm_plot_3b (agent, START, GOAL)

% unpack properties
%--------------------------------------------------------------------------
nag   = size(agent,2);
nenv  = size(agent(1).env,2);
ntime = size(agent(1).env(1).trial(1).MDP,2);
A     = spm_vec(agent(1).env(1).trial(1).MDP(1).A{1}(1,:));
ns    = numel(A);
ni    = sqrt(ns);

for i = 1:nenv
    for j = 1:nag
        A  = spm_vec(agent(j).env(i).trial(end).MDP(end).A{1}(1,:));
        A  = reshape(A,ni,ni);
        subplot(nenv,nag*2,(2*(j-1))+1+(2*nenv*(i-1)))
        imagesc(A, [0,1]), axis image

        hold on
        colormap gray
        % Plot trajectory
        %------------------------------------------------------------------
        for k = 1:ntime
            step    = agent(j).env(i).trial(end).MDP(k).s;
            [l1,m1] = ind2sub([ni,ni],step(1));
            [l2,m2] = ind2sub([ni,ni],step(2));
            plot([m1,m2], [l1,l2],':r');
        end

        [p,q] = ind2sub([ni,ni],START);
        plot(q,p,'.','MarkerSize',32,'Color','g');
        [p,q] = ind2sub([ni,ni],GOAL);
        plot(q,p,'.','MarkerSize',32,'Color','r');
        [p,q] = ind2sub([ni,ni],agent(j).env(i).trial(end).MDP(end).s(2));
        plot(q,p,'.','MarkerSize',32,'Color','b');

        % Plot Likelihood matrix
        %------------------------------------------------------------------
        Q   = agent(j).env(i).trial(end).MDP(end).a{1};
        Q   = Q/diag(sum(Q));
        Q   = Q(1,:);
        a   = reshape(Q(:),ni,ni);
        subplot(nenv,nag*2,(2*(j-1))+2+(2*nenv*(i-1))), imagesc(a, [0,1]), axis image

        hold off

    end
end
return

function [A, As] = spm_gen_process(Aaux, Acp)
% This function generates the observation matrix A of the generative
% process based on the layout of the maze Aaux and the environmental
% concentration parameters Acp
%--------------------------------------------------------------------------
initwh = Acp.open;
initbl = Acp.closed;
As     = initwh*Aaux{1}(1,:)+initbl*Aaux{1}(2,:);
Ap     = repmat(sum(As,1),2,1);

A{1} = As./Ap;    % Expected value of dirichlet is alpha1/(alpha1+alpha2);
A{2} = Aaux{2};

return

function spm_print_agentstrials(trial, START, GOAL,ainit)
%--------------------------------------------------------------------------
ntrial = size(trial,2);
ntime  = size(trial(1).MDP,2);
A  = spm_vec(trial(1).MDP(1).A{1}(1,:));
a  = ainit;
ns = numel(A);
ni = sqrt(ns);

% plot initial maze
%--------------------------------------------------------------------------
A  = spm_vec(trial(1).MDP(1).A{1}(1,:));
A  = reshape(A,ni,ni);
subplot(ntrial+1,2,1), imagesc(A, [0,1]), axis image
title('Initial maze', 'fontsize',14)
hold on
colormap gray
try
    [p,q] = ind2sub([ni,ni],START);
    plot(q,p,'.','MarkerSize',32,'Color','g');
    [p,q] = ind2sub([ni,ni],GOAL);
    plot(q,p,'.','MarkerSize',32,'Color','r');
end
hold off

% plot initial likelihood
%--------------------------------------------------------------------------
Q = a{1};
Q = Q/diag(sum(Q));
Q = Q(1,:);
a = reshape(Q(:),ni,ni);
subplot(ntrial+1,2,2), imagesc(a, [0,1]), axis image
title('Initial Likelihood','fontsize',14);

for j = 1:ntrial
    A = spm_vec(trial(j).MDP(end).A{1}(1,:));
    A = reshape(A,ni,ni);
    subplot(ntrial+1,2,1+(2*j)), imagesc(A, [0,1]), axis image
    title(['Trial ' num2str(j)], 'fontsize',14)
    hold on
    % Plot trajectory
    %----------------------------------------------------------------------
    for k = 1:ntime
        step = trial(j).MDP(k).s;
        [l1,m1]= ind2sub([ni,ni],step(1));
        [l2,m2]= ind2sub([ni,ni],step(2));
        p = plot([m1,m2], [l1,l2],':m');
    end

    % try
    %----------------------------------------------------------------------
    [p,q] = ind2sub([ni,ni],START);
    plot(q,p,'.','MarkerSize',32,'Color','g');
    [p,q] = ind2sub([ni,ni],GOAL);
    plot(q,p,'.','MarkerSize',32,'Color','r');
    [p,q] = ind2sub([ni,ni],trial(j).MDP(end).s(2));
    plot(q,p,'.','MarkerSize',32,'Color','b');

    hold off

    % Plot Likelihood matrix
    %----------------------------------------------------------------------
    Q     = trial(j).MDP(end).a{1};
    Q     = Q/diag(sum(Q));
    Q     = Q(1,:);
    a     = reshape(Q(:),ni,ni);
    subplot(ntrial+1,2,2+(2*j)), imagesc(a, [0,1]), axis image
    title(['Likelihood after trial ' num2str(j)],'fontsize',14);
end

return

author: spm


text:
function clu = updateclu(clu,ix,newclu)
% updates identities in parts of a clu file

%% INPUTS

% clu - N x 1 vector of cluster identities for each extracted waveform
% ix - M x 1 vector of indices where new identities will be assigned
% newclu - integer value for the new cluster identity

%% OUTPUTS
% clu - N x 1 vector of updated cluster identities

% updated 2017/02/24 by David Tingley

for ii=1:length(ix)
    clu(clu==ix(ii))=newclu;
end

% ix = unique(clu);
% for ii=1:length(ix)
%     clu(clu==ix(ii))=ii;
% end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode