text:
classdef maskPaintBrush < handle

    properties
        handles
        position                    %[x y radius]
        Visible = 'on';
        color = [206 82 65]./256;
        nPoints = 20;
        oldWBMF
        oldPTR
    end

    events
    end

    methods
        %Contructor
        function brush = maskPaintBrush(tool)

            %get the parent axes handle
            h = getHandles(tool);
            brush.handles.parent = h.Axes;
            brush.handles.fig = h.fig;
            brush.handles.tool=tool;

            %get the mouse location
            cp = get(brush.handles.parent,'CurrentPoint'); cp=[cp(1,1) cp(1,2)];
            position = [cp(1,1) cp(1,2) brush.handles.tool.brushsize];

            %Create the circle object
            pos = [position(1)-position(3) position(2)-position(3) 2*position(3) 2*position(3)];
            brush.handles.circle = rectangle('Position',pos,'Parent',brush.handles.parent,'EdgeColor',brush.color,'LineWidth',1.5,'Curvature',[1 1],'PickableParts','all');

            %Get the old window button motion function
            brush.oldWBMF = get(brush.handles.fig,'WindowButtonMotionFcn');

            %Get the old pointer style
            brush.oldPTR.ptr = get(brush.handles.fig,'Pointer');
            brush.oldPTR.CData = get(brush.handles.fig,'PointerShapeCData');

            %Make the pointer invisible
            %set(brush.handles.fig,'Pointer','custom','PointerShapeCData',nan(16))

            %Set the new window button motion function
            fun = @(src,evnt) ButtonMotionFunction(src,evnt,brush,'No Click',[]);
            set(brush.handles.fig,'WindowButtonMotionFcn',fun)

            %Set the mouse click function
            fun=@(hObject,eventdata) buttonDownFunction(hObject,eventdata,brush);
            set(brush.handles.circle,'ButtonDownFcn',fun)

            %set the position property of the brush
            brush.position=position;

        end

        function delete(brush)  %destructor
            try
                delete(brush.handles.circle);
                set(brush.handles.fig,'Pointer',brush.oldPTR.ptr,'PointerShapeCData',brush.oldPTR.CData,'WindowButtonMotionFcn',brush.oldWBMF);
            end
        end

        function set.position(brush,position)
            brush.position=position;
            pos = [position(1)-position(3) position(2)-position(3) 2*position(3) 2*position(3)];
            set(brush.handles.circle,'Position',pos);
        end

        function set.Visible(brush,visible)
            set(brush.handles.cirlce,'Visible',visible);
        end

        function mask = getBrushMask(brush)
            %This function creates a binary mask based on the current
            %position of the mask

            t=linspace(0,2*pi,brush.nPoints); %elliptical equation is parameterized by t
            x = brush.position(3)*cos(t)+brush.position(1);
            y= brush.position(3)*sin(t)+brush.position(2);
            S = getImageSize(brush.handles.tool);
            mask = poly2mask(x,y,S(1),S(2));

        end

        function answer = isInAxis(brush)

            xlim=get(brush.handles.parent,'Xlim');
            ylim=get(brush.handles.parent,'Ylim');

            if brush.position(1)>=xlim(1) && brush.position(1)<=xlim(2) && brush.position(2)>=ylim(1) && brush.position(2)<=ylim(2)
                answer = true;
            else
                answer = false;
            end

        end
    end

end

function buttonUpFunction(hObject,eventdata,WBMF_old,WBUF_old,brush)
set(hObject,'WindowButtonMotionFcn',WBMF_old,'WindowButtonUpFcn',WBUF_old);
notify(brush.handles.tool,'maskChanged')

end

function buttonDownFunction(hObject,eventdata,brush)

WBMF_old = get(brush.handles.fig,'WindowButtonMotionFcn');
WBUF_old = get(brush.handles.fig,'WindowButtonUpFcn');
switch get(brush.handles.fig,'SelectionType')
    case 'normal'   %left click (paint on the mask)
        fun = @(src,evnt) ButtonMotionFunction(src,evnt,brush,'Left Click',[]);
        fun2=@(src,evnt) buttonUpFunction(src,evnt,WBMF_old,WBUF_old,brush);
        set(brush.handles.fig,'WindowButtonMotionFcn',fun,'WindowButtonUpFcn',fun2)
        fun([],[]);

    case 'alt'      %right click
        fun = @(src,evnt) ButtonMotionFunction(src,evnt,brush,'Right Click',[]);
        fun2=@(src,evnt) buttonUpFunction(src,evnt,WBMF_old,WBUF_old,brush);
        set(brush.handles.fig,'WindowButtonMotionFcn',fun,'WindowButtonUpFcn',fun2)
        fun([],[]);
    case 'extend'   %center click
        data.bp=get(0,'PointerLocation');
        data.r = brush.position(3);
        fun = @(src,evnt) ButtonMotionFunction(src,evnt,brush,'Middle Click',data);
        fun2=@(src,evnt) buttonUpFunction(src,evnt,WBMF_old,WBUF_old,brush);
        set(brush.handles.fig,'WindowButtonMotionFcn',fun,'WindowButtonUpFcn',fun2)
        fun([],[]);
end

end

function ButtonMotionFunction(src,evnt,brush,tag,data)

switch tag
    case 'No Click'
        if ~isvalid(brush), return; end
        if isempty(brush.position), return; end
        brush.oldWBMF(src,evnt)

        cp = get(brush.handles.parent,'CurrentPoint'); cp=[cp(1,1) cp(1,2)];
        brush.position = [cp(1) cp(2) brush.position(3)];

        %Check if the cursor is within the axis
        if isInAxis(brush)
             set(brush.handles.fig,'Pointer','custom','PointerShapeCData',nan(16))
        else
            set(brush.handles.fig,'Pointer',brush.oldPTR.ptr,'PointerShapeCData',brush.oldPTR.CData)

        end

    case 'Left Click'
        %moves the circle
        cp = get(brush.handles.parent,'CurrentPoint'); cp=[cp(1,1) cp(1,2)];
        brush.position = [cp(1) cp(2) brush.position(3)];

        %get the mask for the new brush position
        mask = getBrushMask(brush);

        %get the current mask of the imtool3D object
        maskOld = getCurrentMaskSlice(brush.handles.tool);

        %Combine the two masks
        mask = mask | maskOld;

        %Update the mask of the tool
        setCurrentMaskSlice(brush.handles.tool,mask)

    case 'Right Click'
        %moves the circle
        cp = get(brush.handles.parent,'CurrentPoint'); cp=[cp(1,1) cp(1,2)];
        brush.position = [cp(1) cp(2) brush.position(3)];

        %get the mask for the new brush position
        mask = getBrushMask(brush);

        %get the current mask of the imtool3D object
        maskOld = getCurrentMaskSlice(brush.handles.tool);

        %Combine the two masks
        mask = maskOld   &   ~(mask & maskOld);

        %Update the mask of the tool
        setCurrentMaskSlice(brush.handles.tool,mask)

    case 'Middle Click'
        %get the current mouse position
        cp = get(0,'PointerLocation');

        %get the difference in the y direction
        d = .25*(data.bp(2)-cp(2));

        rnew = data.r+d;
        if rnew<1
            rnew=1;
        end

        brush.position(3)=rnew;
        brush.handles.tool.brushsize = rnew;

end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [methodinfo,structs,enuminfo,ThunkLibName]=hidapi64mac_proto
%HIDAPI64MAC_PROTO Create structures to define interfaces found in 'hidapi'.

%This function was generated by loadlibrary.m parser version 1.1.6.38 on Thu Jan 23 17:27:34 2014
%perl options:'hidapi.i -outfile=hidapi64mac_proto.m -thunkfile=hidapi64_thunk_maci64.c -header=hidapi.h'
ival={cell(1,0)}; % change 0 to the actual number of functions to preallocate the data.
structs=[];enuminfo=[];fcnNum=1;
fcns=struct('name',ival,'calltype',ival,'LHS',ival,'RHS',ival,'alias',ival,'thunkname', ival);
MfilePath=fileparts(mfilename('fullpath'));
ThunkLibName=fullfile(MfilePath,'hidapi64_thunk_maci64');
% int hid_init ( void );
fcns.thunkname{fcnNum}='int32voidThunk';fcns.name{fcnNum}='hid_init'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}=[];fcnNum=fcnNum+1;
% int hid_exit ( void );
fcns.thunkname{fcnNum}='int32voidThunk';fcns.name{fcnNum}='hid_exit'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}=[];fcnNum=fcnNum+1;
% struct hid_device_info * hid_enumerate ( unsigned short vendor_id , unsigned short product_id );
fcns.thunkname{fcnNum}='voidPtruint16uint16Thunk';fcns.name{fcnNum}='hid_enumerate'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='hid_device_infoPtr'; fcns.RHS{fcnNum}={'uint16', 'uint16'};fcnNum=fcnNum+1;
% void hid_free_enumeration ( struct hid_device_info * devs );
fcns.thunkname{fcnNum}='voidvoidPtrThunk';fcns.name{fcnNum}='hid_free_enumeration'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}=[]; fcns.RHS{fcnNum}={'hid_device_infoPtr'};fcnNum=fcnNum+1;
% hid_device * hid_open ( unsigned short vendor_id , unsigned short product_id , wchar_t * serial_number );
fcns.thunkname{fcnNum}='voidPtruint16uint16voidPtrThunk';fcns.name{fcnNum}='hid_open'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='hid_device_Ptr'; fcns.RHS{fcnNum}={'uint16', 'uint16', 'int32Ptr'};fcnNum=fcnNum+1;
% hid_device * hid_open_path ( const char * path );
fcns.thunkname{fcnNum}='voidPtrcstringThunk';fcns.name{fcnNum}='hid_open_path'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='hid_device_Ptr'; fcns.RHS{fcnNum}={'cstring'};fcnNum=fcnNum+1;
% int hid_write ( hid_device * device , const unsigned char * data , size_t length );
fcns.thunkname{fcnNum}='int32voidPtrvoidPtrulongThunk';fcns.name{fcnNum}='hid_write'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}={'hid_device_Ptr', 'uint8Ptr', 'ulong'};fcnNum=fcnNum+1;
% int hid_read_timeout ( hid_device * dev , unsigned char * data , size_t length , int milliseconds );
fcns.thunkname{fcnNum}='int32voidPtrvoidPtrulongint32Thunk';fcns.name{fcnNum}='hid_read_timeout'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}={'hid_device_Ptr', 'uint8Ptr', 'ulong', 'int32'};fcnNum=fcnNum+1;
% int hid_read ( hid_device * device , unsigned char * data , size_t length );
fcns.thunkname{fcnNum}='int32voidPtrvoidPtrulongThunk';fcns.name{fcnNum}='hid_read'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}={'hid_device_Ptr', 'uint8Ptr', 'ulong'};fcnNum=fcnNum+1;
% int hid_set_nonblocking ( hid_device * device , int nonblock );
fcns.thunkname{fcnNum}='int32voidPtrint32Thunk';fcns.name{fcnNum}='hid_set_nonblocking'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}={'hid_device_Ptr', 'int32'};fcnNum=fcnNum+1;
% int hid_send_feature_report ( hid_device * device , const unsigned char * data , size_t length );
fcns.thunkname{fcnNum}='int32voidPtrvoidPtrulongThunk';fcns.name{fcnNum}='hid_send_feature_report'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}={'hid_device_Ptr', 'uint8Ptr', 'ulong'};fcnNum=fcnNum+1;
% int hid_get_feature_report ( hid_device * device , unsigned char * data , size_t length );
fcns.thunkname{fcnNum}='int32voidPtrvoidPtrulongThunk';fcns.name{fcnNum}='hid_get_feature_report'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}={'hid_device_Ptr', 'uint8Ptr', 'ulong'};fcnNum=fcnNum+1;
% void hid_close ( hid_device * device );
fcns.thunkname{fcnNum}='voidvoidPtrThunk';fcns.name{fcnNum}='hid_close'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}=[]; fcns.RHS{fcnNum}={'hid_device_Ptr'};fcnNum=fcnNum+1;
% int hid_get_manufacturer_string ( hid_device * device , wchar_t * string , size_t maxlen );
fcns.thunkname{fcnNum}='int32voidPtrvoidPtrulongThunk';fcns.name{fcnNum}='hid_get_manufacturer_string'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}={'hid_device_Ptr', 'int32Ptr', 'ulong'};fcnNum=fcnNum+1;
% int hid_get_product_string ( hid_device * device , wchar_t * string , size_t maxlen );
fcns.thunkname{fcnNum}='int32voidPtrvoidPtrulongThunk';fcns.name{fcnNum}='hid_get_product_string'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}={'hid_device_Ptr', 'int32Ptr', 'ulong'};fcnNum=fcnNum+1;
% int hid_get_serial_number_string ( hid_device * device , wchar_t * string , size_t maxlen );
fcns.thunkname{fcnNum}='int32voidPtrvoidPtrulongThunk';fcns.name{fcnNum}='hid_get_serial_number_string'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}={'hid_device_Ptr', 'int32Ptr', 'ulong'};fcnNum=fcnNum+1;
% int hid_get_indexed_string ( hid_device * device , int string_index , wchar_t * string , size_t maxlen );
fcns.thunkname{fcnNum}='int32voidPtrint32voidPtrulongThunk';fcns.name{fcnNum}='hid_get_indexed_string'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32'; fcns.RHS{fcnNum}={'hid_device_Ptr', 'int32', 'int32Ptr', 'ulong'};fcnNum=fcnNum+1;
% const wchar_t * hid_error ( hid_device * device );
fcns.thunkname{fcnNum}='voidPtrvoidPtrThunk';fcns.name{fcnNum}='hid_error'; fcns.calltype{fcnNum}='Thunk'; fcns.LHS{fcnNum}='int32Ptr'; fcns.RHS{fcnNum}={'hid_device_Ptr'};fcnNum=fcnNum+1;
structs.hid_device_.members=struct('');
structs.hid_device_info.members=struct('path', 'cstring', 'vendor_id', 'uint16', 'product_id', 'uint16', 'serial_number', 'int32Ptr', 'release_number', 'uint16', 'manufacturer_string', 'int32Ptr', 'product_string', 'int32Ptr', 'usage_page', 'uint16', 'usage', 'uint16', 'interface_number', 'int32', 'next', 'hid_device_infoPtr');
methodinfo=fcns;

author: petercorke

repo: robotics-toolbox-matlab


text:
function [m] = pm2m(pm)
% Convert length from picometers to meters.
% Chad A. Greene 2012
m = pm*1e-12;

author: Sable

repo: mcbench-benchmarks


text:
%--- help for dsge/forecast_real_time ---
%
%  Forecast from each point in time
%
%  ::
%
%     [ts_fkst,ts_rmse,rmse,Updates]=forecast_real_time(obj)
%     [ts_fkst,ts_rmse,rmse,Updates]=forecast_real_time(obj,varargin)
%
%  Args:
%
%     obj (dsge | svar | rfvar): model object
%     varargin : valid optional inputs coming in pairs. The main inputs
%       of interest for changing the default behavior are:
%
%       - **forecast_rt_nsteps** [integer] : number of periods ahead
%
%  Returns:
%     :
%
%     - **ts_fkst** [struct] : fields are forecasts in the form of ts objects
%       for the different endogenous variables
%
%     - **ts_rmse** [ts|struct] : if only one object is processed, the output
%       is a ts. If instead several objects are processed, fields are RMSEs in
%       the form of ts objects for the different endogenous variables
%
%     - **rmse** [matrix] : RMSEs for the different endogenous variables
%
%     - **Updates** [struct] : fields are the updated (in a filtering sense) in
%       the form of ts objects for the different endogenous variables
%
%  See also:
%     - plot_real_time
%
%

author: jmaih


text:
% Originally Coded By: Niall Mangan
% Copyright 2017, All Rights Reserved
% Code by Niall Mangan for paper "Inferring biological networks by sparse
% identification of nonlinear dynamics"
% by N. M. Mangan S. L. Brunton, J. L. Proctor, and J. N. Kutz
% Original Code: https://github.com/niallmm/iSINDy
%%
% This file is modified to compare the performance of i-SINDy on
% insufficient data.
% Modified by:K, 2019/07/16
% Find equations for yeast glycolysis state variable 4
%%
clc;close all;clear all;

addpath('./utils');
addpath('./bioutils');

% define libarary parameters
laurentorder = 0;
polyorder = 3;
usesine = 0;
dyorder = 1;

% clear variables from other solutions.
clear Theta Thetastring Xi indTheta lambdavec numterms errorv indopt
clear indTheta1 Xi1 numterms1 nT

% Here we load the previously simulated data
load('TrainingData.mat')

% Define the new data length
percent=0.006;new_length=round(percent*length(xt));

% Shuffel the original data
Sequence=randperm(size(xt,1));
xt=xt(Sequence,:);
dxt=dxt(Sequence,:);

% Assign the value to the new variables
Data=xt(1:new_length,:);dData=dxt(1:new_length,:);

% Define the number of states
n=size(Data,2);

% pool Data  (i.e., build library of nonlinear time series)
[Theta, Thetastring] = poolDatady(Data,n,polyorder,usesine, laurentorder, dData(:,4), dyorder);
% %initial lambda value, which is the value used for soft thresholding in ADM

tol = 2e-3;
lambda = 8e-3;

jj = 1; % counter
num= 1; % initialize the number of nonzero terms found for the lambda
errorvec= 0;

MaxIter = 1e3;

% for now calculate null space using null function
nT = null(Theta);

% Define the parameters for plooting
plottag=2;

[indTheta1, Xi1, numterms1] = ADMinitvary(nT,lambda,MaxIter,tol, plottag);
Thetastring(indTheta1)'
n0Xi = Xi1(Xi1~=0);
n0Xi/n0Xi(end)
save('Results/4rd_state_variable.mat')

author: dynamicslab

repo: SINDy-PI


text:
function varargout = lvd_EditElectricalPowerSinksGUI(varargin)
% LVD_EDITELECTRICALPOWERSINKSGUI MATLAB code for lvd_EditElectricalPowerSinksGUI.fig
%      LVD_EDITELECTRICALPOWERSINKSGUI, by itself, creates a new LVD_EDITELECTRICALPOWERSINKSGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITELECTRICALPOWERSINKSGUI returns the handle to a new LVD_EDITELECTRICALPOWERSINKSGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITELECTRICALPOWERSINKSGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITELECTRICALPOWERSINKSGUI.M with the given input arguments.
%
%      LVD_EDITELECTRICALPOWERSINKSGUI('Property','Value',...) creates a new LVD_EDITELECTRICALPOWERSINKSGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditElectricalPowerSinksGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditElectricalPowerSinksGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help lvd_EditElectricalPowerSinksGUI

% Last Modified by GUIDE v2.5 31-Jul-2020 13:47:17

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditElectricalPowerSinksGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditElectricalPowerSinksGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before lvd_EditElectricalPowerSinksGUI is made visible.
function lvd_EditElectricalPowerSinksGUI_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to lvd_EditElectricalPowerSinksGUI (see VARARGIN)

    % Choose default command line output for lvd_EditElectricalPowerSinksGUI
    handles.output = hObject;

    lvdData = varargin{1};
    setappdata(hObject,'lvdData',lvdData);

    populateGUI(handles, lvdData);

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes lvd_EditElectricalPowerSinksGUI wait for user response (see UIRESUME)
    uiwait(handles.lvd_EditElectricalPowerSinksGUI);

function populateGUI(handles, lvdData)
    [powerSinkListBoxStr, powerSinks] = lvdData.launchVehicle.getPowerSinksListBoxStr();
    set(handles.sinkListBox,'String',powerSinkListBoxStr);
    set(handles.sinkListBox,'Value',1);

    numPowerSinks = length(powerSinks);
    if(numPowerSinks <= 0)
        handles.removeSinkButton.Enable = 'off';
    end

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditElectricalPowerSinksGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = [];
    else
        close(handles.lvd_EditElectricalPowerSinksGUI);
    end

% --- Executes on button press in closeButton.
function closeButton_Callback(hObject, eventdata, handles)
% hObject    handle to closeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    uiresume(handles.lvd_EditElectricalPowerSinksGUI);

% --- Executes on selection change in sinkListBox.
function sinkListBox_Callback(hObject, eventdata, handles)
% hObject    handle to sinkListBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns sinkListBox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from sinkListBox
    if(strcmpi(get(handles.lvd_EditElectricalPowerSinksGUI,'SelectionType'),'open'))
        lvdData = getappdata(handles.lvd_EditElectricalPowerSinksGUI,'lvdData');
        lv = lvdData.launchVehicle;

        selPowerSink = get(handles.sinkListBox,'Value');
        powerSink = lv.getPowerSinkForInd(selPowerSink);

        stageState = getStageStateForPowerState(powerSink, lvdData);
        oldSinkState = stageState.getStateForPowerSink(powerSink);
        stageState.removePowerSinkStateForPowerState(powerSink);

        powerSink.openEditDialog();

        stageState = getStageStateForPowerState(powerSink, lvdData);
        newPowerSinkState = powerSink.createDefaultInitialState(stageState);

        if(not(isempty(oldSinkState)))
            newPowerSinkState.setActiveState(oldSinkState.getActiveState());
        end
        stageState.addPowerSinkState(newPowerSinkState);

        set(handles.sinkListBox,'String',lvdData.launchVehicle.getPowerSinksListBoxStr());
    end

function stageState = getStageStateForPowerState(powerSink, lvdData)
        stage = powerSink.stage;
        stageStates = lvdData.initStateModel.stageStates;
        stageStateInd = find([stageStates.stage] == stage,1,'first');
        stageState = stageStates(stageStateInd);

% --- Executes during object creation, after setting all properties.
function sinkListBox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to sinkListBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in addSinkButton.
function addSinkButton_Callback(hObject, eventdata, handles)
% hObject    handle to addSinkButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditElectricalPowerSinksGUI,'lvdData');

    sinkListBoxStr = LaunchVehiclePowerSinkEnum.getListBoxStr();
    [Selection,ok] = listdlgARH('ListString', sinkListBoxStr, ...
                                'SelectionMode','single', ...
                                'Name', 'New Sink Type', ...
                                'PromptString','Select type of new sink:', ...
                                'ListSize',[300,300]);

	if(ok == 1)
        initStage = lvdData.launchVehicle.getStageForInd(1);

        sinkEnum = LaunchVehiclePowerSinkEnum.getEnumForListboxStr(sinkListBoxStr{Selection});
        switch sinkEnum
            case LaunchVehiclePowerSinkEnum.Simple
                powerSink = LaunchVehicleSimplePwrSink(initStage);

            otherwise
                error('Unknown Power Sink Type: %s', class(sinkEnum))
        end

        useTF = powerSink.openEditDialog();

        if(useTF)
            stageStates = lvdData.initStateModel.stageStates;
            stageStateInd = find([stageStates.stage] == powerSink.stage,1,'first');
            stageState = stageStates(stageStateInd);

            newSinkState = powerSink.createDefaultInitialState(stageState);
            stageState.addPowerSinkState(newSinkState);

            set(handles.sinkListBox,'String',lvdData.launchVehicle.getPowerSinksListBoxStr());

            handles.removeSinkButton.Enable = 'on';
        end
	end

% --- Executes on button press in removeSinkButton.
function removeSinkButton_Callback(hObject, eventdata, handles)
% hObject    handle to removeSinkButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditElectricalPowerSinksGUI,'lvdData');
    lv = lvdData.launchVehicle;

    selPowerSink = get(handles.sinkListBox,'Value');
    powerSink = lv.getPowerSinkForInd(selPowerSink);

    tf = powerSink.isInUse();

    if(tf == false)
        stage = powerSink.getAttachedStage();

        stageStates = lvdData.initStateModel.stageStates;
        stageStateInd = find([stageStates.stage] == stage,1,'first');
        stageState = stageStates(stageStateInd);

        stageState.removePowerSinkStateForPowerState(powerSink);

        stage.removePwrSink(powerSink);

        [listBoxStr, powerSinks] = lvdData.launchVehicle.getPowerSinksListBoxStr();
        set(handles.sinkListBox,'String',listBoxStr);

        numPowerSinks = length(powerSinks);
        if(selPowerSink > numPowerSinks)
            set(handles.sinkListBox,'Value',numPowerSinks);
        end

        if(numPowerSinks <= 0)
            set(handles.sinkListBox,'Value',1);

            handles.removeSinkButton.Enable = 'off';
        end
    else
        warndlg(sprintf('Could not delete the electical power sink "%s" because it is in use as part of an event termination condition, event action, objective function, or constraint.  Remove the dependencies before attempting to delete the electical power sink.', powerSink.getName()),'Cannot Delete Power Sink','modal');
    end

% --------------------------------------------------------------------
function copySinkMenu_Callback(hObject, eventdata, handles)
% hObject    handle to copySinkMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditElectricalPowerSinksGUI,'lvdData');
    lv = lvdData.launchVehicle;

    selPowerSink = get(handles.sinkListBox,'Value');
    powerSink = lv.getPowerSinkForInd(selPowerSink);

    newPowerSink = powerSink.copy();
    initStage = newPowerSink.getAttachedStage();

	initStage.addPwrSink(newPowerSink);

    stageStates = lvdData.initStateModel.stageStates;
    stageStateInd = find([stageStates.stage] == initStage,1,'first');
    stageState = stageStates(stageStateInd);

    newSinkState = powerSink.createDefaultInitialState(stageState);
    stageState.addPowerSinkState(newSinkState);

    set(handles.sinkListBox,'String',lvdData.launchVehicle.getPowerSinksListBoxStr());

    handles.removeSinkButton.Enable = 'on';

% --------------------------------------------------------------------
function sinksListboxContextMenu_Callback(hObject, eventdata, handles)
% hObject    handle to sinksListboxContextMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on key press with focus on lvd_EditElectricalPowerSinksGUI or any of its controls.
function lvd_EditElectricalPowerSinksGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditElectricalPowerSinksGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            uiresume(handles.lvd_EditElectricalPowerSinksGUI);
        case 'enter'
            uiresume(handles.lvd_EditElectricalPowerSinksGUI);
        case 'escape'
            uiresume(handles.lvd_EditElectricalPowerSinksGUI);
    end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function weig = get_w2(sv,method)
% w = get_w(sv)
% High value of w means big decrease in the margin.
% Warning: do not return values of w here. Return a score for each feature
% use the max, change max into sum should not change a lot.
% can be used only if same kernel for all base machines
%% method = 1 ->  take the sum
%% method = 2 ->  used for the l0 update, where the get_w is different from rfe

if nargin < 2,
    method = 1;
end;

if method<2,

    temp = sv.child{1};
    ke=temp.child;
    for i=1:length(sv.child),
        temp = sv.child{i};
            if (temp.child==ke)|strcmp(ke.ker,'weighted_linear'),
              alphaTemp(:,i) = temp.alpha;
            else
                error('Get w for one-vs-rest valid only when the kernel is the same for all base machines.');
            return;
            end;
    end;
    Q = size(alphaTemp,2);
    kertmp = ke.kerparam;
    Ktmp = get_kernel(ke,temp.Xsv,temp.Xsv);
    if strcmp(ke.ker,'linear')|strcmp(ke.ker,'weighted_linear'),
  % compute the multi-class margin for each feature removed.
        if strcmp(ke.ker,'weighted_linear'),
         xTemp = get_x(temp.Xsv);
         xTemp = xTemp.*(ones(size(xTemp,1),1)*ke.kerparam);
         Ktmp = xTemp*xTemp';
        end;
        xTemp = get_x(temp.Xsv);
        temp=xTemp'*alphaTemp;
        tmp2=xTemp'*sum(alphaTemp,2);
        weig = sum(temp.^2,2)' + 1/Q*(tmp2.^2)';
        weig = max(weig)-weig;
        return;
    else
        if strcmp(ke.ker,'rbf'),
      % compute the kernel matrix for all components

         xTemp = get_x(temp.Xsv);

        temp=zeros(1,Q);
        for j=1:Q,
          temp(j)=alphaTemp(:,j)'*(Ktmp)*alphaTemp(:,j);
        end;
        temp = sum(temp) + 1/Q*(sum(alphaTemp,2)'*Ktmp*sum(alphaTemp,2));

      % compute the margin when one component is removed
        for i = 1:size(xTemp,2),
           Ki = xTemp(:,i)*ones(1,size(xTemp,1)) - ones(size(xTemp,1),1)*xTemp(:,i)';
           Ki = Ki.^2;
           Ki = Ki/(2*kertmp^2);
           Ki = exp(Ki);
           for j=1:Q,
            tmpi(j)=alphaTemp(:,j)'*(Ktmp.*Ki)*alphaTemp(:,j);
           end;
            tmpi = sum(tmpi) + 1/Q*(sum(alphaTemp,2)'*(Ktmp.*Ki)*sum(alphaTemp,2));
           weig(i) = temp-tmpi;
        end;
        weig = max(weig)-weig;
       return;
    elseif strcmp(ke.ker,'poly'),
        xTemp = get_x(temp.Xsv);

        temp=zeros(1,Q);
      for j=1:Q,
          temp(j)=alphaTemp(:,j)'*(Ktmp)*alphaTemp(:,j);
      end;
      temp = sum(temp) + 1/Q*(sum(alphaTemp,2)'*Ktmp*sum(alphaTemp,2));
        for i = 1:size(xTemp,2),
           Ki = xTemp(:,i)*xTemp(:,i)';
           K_i = (Ktmp - Ki + 1).^(kertmp);
           for j=1:Q,
            tmpi(j)=alphaTemp(:,j)'*(Ki)*alphaTemp(:,j);
           end;
            tmpi = sum(tmpi) + 1/Q*(sum(alphaTemp,2)'*(Ki)*sum(alphaTemp,2));
           weig(i) = temp-tmpi;
        end;
       weig=max(weig)-weig;
       return;
    end;% if strcmp(...,'rbf')
end;% if strcmp(...,'linear')
%%% Otherwise, all classic kernels have been computed. Now use the get_kernel generically
if strcmp(ke.ker,'custom'),
    error('Get w not implemented for CUSTOM kernels.');
    return;
end;
% compute the multi-class margin for each feature removed.

   xTemp = get_x(temp.Xsv);
   temp=zeros(1,Q);
   for j=1:Q,
          temp(j)=alphaTemp(:,j)'*(Ktmp)*alphaTemp(:,j);
   end;
   temp = sum(temp) + 1/Q*(sum(alphaTemp,2)'*Ktmp*sum(alphaTemp,2));

   for i = 1:size(xTemp,2),
       datTemp = data('temp',[xTemp(:,1:(i-1)),xTemp(:,(i+1):size(xTemp,2))],[]);
       Ki = get_kernel(ke,datTemp,datTemp);
       for j=1:Q,
            tmpi(j)=alphaTemp(:,j)'*(Ki)*alphaTemp(:,j);
       end;
       tmpi = sum(tmpi) + 1/Q*(sum(alphaTemp,2)'*(Ki)*sum(alphaTemp,2));
       weig(i) = temp-tmpi;
   end;
   weig=max(weig)-weig;
return;
else, %% if method < 3
    %% method for the multiplicative update
    for i=1:length(sv.child),
        atmp = sv.child{i};
        W(:,i) = abs(get_w(atmp)');
    end;
    weig = sum(W,2);
    weig=weig';
end;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function I = reconstruction(temp_vector, unit, t1,t2, count1)
I_d_temp = zeros(t1,t2);
for ii=1:t1-unit+1
    for jj=1:t2-unit+1
        temp = temp_vector(:,(ii-1)*(t2-unit+1)+jj);
        I_d_temp(ii:(unit+ii-1), jj:(unit+jj-1)) = I_d_temp(ii:(unit+ii-1), jj:(unit+jj-1)) + reshape(temp, [unit unit]);
    end
end
I = I_d_temp./count1;
end

author: thfylsty


text:
% The COBRAToolbox: testpFBA.m
%
% Purpose:
%     - tests the basic functionality of pFBA
%       Tests the basic solution for both minimizing the flux of gene-
%       associated reactions and all rxns, while growing on gluose or lactose
%       minimal media. Does not test the functionality of the map function.
%
% Authors:
%     - Original file: Nathan Lewis 08/30/10
%     - CI integration: Laurent Heirendt February 2017
%
% Note:
%     - The solver libraries must be included separately

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testpFBA'));
cd(fileDir);

%tolerance
tol = 1e-8;

% load models and expected results
model_glc = readCbModel('testpFBAData.mat','modelName','model_glc');
model_lac = readCbModel('testpFBAData.mat','modelName','model_lac');

objGenes = load('testpFBAData.mat', 'GeneClasses_glc2', 'GeneClasses_glc1', 'GeneClasses_glc0', 'GeneClasses_lac2', 'GeneClasses_lac1', 'GeneClasses_lac0');
objRxns = load('testpFBAData.mat', 'RxnClasses_glc2', 'RxnClasses_glc1', 'RxnClasses_glc0', 'RxnClasses_lac2', 'RxnClasses_lac1', 'RxnClasses_lac0');
objModel = load('testpFBAData.mat', 'modelIrrev_glc2', 'modelIrrev_glc1', 'modelIrrev_glc0', 'modelIrrev_lac2', 'modelIrrev_lac1', 'modelIrrev_lac0');

% list of solver packages
solverPkgs = {'gurobi'}; % 'tomlab_cplex', 'glpk'

% create a parallel pool
try
    minWorkers = 2;
    myCluster = parcluster(parallel.defaultClusterProfile);

    if myCluster.NumWorkers >= minWorkers
        poolobj = gcp('nocreate');  % if no pool, do not create new one.
        if isempty(poolobj)
            parpool(minWorkers);  % launch minWorkers workers
        end
    end
catch
    disp('Trying non parallel test')
end

for k = 1:length(solverPkgs)

    fprintf(' -- Running testfindBlockedReaction using the solver interface: %s ... \n', solverPkgs{k});

    solverLPOK = changeCobraSolver(solverPkgs{k}, 'LP', 0);

    if solverLPOK

        % run pFBA
        fprintf('\n*** Test basic pFBA calculations ***\n\n');
        fprintf('\n** Optimal solution - minimize flux not associated with gene: glucose\n');
        [t_objGenes.GeneClasses_glc1 t_objRxns.RxnClasses_glc1 t_objModel.modelIrrev_glc1] = pFBA(model_glc, 'geneoption', 2);

        fprintf('\n** Optimal solution - minimize flux not associated with gene: lactate\n');
        [t_objGenes.GeneClasses_lac1 t_objRxns.RxnClasses_lac1 t_objModel.modelIrrev_lac1] = pFBA(model_lac, 'geneoption', 2);

        fprintf('\n** Optimal solution - minimize gene-associated flux: glucose\n');
        [t_objGenes.GeneClasses_glc1 t_objRxns.RxnClasses_glc1 t_objModel.modelIrrev_glc1] = pFBA(model_glc, 'geneoption', 1);

        fprintf('\n** Optimal solution - minimize gene-associated flux: lactate\n');
        [t_objGenes.GeneClasses_lac1 t_objRxns.RxnClasses_lac1 t_objModel.modelIrrev_lac1] = pFBA(model_lac, 'geneoption', 1);

        fprintf('\n** Optimal solution - minimize all flux: glucose **\n');
        [t_objGenes.GeneClasses_glc0 t_objRxns.RxnClasses_glc0 t_objModel.modelIrrev_glc0] = pFBA(model_glc, 'geneoption', 0);

        fprintf('\n** Optimal solution - minimize all flux: lactate **\n');
        [t_objGenes.GeneClasses_lac0 t_objRxns.RxnClasses_lac0 t_objModel.modelIrrev_lac0] = pFBA(model_lac, 'geneoption', 0);

        t_objGenesf = fieldnames(t_objGenes);
        t_objRxnsf = fieldnames(t_objRxns);
        t_objModelf = fieldnames(t_objModel);

        % testing if gene lists are consistent with expected lists
        t_fg = zeros(40, 1);
        cnt = 0;
        for i = 1:length(t_objGenesf)
            tmp_lists = fieldnames(t_objGenes.(t_objGenesf{i}));
            for j = 1:length(tmp_lists)
                t1 = find(~ismember(t_objGenes.(t_objGenesf{i}).(tmp_lists{j}), objGenes.(t_objGenesf{i}).(tmp_lists{j})));
                t2 = find(~ismember(objGenes.(t_objGenesf{i}).(tmp_lists{j}), t_objGenes.(t_objGenesf{i}).(tmp_lists{j})));
                cnt = cnt + 1;
                if isempty(t1)
                    t_fg(cnt) = 1;
                end
                cnt = cnt + 1;
                if isempty(t2)
                    t_fg(cnt) = 1;
                end
            end
        end

        assert(min(t_fg) == 1)

        % testing if rxn lists are consistent with expected lists
        t_fr = zeros(40, 1);
        cnt = 0;
        for i = 1:length(t_objRxnsf)
            tmp_lists = fieldnames(t_objRxns.(t_objRxnsf{i}));
            for j = 1:length(tmp_lists)
                t1 = find(~ismember(t_objRxns.(t_objRxnsf{i}).(tmp_lists{j}), objRxns.(t_objRxnsf{i}).(tmp_lists{j})));
                t2 = find(~ismember(objRxns.(t_objRxnsf{i}).(tmp_lists{j}), t_objRxns.(t_objRxnsf{i}).(tmp_lists{j})));
                cnt = cnt + 1;
                if isempty(t1)
                    t_fr(cnt) = 1;
                end
                cnt = cnt + 1;
                if isempty(t2)
                    t_fr(cnt) = 1;
                end
            end
        end

        assert(min(t_fr) == 1)

        % testing if flux minima are consistent with expected values
        t_fm = zeros(8, 1);
        cnt = 0;
        for i = 1:length(t_objModelf)
            t1 = t_objModel.(t_objModelf{i}).lb(findRxnIDs(t_objModel.(t_objModelf{i}), 'netFlux')) - objModel.(t_objModelf{i}).lb(findRxnIDs(objModel.(t_objModelf{i}), 'netFlux'));
            t2 = t_objModel.(t_objModelf{i}).ub(findRxnIDs(t_objModel.(t_objModelf{i}), 'netFlux')) - objModel.(t_objModelf{i}).ub(findRxnIDs(objModel.(t_objModelf{i}), 'netFlux'));
            cnt = cnt + 1;
            if t1 < tol
                t_fm(cnt) = 1;
            end
            cnt = cnt + 1;
            if t2 < tol
                t_fm(cnt) = 1;
            end
        end

        assert(min(t_fm) == 1)

        x = min([t_fm; t_fg; t_fr]);

        % output a success message
        fprintf('Done.\n');
    end
end

% change back to the root directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function [xPred, SPred, xPropCenPoints]=sqrtDiscCubKalPred(xPrev,SPrev,f,SQ,xi,w,stateDiffTrans,stateAvgFun,stateTrans)
%%SQRTDISCCUBKALPRED Perform the discrete-time prediction step that comes
%                    with the square root implementation of the cubature
%                    Kalman filter with additive process noise.
%
%INPUTS: xPrev The xDim X 1 state estimate at the previous time-step.
%        SPrev The xDim X xDim lower-triangular square root of the state
%              covariance matrix at the previous time-step.
%            f A function handle for the state transition function that
%              takes the state as its parameter.
%           SQ The xDimXxDim lower-triangular square root of the process
%              noise covariance matrix.
%           xi An xDimXnumCubPoints matrix of cubature points. If this and
%              the next parameter are omitted or empty matrices are passed,
%              then fifthOrderCubPoints(xDim+cDim) is used. It is suggested
%              that xi and w be provided to avoid needless recomputation of
%              the cubature points.
%            w A numCubPointsX1 vector of the weights associated with the
%              cubature points. These must be all positive.
% stateDiffTrans An optional function handle that takes an xDimXN matrix of
%              N differences between states estimates and transforms them
%              however might be necessary. For example, a state continaing
%              angular components will generally need differences between
%              angular components wrapped to the range +/-pi.
%  stateAvgFun An optional function that given an xDimXN matrix of N state
%              estimates and an NX1 vector of weights, provides the
%              weighted average of the state estimates. This is necessary
%              if, for example, states with angular components are
%              averaged.
%   stateTrans An optional function that takes a state estimate and
%              transforms it. This is useful if one wishes the elements of
%              the state to be bound to a certain domain. For example, if
%              an element of the state is an angle, one should generally
%              want to bind it to the region +/-pi. This is not applied to
%              the output of f.
%
%OUTPUTS: xPred The xDim X 1 predicted state estimate.
%         SPred The xDim X xDim lower-triangular square root of the
%               predicted state covariance estimate.
% xPropCenPoints The centered propagated cubature state points. This matrix
%               is needed if the backwards smoothing step is being applied.
%
%The mathematics behind the function sqrtDiscCubKalPred are described in
%more detail in Section IX of [1] and in [2]. The use of stateDiffTrans,
%stateAvgFun, and stateTrans arises when dealing with state space models
%involving, for example, angles. This is analogous to similar issues that
%arise when dealing with angular measurements as described in [3].
%
%REFERENCES:
%[1] David F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems Magazine,
%    vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%[2] I. Arasaratnam and S. Haykin, "Cubature Kalman filters," IEEE
%    Transactions on Automatic Control, vol. 54, no. 6, pp. 1254-1269,
%    Jun. 2009.
%[3] D. F. Crouse, "Cubature/ unscented/ sigma point Kalman filtering with
%    angular measurement models," in Proceedings of the 18th International
%    Conference on Information Fusion, Washington, D.C., 6-9 Jul. 2015.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    xDim=length(xPrev);

    if(nargin<5||isempty(xi))
        [xi,w]=fifthOrderCubPoints(xDim);
    end

    if(nargin<7||isempty(stateDiffTrans))
        stateDiffTrans=@(x)x;
    end

    if(nargin<8||isempty(stateAvgFun))
        stateAvgFun=@(x,w)calcMixtureMoments(x,w);
    end

    if(nargin<9||isempty(stateTrans))
        stateTrans=@(x)x;
    end

    xDim=length(xPrev);
    numCubPoints=length(w);

    xPropPoints=zeros(xDim,numCubPoints);%Allocate space

    %Calculate the cubature state points
    xPoints=stateTrans(transformCubPoints(xi,xPrev,SPrev));

    %Propagate the cubature state points
    for curP=1:numCubPoints
        xPropPoints(:,curP)=f(xPoints(:,curP));
    end

    %Calculate the predicted state
    xPred=stateAvgFun(xPropPoints,w);

    %Centered, propagated cubature state points
    xDiff=stateDiffTrans(bsxfun(@minus,xPropPoints,xPred));
    xPropCenPoints=bsxfun(@times,xDiff,sqrt(w'));

    %The root prediction covariance.
    SPred=tria([xPropCenPoints, SQ]);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function trans = createRotationOx(varargin)
%CREATEROTATIONOX Create the 4x4 matrix of a 3D rotation around x-axis.
%
%   TRANS = createRotationOx(THETA);
%   Returns the transform matrix corresponding to a rotation by the angle
%   THETA (in radians) around the Ox axis. A rotation by an angle of PI/2
%   would transform the vector [0 1 0] into the vector [0 0 1].
%
%   The returned matrix has the form:
%   [1      0            0      0]
%   [0  cos(THETA) -sin(THETA)  0]
%   [0  sin(THETA)  cos(THETA)  0]
%   [0      0            0      1]
%
%   TRANS = createRotationOx(ORIGIN, THETA);
%   TRANS = createRotationOx(X0, Y0, Z0, THETA);
%   Also specifies origin of rotation. The result is similar as performing
%   translation(-X0, -Y0, -Z0), rotation, and translation(X0, Y0, Z0).
%
%   See also
%   transforms3d, transformPoint3d, createRotationOy, createRotationOz
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-02-18
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

% default values
dx = 0;
dy = 0;
dz = 0;
theta = 0;

% get input values
if length(varargin) == 1
    % only one argument -> rotation angle
    theta = varargin{1};

elseif length(varargin) == 2
    % origin point (as array) and angle
    var = varargin{1};
    dx = var(1);
    dy = var(2);
    dz = var(3);
    theta = varargin{2};

elseif length(varargin) == 4
    % origin (x and y) and angle
    dx = varargin{1};
    dy = varargin{2};
    dz = varargin{3};
    theta = varargin{4};
end

% compute coefs
cot = cos(theta);
sit = sin(theta);

% create transformation
trans = [...
    1 0 0 0;...
    0 cot -sit 0;...
    0 sit cot 0;...
    0 0 0 1];

% add the translation part
t = [1 0 0 dx;0 1 0 dy;0 0 1 dz;0 0 0 1];
trans = t * trans / t;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function xyl_header_write ( output_filename, output_unit, point_num, ...
  line_num, line_data_num )

%*****************************************************************************80
%
%% XYL_HEADER_WRITE writes the header of an XYL file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the output file.
%
%    Input, integer OUTPUT_UNIT, the output file unit number.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer LINE_NUM, the number of lines.
%
%    Input, integer LINE_DATA_NUM, the number of line items.
%
  fprintf ( output_unit, '#  %s\n', output_filename );
  fprintf ( output_unit, '#  created by xy_io::xyl_header_write.m\n' );
  fprintf ( output_unit, '#\n');
  fprintf ( output_unit, '#  Number of points =     %d\n', point_num );
  fprintf ( output_unit, '#  Number of lines =      %d\n', line_num );
  fprintf ( output_unit, '#  Number of line items = %d\n', line_data_num );
  fprintf ( output_unit, '#\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function s1 = resize_scales(s0,dim,args)
% Resize scalefactors
% _________________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

dim = [dim ones(1,max(numel(args)-numel(dim),0))];
args1 = cell(1,numel(args));
for i=1:numel(args),
    if max(args{i})>dim(i) || min(args{i})<1,
        error('Index exceeds matrix dimensions (1).');
    end;

    if size(s0,i)==1,
        args1{i} = ones(size(args{i}));
    else
        args1{i} = args{i};
    end;
end;

s1 = s0(args1{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%  Figure 7.13      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% fig7_13.m is a script to generate Figure 7.13

% Impulse response
clf;
f=[0 1;-1 0];
g=[0;1];
h=[1 0];
K=[3 4];
fc=f-g*K;
gc=[4*g [1;0]];
hc=[h;[0 1]; -K/4];
jc=[0 0;0 0; 1 0];
sysCL=ss(fc,gc,hc,jc);
t=0:0.1:7;
y=impulse(sysCL,t);
plot(t,y(:,:,2));
xlabel('Time (sec)');
ylabel('Amplitude');
text(1.7,-.3,'x_2');
text(.5,.1,'u/4');
text(.7,.7,'x_1');
title('Fig. 7.13 Impulse response of the oscillator with full state feedback');
nicegrid;


text:
function c = minus(a,b)
	checkinputs(a,b);
	c = uint64minus(uint64(a),uint64(b));
end


text:
function s = setTerminal()
% This is needed on a MAC to execute some shell commands
% Otherwise the system command executes its own shell that forces the user
% to press Enter key to continue matlab script execution.
s = '';
if ismac()
    s = 'TERM=ansi; ';
end

author: BUNPC


function hashTable = loadHashTable(suite)
    % loads a reference hash table from disk
    hashTable.suite = suite;
    hashTable.contents = struct();
    filename = hashTableName(suite);
    if exist(filename, 'file')
        fid = fopen(filename, 'r');
        finally_fclose_fid = onCleanup(@() fclose(fid));

        data = textscan(fid, '%s : %s');
        if ~isempty(data) && ~all(cellfun(@isempty, data))
            functions = cellfun(@strtrim, data{1},'UniformOutput', false);
            hashes    = cellfun(@strtrim, data{2},'UniformOutput', false);
            for iFunc = 1:numel(functions)
                hashTable.contents.(functions{iFunc}) = hashes{iFunc};
            end
        end
    end
end


text:
function stroud_test22 ( )

%*****************************************************************************80
%
%% TEST22 tests OCTAHEDRON_UNIT_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_ND_INDEX;

  n_max = 3;

  num = function_2d_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST22\n' );
  fprintf ( 1, '  OCTAHEDRON_UNIT_ND approximates integrals in a unit\n' );
  fprintf ( 1, '    octahedron in N dimensions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      F(X)       N = 1         N = 2         N = 3\n' );
  fprintf ( 1, '\n' );

  for i = 1 : num

    FUNC_ND_INDEX = i;

    for n = 1 : n_max
      result(n) = octahedron_unit_nd ( 'function_nd', n );
    end

    fname = function_nd_name ( i );

    fprintf ( 1, '  %s', fname );
    for n = 1 : n_max
      fprintf ( 1, '  %12f', result(n) );
    end
    fprintf ( 1, '\n' );

  end

  return
end


text:
function y = norm_lpdf(x,mu,sigma)
%NORM_LPDF Normal log-probability density function (lpdf).
%
%   Y = NORM_LPDF(X,MU,SIGMA) Returns the log of the normal pdf with
%   mean, MU, and standard deviation, SIGMA, at the values in X.
%
%   The size of Y is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.
%
%   Default values for MU and SIGMA are 0 and 1 respectively.

% Copyright (c) 1998-2004 Aki Vehtari

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

if nargin < 3,
  sigma = 1;
end

if nargin < 2;
  mu = 0;
end

if nargin < 1,
  error('Requires at least one input argument.');
end

y = -0.5 * ((x-mu)./sigma).^2 -log(sigma) -log(2*pi)/2;


text:
function value = mono_000_3d ( n, x )

%*****************************************************************************80
%
%% MONO_000_3D evaluates X**0 Y**0 Z**0.
%
%  Modified:
%
%    25 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, real X(N), the point at which the monomial is to be evaluated.
%
%    Output, real VALUE, the value of the monomial.
%
  value = 1.0E+00;

  return
end


text:
function varargout = process_evt_remove_simult( varargin )
% PROCESS_EVT_REMOVE_SIMULT: Remove occurences of an event A occurring at the same time as an event B.
%
% USAGE:  OutputFiles = process_evt_remove_simult('Run', sProcess, sInputs)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2014-2019

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Remove simultaneous';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Events';
    sProcess.Index       = 49;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/ArtifactsDetect#Remove_simultaneous_blinks.2Fheartbeats';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'raw'};
    sProcess.OutputTypes = {'data', 'raw'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Event name: to remove
    sProcess.options.remove.Comment = 'Remove events named: ';
    sProcess.options.remove.Type    = 'text';
    sProcess.options.remove.Value   = 'cardiac';
    % Event name: target
    sProcess.options.target.Comment = 'When too close to events: ';
    sProcess.options.target.Type    = 'text';
    sProcess.options.target.Value   = 'blink';
    % Maximum duration events
    sProcess.options.dt.Comment = 'Minimum delay between events: ';
    sProcess.options.dt.Type    = 'value';
    sProcess.options.dt.Value   = {0.250, 'ms', 0};
    % Delete original events
    sProcess.options.rename.Comment = 'Rename the events instead of deleting them';
    sProcess.options.rename.Type    = 'checkbox';
    sProcess.options.rename.Value   = 0;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Return all the input files
    OutputFiles = {};

    % ===== GET OPTIONS =====
    % Get the options
    isDelete = ~sProcess.options.rename.Value;
    dt = sProcess.options.dt.Value{1};
    evtA = strtrim(sProcess.options.remove.Value);
    evtB = strtrim(sProcess.options.target.Value);
    % Check that the user gave one event in each box
    if isempty(evtA) || isempty(evtB)
        bst_report('Error', sProcess, [], 'Event name is empty.');
        return;
    elseif any(evtA == ',') || any(evtA == ';')  || any(evtB == ',') || any(evtB == ';')
        bst_report('Error', sProcess, [], 'This event can be used to process only one event type at a time.');
        return;
    end

    % ===== PROCESS ALL FILES =====
    % For each file
    for iFile = 1:length(sInputs)
        % ===== GET FILE DESCRIPTOR =====
        % Load the raw file descriptor
        isRaw = strcmpi(sInputs(iFile).FileType, 'raw');
        if isRaw
            DataMat = in_bst_data(sInputs(iFile).FileName, 'F');
            sEvents = DataMat.F.events;
        else
            DataMat = in_bst_data(sInputs(iFile).FileName, 'Events');
            sEvents = DataMat.Events;
        end
        % If no markers are present in this file
        if isempty(sEvents)
            bst_report('Error', sProcess, sInputs(iFile), 'This file does not contain any event. Skipping File...');
            continue;
        end
        % Call the detection function
        [sEvents, isModified] = Compute(sInputs(iFile), sEvents, evtA, evtB, dt, isDelete);

        % ===== SAVE RESULT =====
        % Only save changes if something was change
        if isModified
            % Report changes in .mat structure
            if isRaw
                DataMat.F.events = sEvents;
            else
                DataMat.Events = sEvents;
            end
            % Save file definition
            bst_save(file_fullpath(sInputs(iFile).FileName), DataMat, 'v6', 1);
        end
        % Return all the input files
        OutputFiles{end+1} = sInputs(iFile).FileName;
    end
end

%% ===== RENAME SIMULTANEOUS EVENTS =====
function [eventsNew, isModified] = Compute(sInput, events, evtA, evtB, dt, isDelete)
    % Initialize returned variables
    isModified = 0;
    eventsNew = events;
    % Get events names
    iEvtA = find(strcmpi({events.label}, evtA));
    iEvtB = find(strcmpi({events.label}, evtB));
    if isempty(iEvtA)
        bst_report('Error', 'process_evt_remove_simult', sInput, ['Event "' evtA '" does not exist.']);
        return;
    elseif isempty(iEvtB)
        bst_report('Error', 'process_evt_remove_simult', sInput, ['Event "' evtB '" does not exist.']);
        return;
    end
    % If events are extended events: error
    if (size(events(iEvtA).times,1) > 1) || (size(events(iEvtB).times,1) > 1)
        bst_report('Error', 'process_evt_remove_simult', sInput, 'Cannot process extended events.');
        return;
    end
    % Look of all the events A that are too close to event B
    iRemoveA = [];
    for iOcc = 1:size(events(iEvtA).times,2)
        if any(abs(events(iEvtA).times(iOcc) - events(iEvtB).times) <= dt)
            iRemoveA(end+1) = iOcc;
        end
    end
    % No events to remove
    if isempty(iRemoveA)
        return;
    end
    % Display info message with the number of events removed
    bst_report('Info', 'process_evt_remove_simult', sInput, sprintf('Removed %d events "%s" that were less than %d ms away from an event "%s".', length(iRemoveA), evtA, round(dt*1000), evtB));
    % If we need to rename the events to a new category
    if ~isDelete
        % Create new event name for the removed occurrences
        newLabel = file_unique([evtA '_rm'], {events.label});
        % Create new event
        iEvtRm = length(eventsNew) + 1;
        eventsNew(iEvtRm) = db_template('event');
        eventsNew(iEvtRm).label    = newLabel;
        eventsNew(iEvtRm).color    = [1 0 0];
        eventsNew(iEvtRm).times    = eventsNew(iEvtA).times(:,iRemoveA);
        eventsNew(iEvtRm).epochs   = eventsNew(iEvtA).epochs(iRemoveA);
        if ~isempty(eventsNew(iEvtA).channels)
            eventsNew(iEvtRm).channels = eventsNew(iEvtA).channels(iRemoveA);
        end
        if ~isempty(eventsNew(iEvtA).notes)
            eventsNew(iEvtRm).notes = eventsNew(iEvtA).notes(iRemoveA);
        end
    end
    % Remove occurrences / remove event
    if isequal(iRemoveA, 1:size(eventsNew(iEvtA).times,2))
        eventsNew(iEvtA) = [];
    else
        eventsNew(iEvtA).times(:,iRemoveA)  = [];
        eventsNew(iEvtA).epochs(iRemoveA)   = [];
        if ~isempty(eventsNew(iEvtA).channels)
            eventsNew(iEvtA).channels(iRemoveA) = [];
        end
        if ~isempty(eventsNew(iEvtA).notes)
            eventsNew(iEvtA).notes(iRemoveA) = [];
        end
    end
    isModified = 1;
end

author: brainstorm-tools


text:
function [A] = dicomrt_rotate180y(matrix)
% dicomrt_rotate180y(matrix)
%
% Rotate 3D matrix accordingly with original patient position.
%
% Example:
%
% [A]=dicomrt_rotate180y(Dose)
%
% rotate Dose 3D matrice of 180 degrees about the Y axis
%
% See also dicomrt_ctcreate, dicomrt_loaddose, dicomrt_rotate180x
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check cases
if ndims(matrix) ~= 3
    error ('Only 3D matrices are supported. Exit now !');
end

% Execute rotation
for i=1:size(matrix,1)
    temp=squeeze(matrix(i,:,:));
    A(i,:,:)=rot90(temp,2);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function alpha = calculateNoiseType(readings)
%this function calculates the dominate noise type at a given tau value.
%variable d is to track number of difference oppperations for identifying
%noise types **************. d should start at 0 and go no higher than 2
%(not sure purpose of d at this point).

%the following is the algorithm for calculating noise type
done = 0;
d = 0;
z = readings;

while done == 0

    r1 = rOneCalculation(z);

    delta = r1 / (r1 + 1);

    if d >= 0 && (delta < .25 || d >= 2)
        p = -2*(delta+d);
        alpha = p + 2; %for time data
        done = 1;
    else
        z = calculateDiffArray(z);
        d = d + 1;
    end
end


text:
function Z = minus(X,Y)
%MINUS Binary subtraction (-) for tenmat.
%
%   See also TENMAT, TENMAT/TENMATFUN.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

fun = @minus;

% One argument is a scalar
if ((prod(size(X)) == 1 || prod(size(Y)) == 1))
  if (prod(size(Y)) == 1) && isa(X,'tenmat')
    Z = X;
    Z.data = fun(Z.data,Y);
  else
    Z = Y;
    Z.data = fun(X,Z.data);
  end
  return;
end

% Both arguments are tenmats
Z = tenmat(Y);
if ~(isequal(size(Y),size(Z)))
  error('Tenmat size mismatch.')
end
Z.data = fun(X.data,Z.data);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function SRout = SR(VAR,SIGN,VARopt)
% =======================================================================
% Compute IRs, VDs, and HDs for a VAR model estimated with VARmodel and
% identified with sign restrictions
% =======================================================================
% SRout = SR(VAR,R,VARopt)
% -----------------------------------------------------------------------
% INPUT
%   - VAR: structure, result of VARmodel function
%   - SIGN: matrix containing the sign restrictions (nvar,nshocks). 1
%       stands for positive, -1 for negative, 0 for unrestricted. For
%       example the following identifies only one shock with positive
%       impact on the VAR1, and negative imnpact on VAR2 and VAR3
%
%             shock1      shock2     shock3
%  SIGN = [     1           0           0          % VAR1
%              -1           0           0          % VAR2
%              -1           0           0];        % VAR3
%
%   - VARopt: options of the VAR (see VARoption from VARmodel)
% -----------------------------------------------------------------------
% OUTPUT
%   - SRout
%       * IRall : 4-D matrix of IRs  (nsteps,nvar,nshocks,ndraws)
%       * IRmed : median of IRall
%       * IRinf : lower bound of IRall
%       * IRsup : upper bound of IRall
%       * IR    : IR based on true B matrix that is closest to median B
%       * (similar structure for VD, B, and HD)
% -----------------------------------------------------------------------
% EXAMPLE
%   - See VARToolbox_Code.m in "../Primer/"
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

%% Check inputs
%==========================================================================
if ~exist('SIGN','var')
    error('You have not provided sign restrictions (SIGN)')
end

if ~exist('VARopt','var')
    error('You need to provide VAR options (VARopt from VARmodel)');
end

%% Retrieve parameters and preallocate variables
%==========================================================================
nvar    = VAR.nvar;
nvar_ex = VAR.nvar_ex;
nsteps  = VARopt.nsteps;
ndraws  = VARopt.ndraws;
nobs    = VAR.nobs;
nlag    = VAR.nlag;
pctg    = VARopt.pctg;

% Initialize empty matrix for the IR draws
IRall  = nan(nsteps,nvar,nvar,ndraws);
VDall = nan(nsteps,nvar,nvar,ndraws);
HDall.shock  = zeros(nobs+nlag,nvar,nvar,ndraws);
HDall.init   = zeros(nobs+nlag,nvar,ndraws);
HDall.const  = zeros(nobs+nlag,nvar,ndraws);
HDall.trend  = zeros(nobs+nlag,nvar,ndraws);
HDall.trend2 = zeros(nobs+nlag,nvar,ndraws);
HDall.endo   = zeros(nobs+nlag,nvar,ndraws);
HDall.exo    = zeros(nobs+nlag,nvar,nvar_ex,ndraws);
Ball = nan(nvar,nvar,ndraws);

%% Sign restriction routine
%==========================================================================
jj = 0; % accepted draws
tt = 0; % total draws
ww = 1; % index for printing on screen
while jj < ndraws

    % Check total number of draws
    if tt>VARopt.sr_draw
        disp('------------------------------------------------------------')
        disp( 'Total number of draws for finding sign restrictions exceeed')
        disp( 'Change the restrictions or increase VARopt.sr_draw');
        disp('------------------------------------------------------------')
        error('See details above')
    end

    % Set up VAR_draw.(draw{j}) for rotations: Only identification uncertainty
    label  = {['draw' num2str(jj)]};
    VAR_draw.(label{1}) = VAR;
    VARopt.ident = 'sign';
    % If selected, set up VAR_draw.(draw{j}) to consider identification + model uncertainty
    if VARopt.sr_mod==1
        % Draw F and sigma from the posterior and
        [sigma_draw, Ft_draw, F_draw, Fcomp_draw] = VARdrawpost(VAR);
        VAR_draw.(label{1}).Ft = Ft_draw;
        VAR_draw.(label{1}).F = F_draw;
        VAR_draw.(label{1}).Fcomp = Fcomp_draw;
        VAR_draw.(label{1}).sigma = sigma_draw;
    end

    % Compute rotated B matrix
    B = SignRestrictions(SIGN,VAR_draw.(label{1}),VARopt);

    % Note: e = (inv(B)*VAR_draw.(draw{j}).resid')';
    % Check orthogonality:
    % round(corr((inv(B)*VAR_draw.(draw{j}).resid')'))

    if ~isempty(B)
        % Store B
        jj = jj+1; tt = tt+1;
        Ball(:,:,jj) = B;

        % Update VAR_draw.(draw{j}) with the rotated B matrix for IR, VD, and HD
        VAR_draw.(label{1}).B = B;

        % Compute and store IR, VD, HD
        [aux_irf, VAR_draw.(label{1})] = VARir(VAR_draw.(label{1}),VARopt);
        IRall(:,:,:,jj)  = aux_irf;
        aux_fevd = VARvd(VAR_draw.(label{1}),VARopt);
        VDall(:,:,:,jj)  = aux_fevd;
        aux_hd = VARhd(VAR_draw.(label{1}),VARopt);
        HDall.shock(:,:,:,jj) = aux_hd.shock;
        HDall.init(:,:,jj)    = aux_hd.init;
        HDall.const(:,:,jj)   = aux_hd.const;
        HDall.trend(:,:,jj)   = aux_hd.trend;
        HDall.trend2(:,:,jj)  = aux_hd.trend2;
        HDall.endo(:,:,jj)    = aux_hd.endo;
        if nvar_ex>0; HDall.exo(:,:,:,jj) = aux_hd.exo; end

        % Display number of loops
        if jj==VARopt.mult*ww
            disp(['Rotation: ' num2str(jj) ' / ' num2str(ndraws)]);
            ww=ww+1;
        end
    else
        tt=tt+1;
    end
end
disp('-- Done!');
disp(' ');

%% Store results
%==========================================================================
% Store all accepted IRs and VDs
SRout.IRall = IRall;
SRout.VDall = VDall;
SRout.Ball  = Ball;
SRout.HDall = HDall;

% Compute and save median B matrix
SRout.Bmed = median(Ball,3);

% Compute and save true B matrix that is closer to median
aux = sum(sum((Ball-SRout.Bmed).^2,1),2);
sel = find(aux==min(aux));
sel = sel(1);
SRout.B = Ball(:,:,sel);

% Compute IR and VD based on all rotations
pctg_inf = (100-pctg)/2;
pctg_sup = 100 - (100-pctg)/2;
SRout.IRmed = median(IRall,4);
aux = prctile(IRall,[pctg_inf pctg_sup],4);
SRout.IRinf = aux(:,:,:,1);
SRout.IRsup = aux(:,:,:,2);
SRout.VDmed = median(VDall,4);
aux = prctile(VDall,[pctg_inf pctg_sup],4);
SRout.VDinf = aux(:,:,:,1);
SRout.VDsup = aux(:,:,:,2);

% Compute IR, VD, and HD based on the VARdraw that is closest to the median
% B matrix
VARopt.ident = 'sign';
label = {['draw' num2str(sel)]};  % Recover the position in VARdraw
VAR = VAR_draw.(label{1});        % Set a VAR based on the selected VARdraw
SRout.IR = VARir(VAR,VARopt);     % Compute IR
SRout.VD = VARvd(VAR,VARopt);     % Compute VD
SRout.HD = VARhd(VAR,VARopt);     % Compute HD

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function [vacc, hacc, vcm, hcm] = mcmcProcessResult_pixels(imsegs, pg)
% same as mcmcProcessResult except that pg is a pixel confidence map

vcm = zeros(3);
hcm = zeros(5);
vtotal = 0;
htotal = 0;
vacc = 0;
hacc = 0;
for f = 1:numel(pg)

    [imh, imw, nc] = size(pg{f});
    pg{f} = reshape(pg{f}, [imh*imw nc]);

    pgv = [pg{f}(:, 1) sum(pg{f}(:, 2:6), 2) pg{f}(:, 7)];
    pgh = pg{f}(:, 2:6);

    vlab = imsegs(f).vert_labels(imsegs(f).segimage(:));
    hlab = imsegs(f).horz_labels(imsegs(f).segimage(:));
    npix = numel(vlab);
    vlab = vlab(:);
    hlab = hlab(:);

    [maxval, vmax] = max(pgv, [], 2);
    [maxval, hmax] = max(pgh, [], 2);

    vacc = vacc + sum(vlab==vmax)/npix;
    vtotal = vtotal + sum(vlab>0)/npix;

    hacc = hacc + sum((hlab==hmax) & (vlab==2))/npix;
    htotal = htotal + sum((hlab>0) & (vlab==2))/npix;

    for s = 1:numel(vmax)
        if vlab(s)~=0
            vcm(vlab(s), vmax(s)) = vcm(vlab(s), vmax(s)) + 1/npix;
        end
        if hlab(s)~=0
            hcm(hlab(s), hmax(s)) = hcm(hlab(s), hmax(s)) + 1/npix;
        end
    end
end

vacc = vacc / vtotal;
hacc = hacc / max(htotal, 1E-9);

vcm = vcm ./ max(repmat(sum(vcm, 2), 1, size(vcm, 2)), 1E-9);
hcm = hcm ./ max(repmat(sum(hcm, 2), 1, size(hcm, 2)), 1E-9);

%disp(['vacc: ' num2str(vacc) '   hacc: ' num2str(hacc)])
%disp(num2str(sum(ny, 1) / sum(totalv)))

author: zouchuhang

repo: LayoutNet


text:
%VL_NNBNORM CNN batch normalisation.
%   Y = VL_NNBNORM(X,G,B) applies batch normalization to the input
%   X. Batch normalization is defined as:
%
%      Y(i,j,k,t) = G(k) * (X(i,j,k,t) - mu(k)) / sigma(k) + B(k)
%
%   where:
%
%      mu(k) = mean_ijt X(i,j,k,t),
%      sigma2(k) = mean_ijt (X(i,j,k,t) - mu(k))^2,
%      sigma(k) = sqrt(sigma2(k) + EPSILON)
%
%   are respectively the per-channel mean, variance, and standard
%   deviation of each feature channel in the data X. The parameters
%   G(k) and B(k) are multiplicative and additive constants use to
%   scale each data channel.
%
%   Means and variances are accumulated across all the data items
%   (images) stored in the 4D tensor X (from which the name batch
%   normalization is derived). The constant EPSILON is used to
%   regularize the computation of sigma(k) and to avoid division by
%   zero.
%
%   [DZDX,DZDG,DZDB] = VL_NNBNORM(X,G,B,DZDY) computes the derviatives
%   of the block projected onto DZDY. DZDX, DZDG, DZDB and DZDY have
%   the same dimensions as X, G, B, and Y respectivey.
%
%   Optionally, [Y,MOMENTS] = VL_NNBNORM(...) and
%   [DZDX,DZDG,DZDB,MOMENTS] = VL_NNBNORM(...,DZDY) return the values
%   of the vectors mu and sigma in the formulas above. Here, MOMENTS
%   is a DEPTH x 2 array [MU, SIGMA].
%
%   VL_NNBNROM(..., 'Option', value) takes the following options:
%
%   `Epsilon`:: 1e-4
%       Specifies the constant EPSILON in the formuals above.
%
%   `Moments`:: unspecified
%       Specifies an array MOMENTS with the values of mu and sigma to
%       use instead of computing them according to the equations
%       above. This is useful to disable batch normalization during
%       testing.
%
%   See also: VL_NNNORMALIZE().

% Copyright (C) 2015 Sébastien Ehrhardt, Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: jiangqy


text:
 hold on;
 len = length(obps);
 for i = 1 : len
     plot3(obps(i,1), obps(i,2), obps(i,3), '*', 'Color', 'r');
 end

 hold off;

author: LenaShengzhen

repo: AerialRobotics


text:
function net = compact(net)

% COMPACT
%
% Remove duplicate support vectors, adjusting Lagrange multipliers of
% remaining support vectors to compensate.
%
% net = compact(net);

%
% File        : @svc/compact.m
%
% Date        : Tuesday 12th September 2000
%
% Author      : Dr Gavin C. Cawley
%
% Description : Remove duplicate support vectors, adjusting Lagrange
%               multipliers of remaining support vectors to compensate.
%               Part of an object-oriented implementation of Vapnik's Support
%               Vector Machine, as described in [1].
%
% References  : [1] V.N. Vapnik,
%                   "The Nature of Statistical Learning Theory",
%                   Springer-Verlag, New York, ISBN 0-387-94559-8,
%                   1995.
%
% History     : 07/07/2000 - v1.00
%               12/09/2000 - v1.01 minor improvements to comments and help
%                                  messages
%               13/09/2000 - v1.10 zeta (pattern replication factors) and C
%                                  fields removed from svc objects
%
% Copyright   : (c) Dr Gavin C. Cawley, September 2000
%
%    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 should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

[sv, i, j] = unique(net.sv, 'rows');

w = zeros(1, length(i));

for i=1:length(j)

   w(j(i)) = w(j(i)) + net.w(i);

end

net.sv   = sv;
net.w    = w;

% bye bye...

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [H] = discrete_mean_curvature(V,F)
  % DISCRETE_MEAN_CURVATURE Compute integrated mean curvature at each vertex.
  %
  % H = discrete_mean_curvature(V,F);
  %
  % Inputs:
  %   V  #V by 3 list of vertex positions
  %   F  #F by 3 list of triangle indices into V
  % Outputs:
  %   H  #V list of integrated mean-curvature values
  %
  % Examples:
  %   H = discrete_mean_curvature(V,F);
  %   M = massmatrix(V,F);
  %   tsurf(F,V,'CData',M\H,'EdgeColor','none',fphong);
  %
  %
  [A,C] = adjacency_dihedral_angle_matrix(V,F);
  [AI,AJ,AV] = find(A);
  [CI,CJ,CV] = find(C);
  assert(isequal(CI,AI));
  assert(isequal(CJ,AJ));
  l = edge_lengths(V,F);
  % index into F(:)
  opp = sub2ind(size(l),CI,CV);
  inc = [ ...
    sub2ind(size(l),CI,mod(CV+1-1,3)+1) ...
    sub2ind(size(l),CI,mod(CV+2-1,3)+1)];
  lV = l(opp);
  % From Keenan's slide: mean_i = ½ ∑{ij} lij φij
  % Extra 0.5 because each edge is counted twice
  % But why extra-extra 0.5?
  H = full(sparse(F(inc),1,0.5*0.5*0.5*repmat((pi-AV).*l(opp),1,2),size(V,1),1));
end

author: alecjacobson

repo: gptoolbox


text:
function [jc,ir,pr] = save_sparse(A)
% convert sparse matrix to jc, ir, pr
%
% Input
%  A - sparse matrix
% Output
%  jc -
%  ir -
%  pr -
%

[ir,j,pr] = find(A);
ir = ir - 1;

jc = zeros(1,size(A,2)+1);
for i = 1:numel(j), jc(j(i)+1) = jc(j(i)+1) + 1;end
jc = cumsum(jc);

%for i = 1:size(A,2)+1,	jc2(i) = sum(j < i);end
%plot([jc(100:200)',jc2(100:200)'])

author: mtex-toolbox

repo: mtex

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


text:
function varargout = rgb2xyz(varargin)
% VL_RGB2XYZ  Convert RGB color space to XYZ
%   J=VL_RGB2XYZ(I) converts the CIE RGB image I to the image J in
%   CIE XYZ format. CIE RGB has a white point of R=G=B=1.0
%
%   VL_RGB2XYZ(I,WS) uses the specified RGB working space WS. The
%   function supports the following RGB working spaces:
%
%   * `CIE'    E illuminant, gamma=2.2
%   * `Adobe'  D65 illuminant, gamma=2.2
%
%   The default workspace is CIE.
%
%   See also: VL_XYZ2RGB(), VL_HELP().
[varargout{1:nargout}] = vl_rgb2xyz(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
%% batman
% Below is a demonstration of the features of the |batman| function

%% Syntax
% |[x,y]=batman(n);|
% |[V]=batman(n);|

%% Description
% The |batman| function implements a particular version of the so called
% batmat-equation, a curve defining the batman logo. The input for this
% function is the number of desired points n. The user may request a sigle
% nx2 output array or two nx1 arrays (x and y coordinates).
%
% This is a MATLAB implementation of the parameterised Batman equation
% presented by Jerome White (http://www.talljerome.com/ @talljerome,
% https://youtu.be/oaIsCJw0QG8), in particular the form presented here:
% https://www.desmos.com/calculator/ajnzwedvql
% Modification: The batman is scaled to be 2 in width.

%% Examples

clear; close all; clc;

%% Example 1:

n=250; %Number of points on curve
[x,y]=batman(n);

%%

cFigure; hold on;
xlabel('x'); ylabel('y');
plot([x;x(1)],[y;y(1)],'k.-','LineWidth',3,'MarkerSize',25);
axis tight; axis equal;
set(gca,'FontSize',25);
grid on; box on;
gdrawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function PIVlab_capture_lensctrl_GUI
fh = findobj('tag', 'lens_control_window');
if isempty(fh)
	hgui=getappdata(0,'hgui');
	mainpos=get(hgui,'Position');
	[focus,aperture,lighting]=get_lens_status;
	lens_control_window = figure('numbertitle','off','MenuBar','none','DockControls','off','Name','Lens control','Toolbar','none','Units','characters','Position',[mainpos(1)+mainpos(3)-35 mainpos(2) 35 15+1.5+4],'tag','lens_control_window','visible','on','KeyPressFcn', @key_press,'resize','off','CloseRequestFcn',@CloseRequestFcn);
	set (lens_control_window,'Units','Characters');

	handles = guihandles; %alle handles mit tag laden und ansprechbar machen
	guidata(lens_control_window,handles)
	setappdata(0,'hlens',lens_control_window);

	parentitem = get(lens_control_window, 'Position');

	margin=1.5;

	panelheight=5;
	handles.aperturepanel = uipanel(lens_control_window, 'Units','characters', 'Position', [1 parentitem(4)-panelheight-1.5 parentitem(3)-2 panelheight],'title','Aperture control','fontweight','bold','tag','aperturepanel');
	handles.focuspanel = uipanel(lens_control_window, 'Units','characters', 'Position', [1 parentitem(4)-panelheight*2-1.5 parentitem(3)-2 panelheight],'title','Focus control','fontweight','bold','tag','focuspanel');
	handles.lightpanel = uipanel(lens_control_window, 'Units','characters', 'Position', [1 parentitem(4)-panelheight*3-1.5 parentitem(3)-2 panelheight],'title','Light control','fontweight','bold','tag','lightpanel');
	handles.anglepanel = uipanel(lens_control_window, 'Units','characters', 'Position', [1 parentitem(4)-panelheight*3-1.5-4 parentitem(3)-2 4],'title','Camera attitude','fontweight','bold','tag','anglepanel');

	%% Load last selected setting & setup Configurations
	warning off
	load ('PIVlab_capture_resources\PIVlab_capture_lensconfig.mat','lens_configurations','selected_lens_config_nr');
	warning on
	if ~exist('selected_lens_config_nr','var')
		selected_lens_config_nr = 2; %set default to zeiss Dimension
	end
	put('selected_lens_config',selected_lens_config_nr)
	save ('PIVlab_capture_resources\PIVlab_capture_lensconfig.mat','lens_configurations','selected_lens_config_nr');
	% New lens configurations can be added to the table by modifying the variable 'lens_configurations' in the file 'PIVlab_capture_lensconfig.mat :
	% Example: lens_configurations=addvars(lens_configurations,[500;2500;500;2500],'NewVariableNames','Generic lens')
	handles.configu = uicontrol(lens_control_window,'Style','popupmenu', 'String',lens_configurations.Properties.VariableNames,'Value',retr('selected_lens_config'),'Units','characters', 'Fontunits','points','Position',[1 parentitem(4)-1.5 parentitem(3)/3*2 1.5],'Tag','configu','TooltipString','Lens configuration. Sets the limits for the servo motors.','Callback',@configu_Callback);
	parentitem=get(lens_control_window, 'Position');
	item=[parentitem(3)/3*2 0 parentitem(3)/3 1.5];
	handles.lens_status = uicontrol(lens_control_window,'Style','edit','units','characters','HorizontalAlignment','center','position',[item(1)+margin parentitem(4)-1.5 item(3)-margin*1.5 item(4)],'String','N/A','tag', 'lens_status','FontName','FixedWidth','BackgroundColor',[1 0 0],'Foregroundcolor',[0 0 0],'Enable','inactive','Fontweight','bold');
	if selected_lens_config_nr ~=4
		all_features_visible='on';
	else
		all_features_visible='off';
	end
	put ('all_features_visible',all_features_visible);

	%% APERTURE
	parentitem=get(handles.aperturepanel, 'Position');
	item=[0 0 0 0];
	item=[parentitem(3)/2*0 item(2)+item(4) parentitem(3)/2*1 1.5];
	handles.aperture_open = uicontrol(handles.aperturepanel,'Style','pushbutton','String','Iris open','Units','characters', 'Fontunits','points','Position',[item(1)+margin parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'Callback', {@aperture_set,'open'} ,'TooltipString','Load image data');

	item=[parentitem(3)/2*1 item(2) parentitem(3)/2*1 1.5];
	handles.aperture_close = uicontrol(handles.aperturepanel,'Style','pushbutton','String','Iris close','Units','characters', 'Fontunits','points','Position',[item(1)+margin parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'Callback', {@aperture_set,'close'} ,'TooltipString','Load image data');

	item=[parentitem(3)/2*0 item(2)+item(4) parentitem(3)/2 1];
	handles.aperture_label = uicontrol(handles.aperturepanel,'Style','text','String','Iris [us]','units','characters','position',[item(1) parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'horizontalalignment','right');

	item=[parentitem(3)/2*1 item(2) parentitem(3)/3 1];
	handles.aperture_edit = uicontrol(handles.aperturepanel,'Style','edit','String',num2str(aperture),'units','characters','position',[item(1) parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'Callback',@aperture_edit_Callback,'tag','aperture_edit');

	%% FOCUS
	parentitem=get(handles.focuspanel, 'Position');
	item=[0 0 0 0];
	item=[0 item(2)+item(4) parentitem(3)/3 1.5];
	handles.focus_close = uicontrol(handles.focuspanel,'Style','pushbutton','String','Near','Units','characters', 'Fontunits','points','Position',[item(1)+margin parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'Callback', {@focus_set,'near'} ,'TooltipString','Load image data');

	item=[parentitem(3)/3*1 item(2) parentitem(3)/3 1.5];
	handles.focus_auto = uicontrol(handles.focuspanel,'Style','pushbutton','String','Auto','Units','characters', 'Fontunits','points','Position',[item(1)+margin parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'Callback', {@focus_set,'auto'} ,'TooltipString','Load image data');

	item=[parentitem(3)/3*2 item(2) parentitem(3)/3 1.5];
	handles.focus_far = uicontrol(handles.focuspanel,'Style','pushbutton','String','Far','Units','characters', 'Fontunits','points','Position',[item(1)+margin parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'Callback', {@focus_set,'far'} ,'TooltipString','Load image data');

	item=[parentitem(3)/2*0 item(2)+item(4) parentitem(3)/2 1];
	handles.focus_label = uicontrol(handles.focuspanel,'Style','text','String','Focus [us]','units','characters','position',[item(1) parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'horizontalalignment','right');

	item=[parentitem(3)/2*1 item(2) parentitem(3)/3 1];
	handles.focus_edit = uicontrol(handles.focuspanel,'Style','edit','String',num2str(focus),'units','characters','position',[item(1) parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'Callback',@focus_edit_Callback,'tag','focus_edit');
	setappdata(lens_control_window,'handle_to_focus_edit_field',handles.focus_edit); %needed so other files can edit the contents.

	%% LIGHT
	parentitem=get(handles.lightpanel, 'Position');
	item=[0 0 0 0];
	item=[parentitem(3)/2*0 item(2)+item(4) parentitem(3)/2*1 1.5];
	handles.light_on = uicontrol(handles.lightpanel,'Style','pushbutton','String','Light on','Units','characters', 'Fontunits','points','Position',[item(1)+margin parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'Callback', {@light_switch,'on'} ,'TooltipString','Load image data');

	item=[parentitem(3)/2*1 item(2) parentitem(3)/2*1 1.5];
	handles.light_off = uicontrol(handles.lightpanel,'Style','pushbutton','String','Light off','Units','characters', 'Fontunits','points','Position',[item(1)+margin parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'Callback', {@light_switch,'off'},'TooltipString','Load image data');

	item=[parentitem(3)/2*0 item(2)+item(4) parentitem(3)/2 1];
	handles.light_label = uicontrol(handles.lightpanel,'Style','text','String','Light status','units','characters','position',[item(1) parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'horizontalalignment','right');
	item=[parentitem(3)/2*1 item(2) parentitem(3)/3 1];
	handles.light_edit = uicontrol(handles.lightpanel,'Style','edit','String',num2str(lighting),'units','characters','position',[item(1) parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'Callback',@focus_edit_Callback,'tag','light_edit','enable','off');

	configu_Callback(handles.configu,[]) %execute callback and set servo limits

	%% Angle feedback
	parentitem=get(handles.anglepanel, 'Position');
	item=[0 0 parentitem(3) 1];
	handles.angle_measure = uicontrol(handles.anglepanel,'Style','checkbox','String','Measure angle','units','characters','position',[item(1) parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'horizontalalignment','left','Callback',@angle_measure_Callback,'tag','angle_measure');

	item=[parentitem(3)/2*0 item(2)+item(4) parentitem(3)/3 1.5];
	handles.pitch = uicontrol(handles.anglepanel,'Style','text','String','Pitch: 0','units','characters','position',[item(1) parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'horizontalalignment','left','tag','pitch');

	item=[parentitem(3)/2*1 item(2) parentitem(3)/3 1.5];
	handles.roll = uicontrol(handles.anglepanel,'Style','text','String','Roll: 0','units','characters','position',[item(1) parentitem(4)-item(4)-margin-item(2) item(3)-margin*2 item(4)],'horizontalalignment','left','tag','roll');
	limit_displayed_features
	find_devices

else %Figure handle does already exist --> bring UI to foreground.
	figure(fh)
end

function limit_displayed_features
all_features_visible=retr('all_features_visible');
aperturepanel=findobj('Tag','aperturepanel');
lightpanel=findobj('Tag','lightpanel');
anglepanel=findobj('Tag','anglepanel');

set(aperturepanel,'Visible',all_features_visible)
set(lightpanel,'Visible',all_features_visible)
set(anglepanel,'Visible',all_features_visible)

function find_devices
hgui = getappdata(0,'hgui');
serpo=getappdata(hgui,'serpo');
try
	serpo.Port; %is there no other way to determine if serialport is working...?
	%configureTerminator(serpo,'CR/LF');
	alreadyconnected=1;
catch
	alreadyconnected=0;
	delete(serpo)
end
if alreadyconnected==1
	flush(serpo)
	writeline(serpo,'WhatIsTheAngle?');
	warning off
	serial_answer=readline(serpo);
	warning on
	handles=gethand;
	lens_available=strfind(serial_answer,'Measured_Roll:');
	if ~isempty(lens_available) &&  lens_available~=0
		set(handles.lens_status, 'Backgroundcolor',[0 1 0])
		set(handles.lens_status, 'String','OK')
	end
end

function angle_measure_Callback (inpt,~)
if inpt.Value == 1 %on
	t = timer;
	t.Tag='lens_timer';
	t.StartDelay = 0.1;
	t.ExecutionMode = 'fixedRate';
	t.Period = 0.5;
	t.TimerFcn = @lens_timer_tick_fcn;
	start(t)
	setappdata(0,'handle_to_lens_timer_checkbox',inpt);
else
	stop(timerfind)
	delete(timerfind)
end

function lens_timer_tick_fcn(~,~)
hgui = getappdata(0,'hgui');
serpo=getappdata(hgui,'serpo');
try
	serpo.Port; %is there no other way to determine if serialport is working...?
	alreadyconnected=1;
catch
	alreadyconnected=0;
	delete(serpo)
end
if alreadyconnected==1
	flush(serpo)
	writeline(serpo,'WhatIsTheAngle?');
	warning off
	serial_answer=readline(serpo);
	warning on
	handles=gethand;
	if ~isempty (serial_answer) & strfind(serial_answer,'Measured_Roll:')==1
		Roll=serial_answer{1}(strfind(serial_answer,'Measured_Roll:')+14:strfind(serial_answer,char(9))-1);
		Pitch=serial_answer{1}(strfind(serial_answer,'Measured_Pitch:')+15:end);
		Roll=str2double(Roll)/100+retr('Roll_Offset');
		Pitch=str2double(Pitch)/100+retr('Pitch_Offset');
		set(handles.pitch,'String',['P: ' num2str(Pitch)])
		set(handles.roll,'String',['R: ' num2str(Roll)])
	else
		set(handles.pitch,'String','No reply')
		set(handles.roll,'String','No reply')
	end
end

function configu_Callback (inpt,~)
load ('PIVlab_capture_lensconfig.mat','lens_configurations');
focus_servo_lower_limit = lens_configurations{1,inpt.Value};
focus_servo_upper_limit = lens_configurations{2,inpt.Value};
aperture_servo_lower_limit = lens_configurations{3,inpt.Value};
aperture_servo_upper_limit = lens_configurations{4,inpt.Value};
Pitch_Offset = lens_configurations{5,inpt.Value};
Roll_Offset = lens_configurations{6,inpt.Value};

put('focus_servo_lower_limit',focus_servo_lower_limit)
put('focus_servo_upper_limit',focus_servo_upper_limit)
put('aperture_servo_lower_limit',aperture_servo_lower_limit)
put('aperture_servo_upper_limit',aperture_servo_upper_limit)
put('Pitch_Offset',Pitch_Offset)
put('Roll_Offset',Roll_Offset)

put('selected_lens_config',inpt.Value)
selected_lens_config_nr=inpt.Value;
save ('PIVlab_capture_resources\PIVlab_capture_lensconfig.mat','lens_configurations','selected_lens_config_nr');

handles=gethand;
focus=retr('focus');
aperture=retr('aperture');
lighting=retr('lighting');
if isempty(focus)
	focus=1500;
end
if isempty(aperture)
	aperture=1500;
end
if isempty(lighting)
	lighting=0;
end
set (handles.aperture_edit,'String',num2str(aperture))
set (handles.focus_edit,'String',num2str(focus))
set (handles.light_edit,'String',num2str(lighting))
%focus_edit_Callback(handles.aperture_edit,[])
%pause(0.2)
%aperture_edit_Callback(handles.focus_edit,[])
if selected_lens_config_nr ~=4
	all_features_visible='on';
else
	all_features_visible='off';
end
put ('all_features_visible',all_features_visible);
limit_displayed_features

function focus_set (~,~,inpt)

focus_step=100;
[focus,aperture,lighting]=get_lens_status;
if strmatch(inpt,'auto')
	put('sharpness_enabled',0);
	put('hist_enabled',0);
	if retr('capturing')==1 %camera is recording
		if retr('autofocus_enabled') == 1 %user pressed button while autofocus is running: Stop autofocus.
			put('autofocus_enabled',0); %toggles the autofocus_enabled variable. That is checked in PIVlab_capture_pco after each frame capture
		else
			put('autofocus_enabled',1);
		end
		%move to lower limit
		PIVlab_capture_lensctrl (retr('focus_servo_lower_limit'), retr('aperture'),retr('lighting'))
	else
		put('autofocus_enabled',0);
	end
end

if strmatch(inpt,'far')
	if focus>=retr('focus_servo_lower_limit')+focus_step
		focus=focus-focus_step;
	else
		focus=retr('focus_servo_lower_limit');
	end
	PIVlab_capture_lensctrl (focus, aperture,lighting)
end
if strmatch(inpt,'near')
	if focus<=retr('focus_servo_upper_limit')-focus_step
		focus=focus+focus_step;
	else
		focus=retr('focus_servo_upper_limit');
	end
	PIVlab_capture_lensctrl (focus, aperture,lighting)
end
%put('focus',focus);

update_edit_fields

function [focus,aperture,lighting]=get_lens_status
try %try to switch of camera angle report
	stop(timerfind)
	delete(timerfind)
	set(getappdata(0,'handle_to_lens_timer_checkbox'),'Value',0)
catch
end
focus=retr('focus');
aperture=retr('aperture');
lighting=retr('lighting');
if isempty(focus)
	focus=1500;
end
if isempty(aperture)
	aperture=1500;
end
if isempty(lighting)
	lighting=0;
end

function aperture_set (~,~,inpt)
[focus,aperture,lighting]=get_lens_status;
aperture_step=100;
if strmatch(inpt,'open')
	if aperture>=retr('aperture_servo_lower_limit')+aperture_step
		aperture=aperture-aperture_step;
	else
		aperture=retr('aperture_servo_lower_limit');
	end
end
if strmatch(inpt,'close')
	if aperture<=retr('aperture_servo_upper_limit')-aperture_step
		aperture=aperture+aperture_step;
	else
		aperture=retr('aperture_servo_upper_limit');
	end
end
%put('aperture',aperture);
PIVlab_capture_lensctrl (focus, aperture,lighting)
update_edit_fields

function aperture_edit_Callback(caller,~)
[focus,~,lighting]=get_lens_status;
aperture=str2double(caller.String);
if aperture > retr('aperture_servo_upper_limit')
	aperture =retr('aperture_servo_upper_limit');
end
if aperture < retr('aperture_servo_lower_limit')
	aperture =retr('aperture_servo_lower_limit');
end
caller.String=num2str(aperture);
%put('aperture',aperture);
PIVlab_capture_lensctrl (focus, aperture,lighting)

function focus_edit_Callback(caller,~)
[~,aperture,lighting]=get_lens_status;
focus=str2double(caller.String);
if focus > retr('focus_servo_upper_limit')
	focus =retr('focus_servo_upper_limit');
end
if focus< retr('focus_servo_lower_limit')
	focus =retr('focus_servo_lower_limit');
end
caller.String=num2str(focus);
%put('focus',focus);
PIVlab_capture_lensctrl (focus, aperture,lighting)

function update_edit_fields(~,~)
[focus,aperture,lighting]=get_lens_status;
handles=gethand;
set(handles.focus_edit ,'String', num2str(focus))
set(handles.aperture_edit ,'String', num2str(aperture))
set(handles.light_edit ,'String', num2str(lighting))

function light_switch (~, ~,inpt)
[focus,aperture,lighting]=get_lens_status;
if strmatch(inpt,'on')
	lighting=1;
end
if strmatch(inpt,'off')
	lighting=0;
end
%put('lighting',lighting);
PIVlab_capture_lensctrl (focus, aperture,lighting)
update_edit_fields

function put(name, what)
hgui=getappdata(0,'hgui');
setappdata(hgui, name, what);

function var = retr(name)
hgui=getappdata(0,'hgui');
var=getappdata(hgui, name);

function handles=gethand
hlens=getappdata(0,'hlens');
handles=guihandles(hlens);

function CloseRequestFcn(hObject, ~, ~)
try
	stop(timerfind)
	delete(timerfind)
catch
end
try
	delete(hObject);
catch
	delete(gcf);
end
put('autofocus_enabled',0);

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
classdef ODE5Integrator < AbstractFixedStepIntegrator
    %ODE5Integrator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        options(1,1) FixedStepSizeIntegratorOptions = FixedStepSizeIntegratorOptions()
        integratorEnum = IntegratorEnum.ODE5;
    end

    methods
        function obj = ODE5Integrator(options)
            if(nargin > 0)
                obj.options = options;
            else
                obj.options = FixedStepSizeIntegratorOptions();
            end
        end

        function [t,y,te,ye,ie] = integrate(obj, odefun, tspan, y0, evtsFunc, odeOutputFun)
%             odeSetOptions = obj.options.getIntegratorOptions();
            odeSetOptions = odeset();
            optionsToUse = odeset(odeSetOptions, 'Events',evtsFunc, 'OutputFcn',odeOutputFun);

            [t,y, te,ye,ie] = AbstractFixedStepIntegrator.integrate(odefun, tspan, y0, optionsToUse);
        end

        function options = getOptions(obj)
            options = obj.options;
        end
    end

    methods(Static,Access=protected)
        function [A,B,C] = getButcherTableauData()
            A = [  1/5,          0,           0,            0,         0
                   3/40,         9/40,        0,            0,         0
                   44/45        -56/15,       32/9,         0,         0
                   19372/6561,  -25360/2187,  64448/6561,  -212/729,   0
                   9017/3168,   -355/33,      46732/5247,   49/176,   -5103/18656];

            B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];

            C = [1/5; 3/10; 4/5; 8/9; 1];

            A = A.';
            B = B(:);
            C = C(:);
        end
    end
end

author: Arrowstar

repo: ksptot


text:
function [kf] = ext_kalman_filter(F, G, H, u, Q, R, time, last_measure, kf, P)

%EXT_KALMAN_FILTER - This is the extended Kalman filter.
%
% Synthax:
%
% Inputs:
%   uu          - state variables, delta, wind
%   kf.x_prev   - previous a posteriori state estimate
%   kf.P_prev   - previous a posteriori estimate covariance
%   last_measure- (time, sensor measure)
%   P           - parameters
%
% Outputs:
%   kf.x        - a priori/a posteriori state estimate
%   kf.P        - a priori/a posteriori estimate covariance
%   kf.z_res    - post-fit residual

% A -> F
% B -> G
% C -> H

% Q - covariance of the zero-mean Gaussian noise on x_dot
% R - covariance of the zero-mean Gaussian noise on y

time_m = last_measure(:,1);
z_m    = last_measure(:,2);

% prediction (a priori)
kf.x = F * kf.x_prev    + G * u;    % a priori state estimate
kf.P = F * kf.P_prev * (F') + Q;    % a priori estimate covariance

% update (a posteriori)
if time - time_m < P.dt
    kf.z_tilde = z_m - H*kf.x;      % innovation (a priori residual)
    S = R + H*kf.P*(H');            % innovation covariance
    K = kf.P*(H')/S;                 % optimal kalman gain
    kf.x = kf.x + K*kf.z_tilde;     % a posteriori state estimate
    kf.P = kf.P - K*S*(K');         % a posteriori estimate covariance
end

kf.z_res = z_m - H*kf.x;            % a posteriori residual

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab