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

    properties
        stateLog LaunchVehicleStateLog
        orbitModel(1,1) AbstractElementSet = KeplerianElementSet.getDefaultElements();

        %Inherit Time
        inheritTime(1,1) logical = false;
        inheritTimeFrom(1,1) InheritStateEnum = InheritStateEnum.InheritFromLastState;
        inheritTimeFromEvent LaunchVehicleEvent = LaunchVehicleEvent.empty(1,0);

        %Inherit State
        inheritPosVel(1,1) logical = false;
        inheritPosVelFrom(1,1) InheritStateEnum = InheritStateEnum.InheritFromLastState;
        inheritPosVelFromEvent LaunchVehicleEvent = LaunchVehicleEvent.empty(1,0);

        %Inherit Stage States
        stageStates(1,:) SetKinematicStateStageState
        engineStates(1,:) SetKinematicStateEngineState
        tankStates(1,:) SetKinematicStateTankState
        epsSinkStates(1,:) SetKinematicStateEPSinkState
        epsSrcStates(1,:) SetKinematicStateEPSrcState
        epsStorageStates(1,:) SetKinematicStateEPStorageState

        %Inherit all other state elements
        inheritStateElems(1,1) logical = true;
        inheritStateElemsFrom(1,1) InheritStateEnum = InheritStateEnum.InheritFromLastState;
        inheritStateElemsFromEvent LaunchVehicleEvent = LaunchVehicleEvent.empty(1,0);

        optVar SetKinematicStateActionVariable
    end

    %deprecated
    properties
        inheritStageStates(1,1) logical = true;
        inheritStageStatesFrom(1,1) InheritStateEnum = InheritStateEnum.InheritFromLastState;
        inheritStageStatesFromEvent LaunchVehicleEvent = LaunchVehicleEvent.empty(1,0);
    end

    properties(Dependent)
        time(1,1) double
        centralBody(1,1) KSPTOT_BodyInfo
    end

    methods
        function obj = SetKinematicStateAction(stateLog, orbitModel)
            if(nargin > 0)
                obj.stateLog = stateLog;
                obj.orbitModel = orbitModel;
            end

            obj.id = rand();
        end

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

        function set.time(obj, newTime)
            obj.orbitModel.time = newTime;
        end

        function time = get.centralBody(obj)
            time = obj.orbitModel.frame.getOriginBody();
        end

        function set.centralBody(obj, newCentralBody)
            obj.orbitModel.frame.setOriginBody(newCentralBody);
        end

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

            % Time
            if(obj.inheritTime)
                if(obj.inheritTimeFrom == InheritStateEnum.InheritFromSpecifiedEvent && ...
                   not(isempty(obj.inheritTimeFromEvent)) && ...
                   not(isempty(obj.inheritTimeFromEvent.getEventNum())) && ...
                   obj.inheritTimeFromEvent.getEventNum() > 0)

                    timeEvtStateLog = obj.stateLog.getLastStateLogForEvent(obj.inheritTimeFromEvent);
                    if(not(isempty(timeEvtStateLog)))
                        newStateLogEntry.time = timeEvtStateLog(end).time;
                        obj.orbitModel.time = newStateLogEntry.time;
                    end
                end
            else
                newStateLogEntry.time = obj.orbitModel.time;
            end

            % Orbit
            if(obj.inheritPosVel)
                if(obj.inheritPosVelFrom == InheritStateEnum.InheritFromSpecifiedEvent && ...
                   not(isempty(obj.inheritPosVelFromEvent)) && ...
                   not(isempty(obj.inheritPosVelFromEvent.getEventNum())) && ...
                   obj.inheritPosVelFromEvent.getEventNum() > 0)

                    posVelEvtStateLog = obj.stateLog.getLastStateLogForEvent(obj.inheritPosVelFromEvent);
                    if(not(isempty(posVelEvtStateLog)))
                        newStateLogEntry.position = posVelEvtStateLog(end).position;
                        newStateLogEntry.velocity = posVelEvtStateLog(end).velocity;
                        newStateLogEntry.centralBody = posVelEvtStateLog(end).centralBody;
                    end
                end
            else
                cartElemSet = obj.orbitModel.convertToCartesianElementSet();

                newStateLogEntry.setCartesianElementSet(cartElemSet);
            end

            for(i=1:length(obj.stageStates)) %#ok<*UNRCH>
                obj.stageStates(i).updateStateLogEntry(newStateLogEntry, obj.stateLog)
            end

            for(i=1:length(obj.engineStates))
                obj.engineStates(i).updateStateLogEntry(newStateLogEntry, obj.stateLog)
            end

            for(i=1:length(obj.tankStates))
                obj.tankStates(i).updateStateLogEntry(newStateLogEntry, obj.stateLog)
            end

            for(i=1:length(obj.epsSinkStates))
                obj.epsSinkStates(i).updateStateLogEntry(newStateLogEntry, obj.stateLog)
            end

            for(i=1:length(obj.epsSrcStates))
                obj.epsSrcStates(i).updateStateLogEntry(newStateLogEntry, obj.stateLog)
            end

            for(i=1:length(obj.epsStorageStates))
                obj.epsStorageStates(i).updateStateLogEntry(newStateLogEntry, obj.stateLog)
            end

            %Comment out this block once done
%             if(obj.inheritStageStates)
%                 if(obj.inheritStageStatesFrom == InheritStateEnum.InheritFromSpecifiedEvent && ...
%                    not(isempty(obj.inheritStageStatesFromEvent)) && ...
%                    not(isempty(obj.inheritStageStatesFromEvent.getEventNum())) && ...
%                    obj.inheritStageStatesFromEvent.getEventNum() > 0)
%
%
%                     stageStateEvtState = obj.stateLog.getLastStateLogForEvent(obj.inheritStageStatesFromEvent);
%                     stageStateEvtState = stageStateEvtState(end);
%                     if(not(isempty(stageStateEvtState)))
%                         lvState = stageStateEvtState.lvState.deepCopy();
%                         newStateLogEntry.lvState = lvState;
%
%                         for(i=1:length(stageStateEvtState.stageStates)) %#ok<NO4LP>
%                             newStateLogEntry.stageStates(i) = stageStateEvtState.stageStates(i).deepCopy(true, lvState);
%                         end
%                     end
%                 end
%
%             else
%                 %do nothing, there is no behavior for this yet
%                 warning('SetKinematicStateAction.inheritStageStates must be set to true, no behavior yet defined when this is set to false.');
%             end

            %Misc states
            if(obj.inheritStateElems)
              if(obj.inheritStateElemsFrom == InheritStateEnum.InheritFromSpecifiedEvent && ...
                   not(isempty(obj.inheritStateElemsFromEvent)) && ...
                   not(isempty(obj.inheritStateElemsFromEvent.getEventNum())) && ...
                   obj.inheritStateElemsFromEvent.getEventNum() > 0)

                    stateElemsEvtState = obj.stateLog.getLastStateLogForEvent(obj.inheritStateElemsFromEvent);
                    stateElemsEvtState = stateElemsEvtState(end);

                    if(not(isempty(stateElemsEvtState)))
                        newStateLogEntry.aero = stateElemsEvtState.aero.deepCopy();
                        newStateLogEntry.thirdBodyGravity = stateElemsEvtState.thirdBodyGravity.copy();
                        newStateLogEntry.stopwatchStates = stateElemsEvtState.stopwatchStates.copy();
                        newStateLogEntry.extremaStates = stateElemsEvtState.extremaStates.copy();
                        newStateLogEntry.steeringModel = stateElemsEvtState.steeringModel; %NEED TO WARN USE TO REINITIALIZE THROTTLE AND STEERING MODELS!
                        newStateLogEntry.throttleModel = stateElemsEvtState.throttleModel;
                    end
              end
            else
                %do nothing, there is no behavior for this yet
                warning('SetKinematicStateAction.inheritStateElems must be set to true, no behavior yet defined when this is set to false.');
            end
        end

        function initAction(obj, initialStateLogEntry)
            %nothing
        end

        function name = getName(obj)
            name = 'Set Kinematic State';
        end

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

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

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

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

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

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

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

        function tf = usesEvent(obj, event)
            tf = false;

            tankEvtsTf = logical([]);
            for(i=1:length(obj.tankStates))
                tankState = obj.tankStates(i);
                if(tankState.inheritTankStateFrom == InheritStateEnum.InheritFromSpecifiedEvent && ...
                   not(isempty(tankState.inheritTankStateFromEvent)) && tankState.inheritTankStateFromEvent == event)
                    tankEvtsTf(end+1) = tankState.inheritTankStateFromEvent == event; %#ok<AGROW>
                end
            end

            epsSrcEvtsTf = logical([]);
            for(i=1:length(obj.epsStorageStates))
                storageState = obj.epsStorageStates(i);

                if(storageState.inheritStorageStateFrom == InheritStateEnum.InheritFromSpecifiedEvent && ...
                   not(isempty(storageState.inheritStorageStateFromEvent)) && storageState.inheritStorageStateFromEvent == event)
                    epsSrcEvtsTf(end+1) = storageState.inheritStorageStateFromEvent == event; %#ok<AGROW>
                end
            end

            if((obj.inheritTime && obj.inheritTimeFrom == InheritStateEnum.InheritFromSpecifiedEvent && not(isempty(obj.inheritTimeFromEvent)) && obj.inheritTimeFromEvent == event) || ...
               (obj.inheritPosVel && obj.inheritPosVelFrom == InheritStateEnum.InheritFromSpecifiedEvent && not(isempty(obj.inheritPosVelFromEvent)) && obj.inheritPosVelFromEvent == event) || ...
               (obj.inheritStateElems && obj.inheritStateElemsFrom == InheritStateEnum.InheritFromSpecifiedEvent && not(isempty(obj.inheritStateElemsFromEvent)) && obj.inheritStateElemsFromEvent == event) || ...
               any(tankEvtsTf) || any(epsSrcEvtsTf))
                tf = true;
            end
        end

        function [tf, vars] = hasActiveOptimVar(obj)
            tf = false;
            vars = AbstractOptimizationVariable.empty(0,1);

            if(not(isempty(obj.optVar)))
                tf = any(obj.optVar.getUseTfForVariable());
                vars(end+1) = obj.optVar;
            end
        end

        function generateLvComponentElements(obj)
            lvdData = obj.event.lvdData;

            if(isempty(obj.stageStates))
                stages = lvdData.launchVehicle.stages;
                for(j=1:length(stages))
                    stageState = SetKinematicStateStageState(stages(j));
                    stageState.inheritStageState = obj.inheritStageStates;
                    stageState.inheritStageStateFrom = obj.inheritStageStatesFrom;
                    stageState.inheritStageStateFromEvent = obj.inheritStageStatesFromEvent;

                    obj.stageStates(j) = stageState;
                end
            end
            addlistener(lvdData.launchVehicle, 'StageAdded',   @(src, evt) stageAddedClbk(obj, src, evt));
            addlistener(lvdData.launchVehicle, 'StageDeleted', @(src, evt) stageDeletedClbk(obj, src, evt));

            if(isempty(obj.engineStates))
                [~, engines] = lvdData.launchVehicle.getEnginesListBoxStr();
                for(j=1:length(engines))
                    engineState = SetKinematicStateEngineState(engines(j));
                    engineState.inheritEngineState = obj.inheritStageStates;
                    engineState.inheritEngineStateFrom = obj.inheritStageStatesFrom;
                    engineState.inheritEngineStateFromEvent = obj.inheritStageStatesFromEvent;

                    obj.engineStates(j) = engineState;
                end
            end
            addlistener(lvdData.launchVehicle, 'EngineAdded',   @(src, evt) engineAddedClbk(obj, src, evt));
            addlistener(lvdData.launchVehicle, 'EngineDeleted', @(src, evt) engineDeletedClbk(obj, src, evt));

            if(isempty(obj.tankStates))
                [~, storages] = lvdData.launchVehicle.getTanksListBoxStr();
                for(j=1:length(storages))
                    tankState = SetKinematicStateTankState(storages(j));
                    tankState.inheritTankState = obj.inheritStageStates;
                    tankState.inheritTankStateFrom = obj.inheritStageStatesFrom;
                    tankState.inheritTankStateFromEvent = obj.inheritStageStatesFromEvent;

                    obj.tankStates(j) = tankState;
                end
            end
            addlistener(lvdData.launchVehicle, 'TankAdded',   @(src, evt) tankAddedClbk(obj, src, evt));
            addlistener(lvdData.launchVehicle, 'TankDeleted', @(src, evt) tankDeletedClbk(obj, src, evt));

            if(isempty(obj.epsSinkStates))
                [~, sinks] = lvdData.launchVehicle.getPowerSinksListBoxStr();
                for(j=1:length(sinks))
                    sinkState = SetKinematicStateEPSinkState(sinks(j));
                    sinkState.inheritSinkState = obj.inheritStageStates;
                    sinkState.inheritSinkStateFrom = obj.inheritStageStatesFrom;
                    sinkState.inheritSinkStateFromEvent = obj.inheritStageStatesFromEvent;

                    obj.epsSinkStates(j) = sinkState;
                end
            end
            addlistener(lvdData.launchVehicle, 'EpsSinkAdded',   @(src, evt) epsSinkAddedClbk(obj, src, evt));
            addlistener(lvdData.launchVehicle, 'EpsSinkDeleted', @(src, evt) epsSinkDeletedClbk(obj, src, evt));

            if(isempty(obj.epsSrcStates))
                [~, srcs] = lvdData.launchVehicle.getPowerSrcsListBoxStr();
                for(j=1:length(srcs))
                    srcState = SetKinematicStateEPSrcState(srcs(j));
                    srcState.inheritSrcState = obj.inheritStageStates;
                    srcState.inheritSrcStateFrom = obj.inheritStageStatesFrom;
                    srcState.inheritSrcStateFromEvent = obj.inheritStageStatesFromEvent;

                    obj.epsSrcStates(j) = srcState;
                end
            end
            addlistener(lvdData.launchVehicle, 'EpsSrcAdded',   @(src, evt) epsSrcAddedClbk(obj, src, evt));
            addlistener(lvdData.launchVehicle, 'EpsSrcDeleted', @(src, evt) epsSrcDeletedClbk(obj, src, evt));

            if(isempty(obj.epsStorageStates))
                [~, storages] = lvdData.launchVehicle.getPowerStoragesListBoxStr();
                for(j=1:length(storages))
                    storageState = SetKinematicStateEPStorageState(storages(j));
                    storageState.inheritStorageState = obj.inheritStageStates;
                    storageState.inheritStorageStateFrom = obj.inheritStageStatesFrom;
                    storageState.inheritStorageStateFromEvent = obj.inheritStageStatesFromEvent;

                    obj.epsStorageStates(j) = storageState;
                end
            end
            addlistener(lvdData.launchVehicle, 'EpsStorageAdded',   @(src, evt) epsStorageAddedClbk(obj, src, evt));
            addlistener(lvdData.launchVehicle, 'EpsStorageDeleted', @(src, evt) epsStorageDeletedClbk(obj, src, evt));

        end
    end

    methods(Access=private)
        function stageAddedClbk(obj, src, evt)
            stage = evt.stage;

            stageState = SetKinematicStateStageState(stage);
            obj.stageStates(end+1) = stageState;
        end

        function stageDeletedClbk(obj, src, evt)
            stage = evt.stage;

            obj.stageStates([obj.stageStates.stage] == stage) = [];
        end

        function engineAddedClbk(obj, src, evt)
            engine = evt.engine;

            engineState = SetKinematicStateEngineState(engine);
            obj.engineStates(end+1) = engineState;
        end

        function engineDeletedClbk(obj, src, evt)
            engine = evt.engine;

            obj.engineStates([obj.engineStates.engine] == engine) = [];
        end

        function tankAddedClbk(obj, src, evt)
            tank = evt.tank;

            tankState = SetKinematicStateTankState(tank);
            obj.tankStates(end+1) = tankState;

            SetKinematicStateActionVariable(tankState); %doesn't need to get stored somewhere, it'll get stored on the tankState
        end

        function tankDeletedClbk(obj, src, evt)
            tank = evt.tank;

            obj.tankStates([obj.tankStates.tank] == tank) = [];
        end

        function epsSinkAddedClbk(obj, src, evt)
            sink = evt.sink;

            sinkState = SetKinematicStateEPSinkState(sink);
            obj.epsSinkStates(end+1) = sinkState;
        end

        function epsSinkDeletedClbk(obj, src, evt)
            sink = evt.sink;

            obj.epsSinkStates([obj.epsSinkStates.sink] == sink) = [];
        end

        function epsSrcAddedClbk(obj, src, evt)
            src = evt.src;

            srcState = SetKinematicStateEPSrcState(src);
            obj.epsSrcStates(end+1) = srcState;
        end

        function epsSrcDeletedClbk(obj, ~, evt)
            src = evt.src;

            obj.epsSrcStates([obj.epsSrcStates.src] == src) = [];
        end

        function epsStorageAddedClbk(obj, src, evt)
            storage = evt.storage;

            storageState = SetKinematicStateEPStorageState(storage);
            obj.epsStorageStates(end+1) = storageState;
        end

        function epsStorageDeletedClbk(obj, src, evt)
            storage = evt.storage;

            obj.epsStorageStates([obj.epsStorageStates.storage] == storage) = [];
        end
    end

    methods(Static)
        function addActionTf = openEditActionUI(action, lv)
            action.generateLvComponentElements();

            output = AppDesignerGUIOutput({false});
            lvd_EditActionSetKinematicStateGUI_App(action, lv.lvdData, output);

            addActionTf = output.output{1};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef ParameterServer < handle
  properties
    params
    paramsRegister
    method
    prefix
    memoryMapFile
    memoryMap
    pinnedMemory
    otherLabs
    inplace
    tflowOpts
  end

  methods
    function obj = ParameterServer(varargin)
      obj.prefix = 'mcn' ;
      obj.params = struct('name', {}, 'shape', {}, ...
                          'dataType', {}, 'deviceType', {}, ...
                          'value', {}) ;
      obj.paramsRegister = struct() ;
      switch computer('arch')
        case 'win64'
          obj.method = 'mmap' ;
        otherwise
          obj.method = 'tflow' ;
      end
      obj.pinnedMemory = true ;
      obj.inplace = true ;
      obj.tflowOpts = { } ;
      obj = vl_argparse(obj, varargin) ;

      % Compute the memory map file from the prefix.  This is used only by
      % the 'mmap' mode.
      obj.memoryMapFile = fullfile(tempdir, sprintf('%s.bin',obj.prefix)) ;
    end

    function register(obj, name, shape, dataType, deviceType)
      x = zeros(shape, dataType) ;
      if strcmp(deviceType,'gpu')
        x = gpuArray(x) ;
      end
      obj.params(end+1).name = name ;
      obj.params(end).shape = shape ;
      obj.params(end).dataType = dataType ;
      obj.params(end).deviceType = deviceType ;
      obj.params(end).isGPU = strcmp(deviceType, 'gpu') ;
      obj.params(end).value =  x ;
      obj.paramsRegister.(name) = numel(obj.params) ;
    end

    function delete(obj)
      obj.stop()
    end

    function start(obj)
      switch obj.method
        case 'mmap'
          if labindex == 1 && exist(obj.memoryMapFile)
            delete(obj.memoryMapFile) ;
          end
          labBarrier() ;
          obj = startWithMMap(obj) ;
        case 'tflow'
          obj = startWithTFlow(obj) ;
      end
    end

    function stop(obj)
      switch obj.method
        case 'mmap', stopWithMMap(obj) ;
        case 'tflow', stopWithTFlow(obj) ;
      end
    end

    function push(obj, name, value)
      p = obj.paramsRegister.(name) ;
      obj.pushWithIndex(p, value) ;
    end

    function value = pull(obj, name)
      p = obj.paramsRegister.(name) ;
      value = obj.pullWithIndex(p) ;
    end

    function pushWithIndex(obj, p, value)
      switch obj.method
        case 'mmap'
          obj.params(p).value = value ;

        case 'tflow'
          if obj.params(p).isGPU && obj.inplace
            obj.params(p).value = value ;
            vl_tflow('push', obj.params(p).name, value, 'inplace', obj.tflowOpts{:}) ;
          else
            vl_tflow('push', obj.params(p).name, value, obj.tflowOpts{:}) ;
          end

        otherwise
          assert(false) ;
      end
    end

    function value = pullWithIndex(obj,p)
      switch obj.method
        case 'mmap'
          for l = obj.otherLabs
            obj.params(p).value = obj.params(p).value ...
                + obj.memoryMap.Data(l).(obj.params(p).name) ;
          end
          value = obj.params(p).value ;

        case 'tflow'
          if obj.params(p).isGPU && obj.inplace
            vl_tflow('pull', obj.params(p).name, 'inplace', obj.tflowOpts{:}) ;
            value = obj.params(p).value ;
          else
            value = vl_tflow('pull', obj.params(p).name, obj.tflowOpts{:}) ;
          end
      end
    end

    function sync(obj)
      switch obj.method
        case 'mmap'
          for p = 1:numel(obj.params)
            if ~obj.pinnedMemory || ~obj.params(p).isGPU
              obj.memoryMap.Data(labindex).(obj.params(p).name) = ...
                  gather(obj.params(p).value) ;
            else
              vl_cudatool('cudaCopyDeviceToHost', ...
                       obj.memoryMap.Data(labindex).(obj.params(p).name), ...
                       obj.params(p).value) ;
            end
          end
          labBarrier() ;
        case 'tflow'
          % nothing to do
      end
    end
  end

  methods (Access = protected)
    function obj = startWithMMap(obj)
      obj.otherLabs = setdiff(1:numlabs, labindex) ;
      format = {} ;
      for i=1:numel(obj.params)
        format(i,1:3) = {obj.params(i).dataType, ...
                         obj.params(i).shape, ...
                         obj.params(i).name} ;
      end
      if labindex == 1
        f = fopen(obj.memoryMapFile,'wb') ;
        for g = 1:numlabs
          for i = 1:size(format,1)
            fwrite(f,zeros(format{i,2},format{i,1}),format{i,1}) ;
          end
        end
        fclose(f) ;
      end
      labBarrier() ;
      obj.memoryMap = memmapfile(obj.memoryMapFile, ...
                                 'Format', format, ...
                                 'Repeat', numlabs, ...
                                 'Writable', true) ;
      if obj.pinnedMemory && any([obj.params.isGPU])
        for i=1:numel(obj.memoryMap.Data)
          names = fieldnames(obj.memoryMap.Data(i))' ;
          for name = names
            vl_cudatool('cudaRegister', obj.memoryMap.Data(i).(char(name))) ;
          end
        end
      end
    end

    function stopWithMMap(obj)
      if ~isempty(obj.memoryMap)
        if obj.pinnedMemory && any([obj.params.isGPU])
          for i=1:numel(obj.memoryMap.Data)
            names = fieldnames(obj.memoryMap.Data(i))' ;
            for name = names
              vl_cudatool('cudaUnregister', obj.memoryMap.Data(i).(char(name))) ;
            end
          end
        end
        obj.memoryMap = [] ;
        if labindex == 1
          delete(obj.memoryMapFile) ;
        end
      end
    end

    function obj = startWithTFlow(obj)
      format = {} ;
      for i=1:numel(obj.params)
        format(i,1:4) = {obj.params(i).dataType, ...
                         obj.params(i).shape, ...
                         obj.params(i).name, ...
                         obj.params(i).deviceType} ;
      end
      vl_tflow('reset', obj.tflowOpts{:}) ;
      vl_tflow('init', format, ...
               labindex, numlabs, ...
               'prefix', obj.prefix, ...
               obj.tflowOpts{:}) ;
    end

    function stopWithTFlow(obj)
      vl_tflow('reset', obj.tflowOpts{:}) ;
    end

  end
end

repo: CREST-Release


text:
% plotAverageWaveforms
%
% Plots average waveforms for all channels and specified units across the
% array. It will prompt for a CMP file to map each channel to its electrode
% representation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use plotAverageWaveforms(NEV, units)
%
% NOTE: All input arguments are optional. Input arguments may be in any order.
%
%   NEV:          NEV structure containing all waveforms.
%                 DEFAULT: Will prompt for NEV.
%
%   'units':      Specify the units to plot (i.e. 2 or 0:5 for all).
%                 DEFAULT: will plot all units.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   USAGE EXAMPLE:
%
%   openNEV(NEV, 0:3);
%
%   In the example above, the waveforms in the NEV structure coming from
%   units 0 through 3 will be plotted.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Kian Torab
%   kian@blackrockmicro.com
%   Blackrock Microsystems
%   Salt Lake City, UT
%
%   Version 2.1.0.0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Version History
%
% 1.0.0.0:
%   - Initial release.
%
% 2.1.0.0:
%   - Added the ability to plot scaled or non-scaled plots. By default, the
%   plots are scaled.
%   - The color of the first unit is always selected at random, so when a
%   channel only has 1 unit, that 1 unit will be a different color
%   (visualization only).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function plotAverageWaveforms(NEV, units, plottype)

mapfile = 1;
colors = repmat({'magenta', 'cyan', 'red', 'green', 'blue', 'black'}, 1, 3);

if ~exist('NEV', 'var')
    NEV = openNEV('read');
end

if ~exist('units', 'var')
    units = 0:5;
end

if ~exist('plottype', 'var')
    plottype = 'scaled';
end

numberOfElectrodes = length(unique(NEV.Data.Spikes.Electrode));

disp('Select a map file for your array. If none, press cancel or ESC.');
mapFile = KTUEAMapFile;
if ~mapFile.isValid
    disp('Map File was not selected.');
    mapfile = 0;
    plotHeight = floor(sqrt(numberOfElectrodes));
    plotWidth = ceil(numberOfElectrodes/floor(sqrt(numberOfElectrodes)));
end

figure = KTFigure;

figure.EnlargeFigure;
figure.MakeBackgroundWhite;
figure.SetActive;
title('Spikes Average Waveforms for Individual Channels');

for chanIDX = 1:numberOfElectrodes

    if mapfile
try
        mapFile.GenerateChannelSubplot(chanIDX);
catch
    disp('e');
end
    else
        subplot(plotHeight, plotWidth, chanIDX);
    end

    hold on;
    initialColor = round(rand(1)*5);
    for unitIDX = units
        [~, allWaveforms] = findEventTimes(NEV, chanIDX, unitIDX);
        if ~isempty(allWaveforms)
            averageWaveform = mean(allWaveforms, 2);
            plot(averageWaveform, 'color', colors{unitIDX+1+initialColor}, 'LineWidth',2);
        end
    end
    hold off;

    if strcmpi(plottype, 'scaled')
        ylim([-1000, 1000]);
    end
%     if size(allWaveforms,2) ~= 0
%         text(1, 0, num2str(size(allWaveforms,2)), 'color', 'red');
%     end
    mapFile.setAxisBackgroundColor('white');
    set(gca, 'XTick', []);
    set(gca, 'YTick', []);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function p1 = gibbs(p1, x)
%GIBBS      Gibbs sampling
%
%           Description
%           P1 = GIBBS(P1, X) Samples one sample for parameters
%           contained in structure P1 at level n and and gets a
%           parameter structure X from level n-1 in hierarchy.
%           Returns a structure P1 containing new samples for
%           parameters in structure P1.
%
%           If P1 contains higher level parameters, e.g P1 =
%           net.p.w{1} and though there is a field P1.p, the
%           function calls itself to sample for those higher level
%           parameters also. So function takes care of the whole
%           hierarchical structure of ARD if called for the lowest
%           level parameters.

% Copyright (c) 1999-2003 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.

%Matlab 6.5
%if any(strcmp(fieldnames(p1),'p')) & ~isempty(p1.p) %isfield(p1,'p')
%  for a = fieldnames(p1.a)';
%    a=a{:};
%    if any(strcmp(fieldnames(p1.p),a)) %isfield(p1.p,a);
%      p2=p1.p.(a);
%      % call the spesified distribution sampler
%      p1.a.(a)=feval(['cond_' p1.f '_' p2.f],a,p1.a,p2.a,x);
%      % sample the upper level parameters
%      p1.p.(a)=gibbs(p2, p1.a.(a));
%    end
%  end
%end

% Pre matlab 6.5
S.type='.';
if any(strcmp(fieldnames(p1),'p')) & ~isempty(p1.p) %isfield(p1,'p')
  for a = fieldnames(p1.a)';
    a=a{:};
    S.subs=a;
    if any(strcmp(fieldnames(p1.p),a)) %isfield(p1.p,a);
      p2=subsref(p1.p,S);
      p1.a=subsasgn(p1.a, S, feval(['cond_' p1.f '_' p2.f],a,p1.a,p2.a,x));
      p1.p=subsasgn(p1.p, S, gibbs(p2, subsref(p1.a,S)));
    end
  end
end


text:
function [L,La,Lb] = points2plucker(A,B)

% POINTS2PLUCKER  Plucker line from two homogeneous points
%   L = POINTS2PLUCKER(A,B) is the Plucker line that passes over two points
%   A and B. The points are specified by their homogeneous coordinates
%   [x;y;z;w].
%
%   The result is a Plucker line expressed as a 6-vector. This vector can
%   be decomosed as L = [a;b], where the 3-vectors a and b admit the
%   following interpretation:
%
%   *   a - is a vector normal to the plane that contains the line and the
%       origin of coordinates.
%
%   *   b - is a vector director of the line, which lies on the plane
%       above.
%
%   *   a and b are orthogonal, i.e. dot(a,b)=0.
%
%   *   the distance from the line to the origin is given by
%       norm(a)/norm(b).
%
%   [L,La,Lb] = ... returns the Jacobians wrt A and B.
%
%   See also PLANES2PLUCKER.

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

L = [cross(A(1:3),B(1:3));A(4)*B(1:3)-B(4)*A(1:3)];

if nargout > 1

    [a1,a2,a3,a4] = split(A);
    [b1,b2,b3,b4] = split(B);

    La = [...
        [   0,  b3, -b2,   0]
        [ -b3,   0,  b1,   0]
        [  b2, -b1,   0,   0]
        [ -b4,   0,   0,  b1]
        [   0, -b4,   0,  b2]
        [   0,   0, -b4,  b3]];

    Lb = [...
        [   0, -a3,  a2,   0]
        [  a3,   0, -a1,   0]
        [ -a2,  a1,   0,   0]
        [  a4,   0,   0, -a1]
        [   0,  a4,   0, -a2]
        [   0,   0,  a4, -a3]];
end

return

%% jac

syms a1 a2 a3 a4 b1 b2 b3 b4 p1 p2 p3 p4 q1 q2 q3 q4 real
A=[a1;a2;a3;a4];
B=[b1;b2;b3;b4];
L=points2plucker(A,B)

La = jacobian(L,A)
Lb = jacobian(L,B)

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function lik = lik_inputdependentweibull(varargin)
%LIK_INPUTDEPENDENTWEIBULL  Create a (right censored) input dependent Weibull likelihood structure
%
%  Description
%    LIK = LIK_INPUTDEPENDENTWEIBULL('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a likelihood structure for a right censored input dependent
%    Weibull survival model in which the named parameters have the
%    specified values. Any unspecified parameters are set to default
%    values.
%
%    LIK = LIK_INPUTDEPENDENTWEIBULL(LIK,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a likelihood structure with the named parameters
%    altered with the specified values.
%
%    Parameters for Weibull likelihood [default]
%      shape       - shape parameter r [1]
%      shape_prior - prior for shape [prior_logunif]
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%    The likelihood is defined as follows:
%                       __ n
%      p(y|f1, f2, z) = || i=1 [ r_i^(1-z_i) exp( (1-z_i)*(-f1_i)
%                           +(1-z_i)*(r_i-1)*log(y_i)
%                           -exp(-f1_i)*y_i^r_i) ]
%
%    where r_i=r*exp(f2_i) is the shape parameter of Weibull distribution.
%    z is a vector of censoring indicators with z = 0 for uncensored event
%    and z = 1 for right censored event. Here the second latent variable f2
%    implies the input dependance to the shape parameter in the original
%    Weibull likelihood.
%
%    When using the Weibull likelihood you can give the
%    vector z as an extra parameter to each function that requires
%    also y. For example, you can call gp_optim as follows:
%      gp_optim(gp, x, y, 'z', z)
%    If z is not given or it is empty, then usual likelihood for
%    uncensored data is used
%
%  See also
%    GP_SET, LIK_*, PRIOR_*
%
% Copyright (c) 2011 Jaakko Riihimäki
% Copyright (c) 2011 Aki Vehtari
% Copyright (c) 2012 Ville Tolvanen

% 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.

  ip=inputParser;
  ip.FunctionName = 'LIK_INPUTDEPENDENTWEIBULL';
  ip.addOptional('lik', [], @isstruct);
  ip.addParamValue('shape',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('shape_prior',prior_logunif(), @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.nondiagW=true;
    lik.type = 'Inputdependent-Weibull';
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'Weibull')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('shape',ip.UsingDefaults)
    lik.shape = ip.Results.shape;
  end
  % Initialize prior structure
  if init
    lik.p=[];
  end
  if init || ~ismember('shape_prior',ip.UsingDefaults)
    lik.p.shape=ip.Results.shape_prior;
  end

  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_inputdependentweibull_pak;
    lik.fh.unpak = @lik_inputdependentweibull_unpak;
    lik.fh.lp = @lik_inputdependentweibull_lp;
    lik.fh.lpg = @lik_inputdependentweibull_lpg;
    lik.fh.ll = @lik_inputdependentweibull_ll;
    lik.fh.llg = @lik_inputdependentweibull_llg;
    lik.fh.llg2 = @lik_inputdependentweibull_llg2;
    lik.fh.llg3 = @lik_inputdependentweibull_llg3;
    lik.fh.invlink = @lik_inputdependentweibull_invlink;
    lik.fh.predy = @lik_inputdependentweibull_predy;
    lik.fh.recappend = @lik_inputdependentweibull_recappend;
  end

end

function [w,s,h] = lik_inputdependentweibull_pak(lik)
%LIK_INPUTDEPENDENTWEIBULL_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_INPUTDEPENDENTWEIBULL_PAK(LIK) takes a likelihood structure LIK and
%    combines the parameters into a single row vector W. This is a
%    mandatory subfunction used for example in energy and gradient
%    computations.
%
%       w = log(lik.shape)
%
%   See also
%   LIK_INPUTDEPENDENTWEIBULL_UNPAK, GP_PAK

  w=[];s={};h=[];
  if ~isempty(lik.p.shape)
    w = log(lik.shape);
    s = [s; 'log(weibull.shape)'];
    h = [h 0];
    [wh, sh, hh] = lik.p.shape.fh.pak(lik.p.shape);
    w = [w wh];
    s = [s; sh];
    h = [h hh];
  end
end

function [lik, w] = lik_inputdependentweibull_unpak(lik, w)
%LIK_INPUTDEPENDENTWEIBULL_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    [LIK, W] = LIK_INPUTDEPENDENTWEIBULL_UNPAK(W, LIK) takes a likelihood
%    structure LIK and extracts the parameters from the vector W
%    to the LIK structure. This is a mandatory subfunction used
%    for example in energy and gradient computations.
%
%   Assignment is inverse of
%       w = log(lik.shape)
%
%   See also
%   LIK_INPUTDEPENDENTWEIBULL_PAK, GP_UNPAK

  if ~isempty(lik.p.shape)
    lik.shape = exp(w(1));
    w = w(2:end);
    [p, w] = lik.p.shape.fh.unpak(lik.p.shape, w);
    lik.p.shape = p;
  end
end

function lp = lik_inputdependentweibull_lp(lik, varargin)
%LIK_INPUTDEPENDENTWEIBULL_LP  log(prior) of the likelihood parameters
%
%  Description
%    LP = LIK_INPUTDEPENDENTWEIBULL_LP(LIK) takes a likelihood structure LIK and
%    returns log(p(th)), where th collects the parameters. This
%    subfunction is needed when there are likelihood parameters.
%
%  See also
%    LIK_INPUTDEPENDENTWEIBULL_LLG, LIK_INPUTDEPENDENTWEIBULL_LLG3, LIK_INPUTDEPENDENTWEIBULL_LLG2, GPLA_E

% If prior for shape parameter, add its contribution
  lp=0;
  if ~isempty(lik.p.shape)
    lp = lik.p.shape.fh.lp(lik.shape, lik.p.shape) +log(lik.shape);
  end

end

function lpg = lik_inputdependentweibull_lpg(lik)
%LIK_INPUTDEPENDENTWEIBULL_LPG  d log(prior)/dth of the likelihood
%                parameters th
%
%  Description
%    E = LIK_INPUTDEPENDENTWEIBULL_LPG(LIK) takes a likelihood structure LIK and
%    returns d log(p(th))/dth, where th collects the parameters.
%    This subfunction is needed when there are likelihood parameters.
%
%  See also
%    LIK_INPUTDEPENDENTWEIBULL_LLG, LIK_INPUTDEPENDENTWEIBULL_LLG3, LIK_INPUTDEPENDENTWEIBULL_LLG2, GPLA_G

  lpg=[];
  if ~isempty(lik.p.shape)
    % Evaluate the gprior with respect to shape
    ggs = lik.p.shape.fh.lpg(lik.shape, lik.p.shape);
    lpg = ggs(1).*lik.shape + 1;
    if length(ggs) > 1
      lpg = [lpg ggs(2:end)];
    end
  end
end

function ll = lik_inputdependentweibull_ll(lik, y, ff, z)
%LIK_INPUTDEPENDENTWEIBULL_LL  Log likelihood
%
%  Description
%    LL = LIK_INPUTDEPENDENTWEIBULL_LL(LIK, Y, F, Z) takes a likelihood
%    structure LIK, survival times Y, censoring indicators Z, and
%    latent values F. Returns the log likelihood, log p(y|f,z).
%    This subfunction is needed when using Laplace approximation
%    or MCMC for inference with non-Gaussian likelihoods. This
%    subfunction is also used in information criteria (DIC, WAIC)
%    computations.
%
%  See also
%    LIK_INPUTDEPENDENTWEIBULL_LLG, LIK_INPUTDEPENDENTWEIBULL_LLG3, LIK_INPUTDEPENDENTWEIBULL_LLG2, GPLA_E

  if numel(z)==0
    % no censoring by default
    z=0;
  end

  f=ff(:);
  n=size(y,1);
  f1=f(1:n);
  f2=f((n+1):2*n);
  expf2=exp(f2);
  expf2(isinf(expf2))=realmax;
  a = lik.shape;
  ll = sum((1-z).*(log(a*expf2) + (a*expf2-1).*log(y)-f1) - exp(-f1).*y.^(a*expf2));

end

function llg = lik_inputdependentweibull_llg(lik, y, ff, param, z)
%LIK_INPUTDEPENDENTWEIBULL_LLG  Gradient of the log likelihood
%
%  Description
%    LLG = LIK_INPUTDEPENDENTWEIBULL_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, survival times Y, censoring indicators Z and
%    latent values F. Returns the gradient of the log likelihood
%    with respect to PARAM. At the moment PARAM can be 'param' or
%    'latent'. This subfunction is needed when using Laplace
%    approximation or MCMC for inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_INPUTDEPENDENTWEIBULL_LL, LIK_INPUTDEPENDENTWEIBULL_LLG2, LIK_INPUTDEPENDENTWEIBULL_LLG3, GPLA_E

  if numel(z)==0
    % no censoring by default
    z=0;
  end

  f=ff(:);
  n=size(y,1);
  f1=f(1:n);
  f2=f((n+1):2*n);
  expf2=exp(f2);
  expf2(isinf(expf2))=realmax;
  a = lik.shape;
  switch param
    case 'param'
      llg = sum((1-z).*(1./a + expf2.*log(y)) - exp(-f1).*y.^(a.*expf2).*log(y).*expf2);
      % correction for the log transformation
      llg = llg.*lik.shape;
    case 'latent'
      llg1 = -(1-z) + exp(-f1).*y.^(a.*expf2);
      llg2 = (1-z).*(1 + a.*expf2.*log(y)) - exp(-f1).*y.^(a.*expf2).*log(y).*a.*expf2;
      llg = [llg1; llg2];
  end
end

function llg2 = lik_inputdependentweibull_llg2(lik, y, ff, param, z)
%LIK_INPUTDEPENDENTWEIBULL_LLG2  Second gradients of the log likelihood
%
%  Description
%    LLG2 = LIK_INPUTDEPENDENTWEIBULL_LLG2(LIK, Y, F, PARAM) takes a
%    likelihood structure LIK, survival times Y, censoring indicators Z,
%    and latent values F. Returns the hessian of the log likelihood with
%    respect to PARAM. LLG2 is a (2*length(y)) x 2 matrix containing the
%    non-zero elements of the Hessian matrix. This subfunction is needed
%    when using Laplace approximation or EP for  inference with
%    non-Gaussian likelihoods.
%
%  See also
%    LIK_INPUTDEPENDENTWEIBULL_LL, LIK_INPUTDEPENDENTWEIBULL_LLG,
%    LIK_INPUTDEPENDENTWEIBULL_LLG3, GPLA_E

  if numel(z)==0
    % no censoring by default
    z=0;
  end

  a = lik.shape;
  f=ff(:);
  n=size(y,1);
  f1=f(1:n);
  f2=f((n+1):2*n);
  expf2=exp(f2);
  expf2(isinf(expf2))=realmax;
  switch param
    case 'param'

    case 'latent'
      t1=exp(-f1).*y.^(a.*expf2);
      t2=log(y).*a.*expf2;
      t3=t1.*t2;

      llg2_11 = -t1;
      llg2_12 = t3;
      llg2_22 = (1-z).*t2 - (t2 + 1).*t3;

      llg2 = [llg2_11 llg2_12; llg2_12 llg2_22];

    case 'latent+param'
      t1=expf2.*log(y);
      t2=exp(-f1).*y.^(a.*expf2);
      t3=t1.*t2;
      llg2_1 = t3;
      llg2_2 = (1-z).*t1 - (t1.*a + 1).*t3;
      llg2 = [llg2_1; llg2_2];
      % correction due to the log transformation
      llg2 = llg2.*lik.shape;
  end
end

function llg3 = lik_inputdependentweibull_llg3(lik, y, ff, param, z)
%LIK_INPUTDEPENDENTWEIBULL_LLG3  Third gradients of the log likelihood
%
%  Description
%    LLG3 = LIK_INPUTDEPENDENTWEIBULL_LLG3(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, survival times Y, censoring indicators Z and
%    latent values F and returns the third gradients of the log
%    likelihood with respect to PARAM. At the moment PARAM can be
%    only 'latent'. LLG3 is a vector with third gradients. This
%    subfunction is needed when using Laplace approximation for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_INPUTDEPENDENTWEIBULL_LL, LIK_INPUTDEPENDENTWEIBULL_LLG, LIK_INPUTDEPENDENTWEIBULL_LLG2, GPLA_E, GPLA_G

  if numel(z)==0
    % no censoring by default
    z=0;
  end

  a = lik.shape;
  f=ff(:);
  n=size(y,1);
  f1=f(1:n);
  f2=f((n+1):2*n);
  expf2=exp(f2);
  expf2(isinf(expf2))=realmax;
  switch param
    case 'param'

    case 'latent'
      t1=a.*expf2.*log(y);
      t2=exp(-f1).*y.^(a.*expf2);
      t3=t2.*t1;
      t4=t3.*t1;

      nl=2;
      llg3=zeros(nl,nl,nl,n);

      llg3(1,1,1,:) = t2;

      llg3(2,2,1,:) = t4 + t3;
      llg3(2,1,2,:) = llg3(2,2,1,:);
      llg3(1,2,2,:) = llg3(2,2,1,:);

      llg3(2,1,1,:) = -t3;
      llg3(1,2,1,:) = llg3(2,1,1,:);
      llg3(1,1,2,:) = llg3(2,1,1,:);

      llg3(2,2,2,:) = (1-z).*t1 - t4.*t1 - 3.*t4 - t3;

    case 'latent2+param'
      t1 = log(y).*expf2;
      t2 = exp(-f1).*y.^(a*expf2);
      t3 = t2.*t1;
      t4 = t3.*t1;
      llg3_11 = -t3;
      llg3_12 = a.*t4 + t3;
      llg3_22 = (1-z).*t1 - a.^2.*t4.*t1 - 3.*a.*t4 - t3;

      llg3 = [diag(llg3_11) diag(llg3_12); diag(llg3_12) diag(llg3_22)];
      % correction due to the log transformation
      llg3 = llg3.*lik.shape;
  end
end

function [lpy, Ey, Vary] = lik_inputdependentweibull_predy(lik, Ef, Varf, yt, zt)
%LIK_INPUTDEPENDENTWEIBULL_PREDY  Returns the predictive mean, variance and density of y
%
%  Description
%    LPY = LIK_INPUTDEPENDENTWEIBULL_PREDY(LIK, EF, VARF YT, ZT)
%    Returns logarithm of the predictive density PY of YT, that is
%        p(yt | zt) = \int p(yt | f, zt) p(f|y) df.
%    This requires also the survival times YT, censoring indicators ZT.
%    This subfunction is needed when computing posterior predictive
%    distributions for future observations.
%
%    [LPY, EY, VARY] = LIK_INPUTDEPENDENTWEIBULL_PREDY(LIK, EF, VARF) takes a
%    likelihood structure LIK, posterior mean EF and posterior
%    Variance VARF of the latent variable and returns the
%    posterior predictive mean EY and variance VARY of the
%    observations related to the latent variables. This subfunction
%    is needed when computing posterior predictive distributions for
%    future observations.
%
%
%  See also
%    GPLA_PRED, GPEP_PRED, GPMC_PRED

  if numel(zt)==0
    % no censoring
    zt=zeros(size(yt));
  end

  yc = 1-zt;
  r = lik.shape;

  Ef = Ef(:);
  ntest = 0.5*size(Ef,1);
  Ef1=Ef(1:ntest); Ef2=Ef(ntest+1:end);
  %     Varf1=squeeze(Varf(1,1,:)); Varf2=squeeze(Varf(2,2,:));
  if size(Varf,2) == size(Varf,1)
    Varf1=diag(Varf(1:ntest,1:ntest));Varf2=diag(Varf(ntest+1:end,ntest+1:end));
  else
    Varf1=Varf(:,1); Varf2=Varf(:,2);
  end

  Ey=[];
  Vary=[];

  % Evaluate the posterior predictive densities of the given observations
  lpy = zeros(length(yt),1);

  for i2=1:ntest
    m1=Ef1(i2); m2=Ef2(i2);
    s1=sqrt(Varf1(i2)); s2=sqrt(Varf2(i2));
    % Function handle for Weibull * Gaussian_f1 * Gaussian_f2
    pd=@(f1,f2) exp(yc(i2).*((log(r) + f2) + (r.*exp(f2)-1).*log(yt(i2))-f1) - exp(-f1).*yt(i2).^(r*exp(f2))) ...
      .*norm_pdf(f1,Ef1(i2),sqrt(Varf1(i2))).*norm_pdf(f2,Ef2(i2),sqrt(Varf2(i2)));
    % Integrate over latent variables
    lpy(i2) = log(dblquad(pd, m1-6.*s1, m1+6.*s1, m2-6.*s2, m2+6.*s2));
  end

end

function p = lik_inputdependentweibull_invlink(lik, f, z)
%LIK_INPUTDEPENDENTWEIBULL Returns values of inverse link function
%
%  Description
%    P = LIK_INPUTDEPENDENTWEIBULL_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values of inverse link function P.
%    This subfunction is needed when using function gp_predprctmu.
%
%     See also
%     LIK_INPUTDEPENDENTWEIBULL_LL, LIK_INPUTDEPENDENTWEIBULL_PREDY

p = exp(f);
end

function reclik = lik_inputdependentweibull_recappend(reclik, ri, lik)
%RECAPPEND  Append the parameters to the record
%
%  Description
%    RECLIK = LIK_INPUTDEPENDENTWEIBULL__RECAPPEND(RECLIK, RI, LIK) takes a
%    likelihood record structure RECLIK, record index RI and
%    likelihood structure LIK with the current MCMC samples of
%    the parameters. Returns RECLIK which contains all the old
%    samples and the current samples from LIK. This subfunction
%    is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC

  if nargin == 2
    % Initialize the record
    reclik.type = 'Inputdependent-Weibull';
    reclik.nondiagW=true;

    % Initialize parameter
%     reclik.shape = [];

    % Set the function handles
    reclik.fh.pak = @lik_inputdependentweibull_pak;
    reclik.fh.unpak = @lik_inputdependentweibull_unpak;
    reclik.fh.lp = @lik_t_lp;
    reclik.fh.lpg = @lik_t_lpg;
    reclik.fh.ll = @lik_inputdependentweibull_ll;
    reclik.fh.llg = @lik_inputdependentweibull_llg;
    reclik.fh.llg2 = @lik_inputdependentweibull_llg2;
    reclik.fh.llg3 = @lik_inputdependentweibull_llg3;
    reclik.fh.invlink = @lik_inputdependentweibull_invlink;
    reclik.fh.predy = @lik_inputdependentweibull_predy;
    reclik.fh.recappend = @lik_inputdependentweibull_recappend;
    reclik.p=[];
    reclik.p.shape=[];
    if ~isempty(ri.p.shape)
      reclik.p.shape = ri.p.shape;
    end
  else
    % Append to the record
    reclik.shape(ri,:)=lik.shape;
    if ~isempty(lik.p.shape)
      reclik.p.shape = lik.p.shape.fh.recappend(reclik.p.shape, ri, lik.p.shape);
    end
  end
end

author: gpstuff-dev

repo: gpstuff

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


text:
clear all;
close all;
clc;
spx.cluster.ssc.util.bench_subspace_preservation(@ssc_nn_omp, 'ssc_nn_omp');

author: indigits

repo: sparse-plex


text:
function v = cnn_vecms_sketch(im, net, emodel, scales, mirror, aggregate)
% CNN_VECMS_SKETCH computes a multi-scale D-dimensional CNN sketch vector for an image.
%
%   V = cnn_vecms_sketch(IM, NET, EMODEL, SCALES, MIRROR, AGGREGATE)
%
%   IM        : Input image, or input image path as a string
%   NET       : CNN network to evaluate on image IM
%   EMODEL    : Edge detector model, 0 means do not extract edges because input is already sketch
%   SCALES    : Vector of scales to resize image prior to the descriptor computation (DEFAULT: [1])
%   MIRROR    : Mirror (horizontal flip) of image for the descriptor computation (DEFAULT: 0)
%   AGGREGATE : Aggregate over scales and mirror if 1, otherwise return one vector per image scale and mirror (DEFAULT: 1)
%
%   V         : Multi-scale output, D x 1 if AGGREGATE=1, D x numel(SCALES) x (MIRROR+1) if AGGREGATE=0
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2018.

    if ~exist('scales'), scales = 1; end
    if ~exist('mirror'), mirror = 0; end
    if ~exist('aggregate'), aggregate = 1; end

    padsize = 30;

    v = [];
    for i = 1:(mirror+1)*numel(scales)
        if i == numel(scales)+1, im = fliplr(im); end;
        imr = imresize(im, scales(mod(i-1,numel(scales))+1));

        if ~isstruct(emodel)
            %% Treated as SKETCH
            if size(imr, 3) > 1; imr = rgb2gray(imr);  end
            % making sure its a binary sketch image
            imr = single(imr); imr = 255 * ((imr - min(imr(:))) ./ (max(imr(:)) - min(imr(:))));
            imr = single(imr < 0.8 * max(imr(:)));
            imr = padarray(imr, [padsize padsize], 0, 'both');
            % unifying sketches:
            % morphological thinning followed by dilation
            imr = single(bwmorph(imr, 'thin', inf));
            imr = single(imdilate(imr,strel('disk',1,0)));
            imr = single(bwmorph(imr, 'thin', inf));
            imr = single(imdilate(imr,strel('disk',2,0)));
            % final sketch2edgemap
            e = imr;
        else
            %% Treated as IMAGE
            if size(imr, 3) == 1; imr = repmat(imr, 1, 1, 3); end
            [edg, ori] = edgesDetect(imr, emodel);
            e = padarray(edg, [padsize padsize], 0, 'both');
        end

        v = [v, cnn_vec_edge(e, net)];
    end

    if aggregate
        v = sum(v, 2);
        v = vecpostproc(v);
    end

function v = cnn_vec_edge(im, net)
% CNN_VEC_EDGE computes a D-dimensional CNN edge vector for a given edge-map.
%
%   V = cnn_vec_edge(IM, NET)
%
%   IM  : Input edge-map
%   NET : CNN network to evaluate on image IM
%
%   V   : D-dimensional vector
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2018.

    minsize = 67;
    net.mode = 'test';
    descvarname = 'l2descriptor';

    use_gpu = strcmp(net.device, 'gpu');
    if use_gpu
        gpuarrayfun = @(x) gpuArray(x);
        gatherfun = @(x) gather(x);
    else
        gpuarrayfun = @(x) x; % do not convert to gpuArray
        gatherfun = @(x) x; % do not gather
    end

    if isstr(im)
        im = imread(im);
    end

    im = single(im) - mean(net.meta.normalization.averageImage(:));
    if min(size(im, 1), size(im, 2)) < minsize
        im = pad2minsize(im, minsize, 0);
    end

    net.eval({'input', gpuarrayfun(reshape(im, [size(im), 1]))});

    v = gatherfun(squeeze(net.getVar(descvarname).value));

author: filipradenovic

repo: cnnimageretrieval


text:
classdef (InferiorClasses = {?rotation,?quaternion}) DelaunaySO3 < orientation
%
% Syntax
%
%   % define a Delaunay trinangulation from a list of orientations
%   DSO3 = DelaunaySO3(ori)
%
  properties

    tetra          % list of vertices of the tetrahegons
    tetraNeighbour % neigbouring tetraeders orderd as faces
    lookup         % lookup table orientation -> tetrahedrons

  end

  methods

    function DSO3 = DelaunaySO3(varargin)

      DSO3 = DSO3@orientation(varargin{:});

      if nargin == 0, return; end

      % pertube data a bit
      % it would be better if this would be needed only for the
      % triangulation, but for some reason we cant skip it
      [DSO3.a,DSO3.b,DSO3.c,DSO3.d] = double(perturbe(DSO3,0.05*degree));

      % compute tetrahegons
      DSO3.tetra = calcDelaunay2(DSO3);

      % compute neighbouring list
      DSO3.tetraNeighbour = calcNeighbour(DSO3.tetra);

      % compute lookup table
      res = 40*degree;
      for i = 1:4
        res = res / 2;
        DSO3.lookup = calcLookUp(DSO3,res);
      end
    end

end

end

author: mtex-toolbox

repo: mtex

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


text:
%
% Copyright (c) 2015, Yarpiz (www.yarpiz.com)
% All rights reserved. Please read the "license.txt" for license terms.
%
% Project Code: YPML110
% Project Title: Implementation of DBSCAN Clustering in MATLAB
% Publisher: Yarpiz (www.yarpiz.com)
%
% Developer: S. Mostapha Kalami Heris (Member of Yarpiz Team)
%
% Contact Info: sm.kalami@gmail.com, info@yarpiz.com
%

function PlotClusterinResult(X, IDX)
    nDim = size(X, 2);
    k=max(IDX);

    Colors=hsv(k);

    Legends = {};
    for i=0:k
        Xi=X(IDX==i,:);
        if i~=0
            Style = 'x';
            MarkerSize = 8;
            Color = Colors(i,:);
            Legends{end+1} = ['Cluster #' num2str(i)];
        else
            Style = 'o';
            MarkerSize = 6;
            Color = [0 0 0];
            if ~isempty(Xi)
                Legends{end+1} = 'Noise';
            end
        end
        if ~isempty(Xi)
            str = sprintf('Id=%02d, Len=%02d', i-1, length(Xi(:, 1)));
            pt = Xi(1, :);
            if nDim == 2
                text(pt(1), pt(2), str, 'FontSize', 12 );
                plot(Xi(:,1),Xi(:,2),Style,'MarkerSize',MarkerSize,'Color',Color);
            else
                text(pt(1), pt(2), pt(3), str, 'FontSize', 12 );
                plot3(Xi(:,1),Xi(:,2),Xi(:, 3), Style,'MarkerSize',MarkerSize,'Color',Color);
            end
        end
        hold on;
    end
    hold off;
    axis equal;
    grid on;
    % legend(Legends);
    % legend('Location', 'NorthEastOutside');

end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function geneScoreStructure=mapCompartments(geneScoreStructure,varargin)
% mapCompartments
%   Maps compartments in the geneScoreStructure. This is used if you do not
%   want a models that uses all of the compartment from the predictor. This
%   function will then let you define rules on how the compartments should
%   be merged.
%
%   Any number of rules could be defined as consecutive strings or in a cell array.
%   'comp1'             comp1 should be kept in the structure

%   'comp1=comp2'       The scores in comp2 are merged to comp1 and comp2 is
%                       removed from the structure. This automatically
%                       keeps comp1 in the structure

%   'comp1=comp2 comp3' The scores in comp2 and comp3 are merged to comp1
%                       and comp2 & comp2 are removed from the structure.
%                       This automatically keeps comp1 in the structure

%   'comp1 comp2=comp3' The scores in comp3 are split between comp1 and
%                       comp2. This automatically keeps comp1 and comp2 in
%                       the structure

%   'comp1=other'       The scores in any compartment not included are
%                       merged to comp1. This is applied after all other
%                       rules.
%
%   When one compartment is merged to another the resulting scores will be
%   the best for each gene in either of the compartments. In the case where
%   one compartment is split among several, the scores for the compartment
%   to be merged is weighted with the number of compartments to split to.
%
%   Example: The predictor you use give prediction for Extracellular,
%   Cytosol, Nucleus, Peroxisome, Mitochondria, ER, and Lysosome. You want to
%   have a model with Extracellular, Cytosol, Mitochondria, and Peroxisome
%   where Lysosome is merged with Peroxisome and all other compartments
%   are merged to the Cytosol.
%
%   GSS=mapCompartments(GSS,'Extracellular','Mitochondria','Peroxisome=Lyso
%   some','Cytosol=other');
%
%   geneScoreStructure  a structure to be used in predictLocalization
%
%   Usage: geneScoreStructure=mapCompartments(geneScoreStructure,varargin)

varargin=upper(varargin);

%First find the compartment that will end up in the final structure. They
%are the ones that stand alone or are to the left of some '='
toKeep={};
toMerge={};
I=regexp(varargin,'=','split');
for i=1:numel(varargin)
    if numel(I{i})==1
        toKeep=[toKeep;I{i}];
    else
        J=regexp(I{i}(1),' ','split');
        K=regexp(I{i}(2),' ','split');
        toKeep=[toKeep;J{1}(:)];
        toMerge=[toMerge;K{1}(:)];
    end
end

%Check that there are no compartment that should both be merged and kept
if ~isempty(intersect(toKeep,toMerge))
    EM='There are inconsistencies where one or more compartment(s) should be both kept and merged to another';
    dispEM(EM);
end

%Check that there are no compartments in the rules that are not in the
%geneScoreStructure.
uComps=upper(geneScoreStructure.compartments);
J=[uComps;{'OTHER'}];

if ~isempty(setdiff([toKeep;toMerge],J))
    EM='There are compartment in the rules that are not in geneScoreStructure.compartments';
    dispEM(EM);
end

%Loop through it again and do the mapping
otherIndex=[]; %This stores the rule which maps 'other'.

for i=1:numel(I)
    if numel(I{i})>1
        %Get the compartment indexes that should be mapped
        J=regexp(I{i}(2),' ','split');
        if strcmpi(J{1},'other')
            otherIndex=i;
            continue;
        end
        [k, K]=ismember(J{1},uComps);

        %And to where they should be mapped
        J=regexp(I{i}(1),' ','split');
        [l, L]=ismember(J{1},uComps);

        %It's not allowed to have rules like A B=C D
        if numel(K)>1 && numel(L)>1
            EM='It is not allowed to have rules like "A B=C D" (map more than one compartment to more than one compartment)';
            dispEM(EM);
        end

        if ~all(k) || ~all(l)
            EM='Error in mapping. This most likely means that some compartment(s) are mapped to different compartments in different rules. Use A B=C if you want to map C to several compartments';
            dispEM(EM);
        end

        %Get the sum of the scores for the compartments that should be
        %merged to something else
        S=max(geneScoreStructure.scores(:,K),[],2);
        for j=1:numel(L)
            %If the scores are mapped to several different compartments
            %then split the scores between them
            geneScoreStructure.scores(:,L(j))=max(geneScoreStructure.scores(:,L(j)),S./numel(L));
        end

        %Remove the comparement that were merged
        geneScoreStructure.compartments(K)=[];
        geneScoreStructure.scores(:,K)=[];
        uComps(K)=[];
    end
end

%Then check if there are remaining compartments that should be removed or
%mapped as 'other'
J=find(~ismember(uComps,toKeep));
if any(J)
    if any(otherIndex)
        K=regexp(I{otherIndex}(1),' ','split');
        [l, L]=ismember(K{1},uComps);
        if l==1 && numel(l)==1
            S=max(geneScoreStructure.scores(:,J),[],2);
            geneScoreStructure.scores(:,L)=max(geneScoreStructure.scores(:,L),S);
        else
            EM='Could not map "other" to more than one compartment';
            dispEM(EM);
        end
    else
        EM='There are compartments that are not defined if they should be kept or removed. Use "A=other" or define more rules if you do not want them to be deleted';
        dispEM(EM,false);
    end

    %Remove the comparement that were merged
    geneScoreStructure.compartments(J)=[];
    geneScoreStructure.scores(:,J)=[];
end

%Renormalize
I=max(geneScoreStructure.scores,[],2);
geneScoreStructure.scores=bsxfun(@times, geneScoreStructure.scores, 1./I);

%If there are genes that have score 0 in all compartments, remove them and
%print a warning.
I=find(isnan(geneScoreStructure.scores(:,1))); %Only looks a the first colum as it will be the same for the other ones
if any(I)
    EM='The following genes had score 0.0 in all compartments. They have been removed from the structure. Consider using more rules or "A=other" in order to prevent this:';
    dispEM(EM,false,geneScoreStructure.genes(I));
    geneScoreStructure.scores(I,:)=[];
    geneScoreStructure.genes(I)=[];
end
end

author: SysBioChalmers

repo: RAVEN


text:
function l = luma(img)
%
%       l = luma(img)
%
%       This function calculates the luma
%
%
%       input:
%           img: an RGB image
%
%       output:
%           l: luminance as XYZ color
%
%     Copyright (C) 2016  Francesco Banterle
%
%     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/>.
%

col = size(img, 3);

switch col
    case 1
        l = img;

    case 3
        l = 0.299 * img(:,:,1) + 0.587 * img(:,:,2) + 0.114 * img(:,:,3);

    otherwise
        l = mean(img, 3);
        disp('Mean of channels was computed; the input image is not an RGB or luminance image!');
end

end

author: banterle


text:
function test_nmtf(varargin)
% function test_nmtf()
%
% demonstration file for NMFLibrary.
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on June 14, 2022

    if nargin < 1
        clc;
        clear;
        close all;
        rng('default')

        m = 500;
        n = 100;
        V = rand(m,n);
        rank = 20;
        options = [];
        options.verbose = 2;
        options.max_epoch = 100;
        health_check_mode = false;
    else
        V = varargin{1};
        rank = varargin{2};
        options = varargin{3};
        health_check_mode = true;
    end

    %% perform factroization
    [x, info] = sep_symm_nmtf(V, rank, options);

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function p = predict(theta, X)
%PREDICT Predict whether the label is 0 or 1 using learned logistic
%regression parameters theta
%   p = PREDICT(theta, X) computes the predictions for X using a
%   threshold at 0.5 (i.e., if sigmoid(theta'*x) >= 0.5, predict 1)

m = size(X, 1); % Number of training examples

% You need to return the following variables correctly
p = zeros(m, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters.
%               You should set p to a vector of 0's and 1's
%

p = round(sigmoid(X * theta))

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

end

author: Ayatans

repo: Machine-Learning-homework


text:
%% OpenMAS MONTE-CARLO SIMULATION SETUP (monteCarlo_setup.m) %%%%%%%%%%%%%%
% This function is designed to build a monte-carlo style analysis of a
% given scenario.

% Author: James A. Douthwaite 26/03/18

% ADD THE PROGRAM PATHS
clear all; close all;
addpath('environment');
addpath('objects');
addpath('scenarios');
addpath('toolboxes/Intlab_V7.1');

%If intlab needs to be reloaded
try
    wrkDir = pwd;
    test = infsup(0,1);
    clearvars test
    IntDir = strcat(pwd,'\Intlab_V7.1\startupJD.m');
catch
    run('startup.m');
    cd(wrkDir)
end

%% //////////////////// AGENT & SCENARIO PARAMETERS ///////////////////////
fprintf('[SETUP]\tInitialising Monte-Carlo setup script.\n');

%% SIMULATION PARAMETERS
[~, userdir]   = system('echo %USERPROFILE%'); % Get desktop path
sim_outputPath = strcat(userdir,'\desktop\OpenMAS_data');
sim_warningDistance = 10;
sim_threadPool = 0;
sim_maxDuration = 30;
sim_timeStep = 0.25;    % RVO timestep (0.25s)
sim_vebosity = 1;
% sim_figureSet ={'all'};
sim_figureSet = {'plan','closest','inputs','gif'}; %{'fig','gif','inputs','separations'};

% SCENARIO PARAMETERS
sim_plotScenario = 1;
sim_agentNumber = 8;    % TOTAL NUMBER OF AGENTS
sim_agentOrbit  = 10;
sim_waypointOrbit = 10;
sim_agentVelocity = 0;
sim_noiseSigma = 0.2;

% ////////////////////////// INITIALISE AGENTS ////////////////////////////
fprintf('[SETUP]\tAssigning agent definitions:\n');
for index = 1:sim_agentNumber
    % BASIC CLASSES
%     agentIndex{index} = agent_VO();
%     agentIndex{index} = agent_RVO();
%     agentIndex{index} = agent_HRVO();
%     agentIndex{index} = agent_vectorSharing();
%     agentIndex{index} = agent_vectorSharing_interval();

%     agentIndex{index} = agent_2D_VO();
    agentIndex{index} = agent_2D_RVO();
%     agentIndex{index} = agent_2D_HRVO();
%     agentIndex{index} = agent_2D_RVO2();
%     agentIndex{index} = agent_2D_vectorSharing();
%     agentIndex{index} = agent_2D_vectorSharing_interval();
end

fprintf('[SETUP]\t Intialising the agents in the given scenario.\n');
% CALL THE FUNCTION THAT GENERATES THE UKACC SCENARIO
[ objectIndex ] = getScenario_UKACC_concentric('agents',agentIndex,'agentOrbit',sim_agentOrbit,'agentVelocity',sim_agentVelocity,'waypointOrbit',sim_waypointOrbit,'offsetAngle',pi/2,'plot',sim_plotScenario);

% ////////////// INITIALISE THE SIMULATION WITH THE OBJECT INDEX /////////////////
[DATA,META] = OMAS_initialise('objects',objectIndex,...
                             'duration',sim_maxDuration,...
                                   'dt',sim_timeStep,...
                              'figures',sim_figureSet,...
                      'warningDistance',sim_warningDistance,...
                            'verbosity',sim_vebosity,...
                           'threadPool',sim_threadPool,...
                           'outputPath',sim_outputPath);

clearvars -except DATA META
load(strcat(META.outputPath,'META.mat'));
load(strcat(META.outputPath,'EVENTS.mat'));

% ===================== MONTE-CARLO SIMULATION ============================
run_monteCarlo = false;
if ~run_monteCarlo
   return
end

fprintf('[SETUP]\tInitialising Monte-Carlo setup script.\n');

% /////////////////// AGENT & SCENARIO PARAMETERS /////////////////////////
sim_agentOrbit  = 30;       % [FIXED]
sim_agentVelocity = 0;      % [FIXED]
sim_waypointOrbit = 35;     % [FIXED]

% //////////////////////// OMAS CONFIGURATION /////////////////////////////
sim_warningDistance = 10;   % [FIXED]
sim_maxDuration = 120;      % [FIXED]
sim_timeStep = 0.25;     	% [FIXED]
sim_verbosity = 0;          % [FIXED]

% ///////////////////// MONTE-CARLO CONFIGURATION /////////////////////////
monteCarlo_algs = {'agent_2D_VO',...
                   'agent_2D_RVO',...
                   'agent_2D_HRVO',...
                   'agent_2D_RVO2'};
% monteCarlo_algs = {'agent_2D_RVO2'};
monteCarlo_agentN = [2,5,10,20];
monteCarlo_noiseSigma = 0.2;
monteCarlo_outputDir = strcat(pwd,'\data\UKACC_realisticSensing');
monteCarlo_cycles = 1000;   % [FIXED]
monteCarlo_parallel = 1;
monteCarlo_offOnComplete = 0;

% FOR EACH NUMBER OF AGENTS IN THE GIVEN SCENARIO
for n = 1:numel(monteCarlo_agentN)                                         % Differing agent numbers
    % FOR EACH ALGORITHM TO BE EVALUATED
    for index = 1:numel(monteCarlo_algs)                                   % Differing agent algorithms
        % CREATE THE AGENT SET
        agentVector = cell(1,monteCarlo_agentN(n));                        % Agent container
        for i = 1:monteCarlo_agentN(n)
            agentVector{1,i} = eval(monteCarlo_algs{index});               % Assign agents
        end

        % ///////////////////// SCENARIO CONFIGURATION ////////////////////
        initialisedObjects = getScenario_UKACC(...                         % Initialise the objects in the UKACC scenario
            'agents',agentVector,...
            'agentOrbit',sim_agentOrbit,...
            'agentVelocity',sim_agentVelocity,...
            'waypointOrbit',sim_waypointOrbit,...
            'offsetAngle',0);
        clear agentVector

        % UNIQUE OUTPUT FILE
        studyOutputDir = strcat(monteCarlo_outputDir,...
                                '\',monteCarlo_algs{index});

        % ///////////////// DEFINE THE MONTE-CARLO INSTANCE ///////////////
        fprintf('[SETUP]\t Defining Monte-Carlo Simulation Series...\n');
        [monteCarlo] = OMAS_monteCarlo(...
                        'objects',initialisedObjects,...
                        'duration',sim_maxDuration,...
                        'dt',sim_timeStep,...
                        'warningDistance',sim_warningDistance,...
                        'verbosity',sim_verbosity,...
                        'outputPath',studyOutputDir,...
                        'positionSigma',monteCarlo_noiseSigma,...
                        'threadPool',monteCarlo_parallel,...
                        'shutDownOnComplete',monteCarlo_offOnComplete,...
                        'cycles',monteCarlo_cycles);

        % ///////////////// EXECUTE CYCLE EVALUATION //////////////////////
        [summaryData,monteCarlo] = monteCarlo.evaluateCycles();            % Process the defined cycles
        clear summaryData monteCarlo
    end
end

repo: OpenMAS


text:
function [AC,Rank,fS] = UpdateConvergenceArchive(AC,Population,NC,Z,Xic,sigma_niche,Problem)
% Update the Convergence Archive

%------------------------------- Copyright --------------------------------
% Copyright 2017-2018 Yiping Liu
% Please contact {yiping0liu@gmail.com} if you have any problem.
%--------------------------------------------------------------------------

    Population = [Population,AC];
    PopObj = Population.objs;
    PopDec = Population.decs;
    N = size(PopObj,1);
    Rank   = inf(1,N);      % Rank of each solution
    nrank  = 1;             % Current rank

    %% Normalize
    PopObj = PopObj - repmat(Z,N,1);
    PopDec = (PopDec - repmat(Problem.lower,N,1))./repmat(Problem.upper-Problem.lower,N,1);

    %% Convergence indicator
    fS = mean(PopObj');

    %% Calculate distance between every two solutions in the IC decsion subspace
    d = pdist2(PopDec(:,Xic),PopDec(:,Xic),'chebychev');

    %% Rank
    Choose = false(1,N);
    Q = true(1,N);
    Q1 = false(1,N);
    while sum(Choose) < NC
        if sum(Q) == 0
           Q = Q1;
           Q1 = false(1,N);
           nrank = nrank+1;
        end
        % Choose x with min FS
        temp1 = fS == min(fS(Q));
        xmin = find(and(temp1,Q));
        xmin = xmin(1);
        Rank(xmin) = nrank;
        Choose(xmin) = true;
        Q(xmin) = false;
        % Delete solution near x_min
        temp3=d(xmin,:);
        temp2=temp3<sigma_niche;
        Delete = and(temp2,Q);
        Q(Delete) = false;
        Q1(Delete) = true;
    end
    AC = Population(Choose);
    Rank = Rank(Choose);
    fS = fS(Choose);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function ci = cilssvm(model,alpha,conftype)
%
% Construction of bias corrected 100(1-\alpha)% pointwise or
% simultaneous confidence intervals
%
% >> ci = cilssvm({X,Y,type,gam,kernel_par,kernel,preprocess},alpha,conftype)
% >> ci = cilssvm(model,alpha,conftype)
%
% This function calculates bias corrected 100(1-\alpha)% pointwise or
% simultaneous confidence intervals. The procedure support homoscedastic
% data sets as well heteroscedastic data sets. The construction of the
% confidence intervals are based on the central limit theorem for linear
% smoothers combined with bias correction and variance estimation.
%
% 1. Using the functional interface:
%
%
% >> ci = cilssvm({X,Y,type,gam,kernel_par,kernel,preprocess})
% >> ci = cilssvm({X,Y,type,gam,kernel_par,kernel,preprocess}, alpha)
% >> ci = cilssvm({X,Y,type,gam,kernel_par,kernel,preprocess}, alpha, conftype)
%
%
%      Outputs
%        ci            : N x 2 matrix containing the lower and upper confidence intervals
%
%      Inputs
%        X             : N x d matrix with the inputs of the training data
%        Y             : N x 1 vector with the outputs of the training data
%        type          : 'function estimation' ('f') or 'classifier' ('c')
%        gam           : Regularization parameter
%        sig2          : Kernel parameter(s) (bandwidth in the case of the 'RBF_kernel')
%        kernel(*)     : Kernel type (by default 'RBF_kernel')
%        preprocess(*) : 'preprocess'(*) or 'original'
%        alpha(*)      : Significance level (by default 5%)
%        conftype(*)   : Type of confidence interval 'pointwise' or 'simultaneous' (by default 'simultaneous')
%
% 2. Using the object oriented interface:
%
%
% >> ci = cilssvm(model)
% >> ci = cilssvm(model, alpha)
% >> ci = cilssvm(model, alpha, conftype)
%
%
%      Outputs
%        ci          : N x 2 matrix containing the lower and upper confidence intervals
%
%      Inputs
%        model       : Object oriented representation of the LS-SVM model
%        alpha       : Significance level (by default 5%)
%        conftype    : Type of confidence interval 'pointwise' or 'simultaneous' (by default 'simultaneous')
%
%
%  See also:
%    trainlssvm, simlssvm, predlssvm

% Copyright (c) 2011,  KULeuven-ESAT-SCD, License & help @ http://www.esat.kuleuven.be/sista/lssvmlab

if iscell(model)
    model = initlssvm(model{:});
end

if nargin <= 1
    alpha = 0.05;
    conftype = 'simul';

elseif nargin <= 2
    conftype = 'simul';
end

if model.preprocess(1)=='p'
    error('Please use original data to compute confidence intervals...')
end

x = model.xtrain; y = model.ytrain;

% train model
if isempty(model.gam) && isempty(model.kernel.pars)
    error('Please tune model first with ''tunelssvm'' to obtain tuning parameters');
end
model = trainlssvm(model);

s = smootherlssvm(model);
Yhat = simlssvm(model,x);

% bias: double smoothing with fourt order kernel RBF4
modelb = initlssvm(x,y,'f',[],[],'RBF4_kernel','o');
modelb = tunelssvm(modelb,'simplex','crossvalidatelssvm',{10,'mse'});
modelb = trainlssvm(modelb);

biascorr = (s-eye(size(x,1)))*simlssvm(modelb,x);

% construct approximate 100(1-alpha)% confidence interval
%1) estimate variance nonparametrically
sigma2 = varest(model);

%2) calculate var-cov matrix
s = s*diag(sigma2)*s';

%2b) find standardized absolute maxbias
delta = max(abs(biascorr./sqrt(diag(s))));

%3) pointwise or simultaneous?
if conftype(1)=='s'
    z = tbform(model,alpha) + delta;
elseif conftype(1)=='p'
    z = norminv(alpha/2);
    Yhat = Yhat - biascorr;
else
    error('Wrong type of confidence interval. Please choose ''pointwise'' or ''simultaneous''');
end

ci = [Yhat+z*sqrt(diag(s)) Yhat-z*sqrt(diag(s))];

function [var,modele] = varest(model)

% if preprocessed data, construct original data
if model.preprocess(1)=='p'
    [x,y] = postlssvm(model,model.xtrain,model.ytrain);
else
    x = model.xtrain; y = model.ytrain;
end

model = trainlssvm(model);

Yh = simlssvm(model,x);

% Squared normalized residuals
e2 = (y-Yh).^2;

% Make variance model
if model.nb_data <= 200
    costfun = 'leaveoneoutlssvm'; costargs = {'mae'};
else
    costfun = 'crossvalidatelssvm'; costargs = {10,'mae'};
end
modele = initlssvm(x,e2,'f',[],[],'RBF_kernel');
modele = tunelssvm(modele,'simplex',costfun,costargs);
modele = trainlssvm(modele);

% variance model
var = max(simlssvm(modele,x),0);

% make estimate of var unbiased in homoscedastic case if regression
% estimate is unbiased
L = smootherlssvm(model);
S = smootherlssvm(modele);

var = var./(ones(size(x,1),1)+S*diag(L*L'-L-L'));

author: peterhcharlton

repo: RRest


text:
%% Example of using KST class for interfacing with KUKA iiwa robots

% soft real-time control of the KUKA robot with
% impedence
% Moving the joints of the robot using a sinusoidal function

% The external torques are plotted in real-time during the test
% the feedback from the measured torques can be used to for a closed loop
% control

% First start the server on the KUKA iiwa controller
% Then run this script using Matlab

% This code was tested using Matlab 2013b

% This example works with Sunrise application version KST_1.7  and higher.
% Copyright Mohammad SAFEEA, 26th-June-2018

close all;clear;clc;
warning('off')
%% Create the robot object
ip='172.31.1.147'; % The IP of the controller
arg1=KST.LBR7R800; % choose the robot iiwa7R800 or iiwa14R820
arg2=KST.Medien_Flansch_elektrisch; % choose the type of flange
Tef_flange=eye(4); % transofrm matrix of EEF with respect to flange
iiwa=KST(ip,arg1,arg2,Tef_flange); % create the object

%% Start a connection with the server
flag=iiwa.net_establishConnection();
if flag==0
return;
end
pause(1);
disp('Doing some stuff')

%% Move point to point to an initial position
jPos={0,0,0,-pi/2,0,pi/2,0};
relVel=0.15;
iiwa.movePTPJointSpace(jPos, relVel); % move to initial configuration

%% Pause for 3 seocnds
pause(3);

%% Tool/Impedance parameters
massOfTool=0.5; % the mass of the tool attached to flange in Kg
cOMx=0; % X coordinate of the center of mass of the tool in (mm)
cOMy=0; % Y coordinate of the center of mass of the tool in (mm)
cOMz=40; % Z coordinate of the center of mass of the tool in (mm)
cStiness=900; % cartizian stifness
rStifness=80; % rotational stifness
nStifness=50; % null space stifness

% Start the soft realtime control with impedance
iiwa.realTime_startImpedanceJoints(massOfTool,cOMx,cOMy,cOMz,...
cStiness,rStifness,nStifness);

w=0.6; % motion constants, frequency rad/sec
A=0.2; % motion constants, amplitude of motion

a=datevec(now);
t0=a(6)+a(5)*60+a(4)*60*60; % calculate initial time

dt=0;
tstart=t0;
counter=0;
duration=0.5*60; %0.5 minutes

%% Control loop
% real time plot handles
 colors={'k','b','r','g','c','m','y'};
 figureHandle=figure('Units','inches','Position',[0 0 5 3.75]);
 numOfSamples=120; % number of samples to show in the plot
 timeVec=1:numOfSamples;
 tawVec=zeros(7,numOfSamples);
 plotHandle=[];
 for i=1:7
    plotHandle=[plotHandle,plot(timeVec,tawVec(i,:),colors{i},'LineWidth',2)];
    hold on;
 end
% format the plot
ylim([0,8]);
temp=title('External torques');
set(temp,'FontSize',16);
temp=xlabel('Time (seconds)');
set(temp,'FontSize',14);
temp=ylabel('Distance (m)');
set(temp,'FontSize',14);

while(dt<duration)
 %% perform trajectory calculation here
  a=datevec(now);
  time=a(6)+a(5)*60+a(4)*60*60;
  dt=time-t0;

  temp=A*(1-cos(w*dt));
    for dacount=1:7
        jPosCommand{dacount}=jPos{dacount}+temp;
    end
  counter=counter+1;
  %% Send joint positions to robot
  taw=iiwa.sendJointsPositionsExTorque(jPosCommand);

  tawTemp=zeros(7,1);
  for i=1:7
      tawTemp(i)=taw{i};
      torque_array(i,counter)=taw{i};
  end
  % delete first measurment
  tawVec(:,1)=[];
  % add new measurment to end of arrays
  tawVec=[tawVec,tawTemp];
  % update plots data
  for i=1:6
    set(plotHandle(i),'YData',tawVec(i,:));
  end
  set(plotHandle(7),'YData',tawVec(i,:));
  % show plot data
  set(figureHandle,'Visible','on');

end

%% Stop the realtime control with impedence motion
iiwa.realTime_stopImpedanceJoints();
fprintf('\n Motion stopped \n');
pause(2);

%% turn off the server
iiwa.net_turnOffServer();

warning('on')

repo: KST-Kuka-Sunrise-Toolbox


text:
% NOTE: This srcipt will not run as expected unless you fill in proper
% code in trajhandle and controlhandle
% You should not modify any part of this script except for the
% visualization part
%
% ***************** MEAM 620 QUADROTOR SIMULATION *****************
close all
clear all
addpath('utils')
addpath('trajectories')

% You can change trajectory here

% trajectory generator
% trajhandle = @step;
% trajhandle = @circle;
trajhandle = @diamond;

% controller
controlhandle = @controller;

% real-time
real_time = true;

% *********** YOU SHOULDN'T NEED TO CHANGE ANYTHING BELOW **********
% number of quadrotors
nquad = 1;

% max time
time_tol = 25;

% parameters for simulation
params = crazyflie();

%% **************************** FIGURES *****************************
fprintf('Initializing figures...\n')
h_fig = figure;
h_3d = gca;
axis equal
grid on
view(3);
xlabel('x [m]'); ylabel('y [m]'); zlabel('z [m]')
quadcolors = lines(nquad);

set(gcf,'Renderer','OpenGL')

%% *********************** INITIAL CONDITIONS ***********************
fprintf('Setting initial conditions...\n')
max_iter  = 5000;      % max iteration
starttime = 0;         % start of simulation in seconds
tstep     = 0.01;      % this determines the time step at which the solution is given
cstep     = 0.05;      % image capture time interval
nstep     = cstep/tstep;
time      = starttime; % current time
err = []; % runtime errors
for qn = 1:nquad
    % Get start and stop position
    des_start = trajhandle(0, qn);
    des_stop  = trajhandle(inf, qn);
    stop{qn}  = des_stop.pos;
    x0{qn}    = init_state( des_start.pos, 0 );
    xtraj{qn} = zeros(max_iter*nstep, length(x0{qn}));
    ttraj{qn} = zeros(max_iter*nstep, 1);
end

x         = x0;        % state

pos_tol   = 0.01;
vel_tol   = 0.01;

%% ************************* RUN SIMULATION *************************
OUTPUT_TO_VIDEO = 1;
if OUTPUT_TO_VIDEO == 1
    v = VideoWriter('diamond.avi');
    open(v)
end

fprintf('Simulation Running....')
% Main loop
for iter = 1:max_iter
    iter;
    timeint = time:tstep:time+cstep;

    tic;
    % Iterate over each quad
    for qn = 1:nquad
        % Initialize quad plot
        if iter == 1
            QP{qn} = QuadPlot(qn, x0{qn}, 0.1, 0.04, quadcolors(qn,:), max_iter, h_3d);
            desired_state = trajhandle(time, qn);
            QP{qn}.UpdateQuadPlot(x{qn}, [desired_state.pos; desired_state.vel], time);
            h_title = title(sprintf('iteration: %d, time: %4.2f', iter, time));
        end

        % Run simulation
        [tsave, xsave] = ode45(@(t,s) quadEOM(t, s, qn, controlhandle, trajhandle, params), timeint, x{qn});
        x{qn}    = xsave(end, :)';

        % Save to traj
        xtraj{qn}((iter-1)*nstep+1:iter*nstep,:) = xsave(1:end-1,:);
        ttraj{qn}((iter-1)*nstep+1:iter*nstep) = tsave(1:end-1);

        % Update quad plot
        desired_state = trajhandle(time + cstep, qn);
        QP{qn}.UpdateQuadPlot(x{qn}, [desired_state.pos; desired_state.vel], time + cstep);
        set(h_title, 'String', sprintf('iteration: %d, time: %4.2f', iter, time + cstep))
        if OUTPUT_TO_VIDEO == 1
            im = frame2im(getframe(gcf));
            writeVideo(v,im);
        end
    end
    time = time + cstep; % Update simulation time
    t = toc;
    % Check to make sure ode45 is not timing out
    if(t> cstep*50)
        err = 'Ode45 Unstable';
        break;
    end

    % Pause to make real-time
    if real_time && (t < cstep)
        pause(cstep - t);
    end

    % Check termination criteria
    if terminate_check(x, time, stop, pos_tol, vel_tol, time_tol)
        break
    end
end

if OUTPUT_TO_VIDEO == 1
    close(v);
end

%% ************************* POST PROCESSING *************************
% Truncate xtraj and ttraj
for qn = 1:nquad
    xtraj{qn} = xtraj{qn}(1:iter*nstep,:);
    ttraj{qn} = ttraj{qn}(1:iter*nstep);
end

% Plot the saved position and velocity of each robot
for qn = 1:nquad
    % Truncate saved variables
    QP{qn}.TruncateHist();
    % Plot position for each quad
    h_pos{qn} = figure('Name', ['Quad ' num2str(qn) ' : position']);
    plot_state(h_pos{qn}, QP{qn}.state_hist(1:3,:), QP{qn}.time_hist, 'pos', 'vic');
    plot_state(h_pos{qn}, QP{qn}.state_des_hist(1:3,:), QP{qn}.time_hist, 'pos', 'des');
    % Plot velocity for each quad
    h_vel{qn} = figure('Name', ['Quad ' num2str(qn) ' : velocity']);
    plot_state(h_vel{qn}, QP{qn}.state_hist(4:6,:), QP{qn}.time_hist, 'vel', 'vic');
    plot_state(h_vel{qn}, QP{qn}.state_des_hist(4:6,:), QP{qn}.time_hist, 'vel', 'des');
end
if(~isempty(err))
    error(err);
end

fprintf('finished.\n')

author: yrlu

repo: quadrotor

save_path: github-repos/MATLAB/yrlu-quadrotor


text:
% make histogram of average within-cluster distance for all clusters,
% pooled for all fish

clear all; close all; clc;

outputDir = GetOutputDataDir;

%% Init load
hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

%% Load fish
range_fish = 1:18;

M_Count = cell(length(range_fish),1);

%%
% i_set = 1;
for i_fish = range_fish
    %% load fish
    ClusterIDs = [6,1];
    [cIX_load,gIX_load] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs);
    CellXYZ_norm = getappdata(hfig,'CellXYZ_norm');
    absIX = getappdata(hfig,'absIX');
    cIX_abs = absIX(cIX_load);
    M_xyz_norm = CellXYZ_norm(cIX_abs,:);

    %% compute average within-cluster-distance for each cluster
    numClus = length(unique(gIX_load));
    D = zeros(numClus,1);
    for i_clus = 1:numClus
        IX = find(gIX_load == i_clus);
        XYZ_clus = M_xyz_norm(IX,:);
       D(i_clus) = mean(pdist(XYZ_clus,'euclidean'));
    end

    M_Count{i_fish} = D;
end

%% fig4f: Distribution of within-cluster distance
figure('Position',[500,500,150,120]);hold on;
xbins = 0:50:400;

% pool
N = zeros(length(range_fish),length(xbins)-1);
for i_fish = range_fish
    [N(i_fish,:),edges,bin] = histcounts(M_Count{i_fish},xbins);
end

for i = 1:size(N,2)
    meanN = mean(N(:,i));
    semN = std(N(:,i))/sqrt(length(range_fish));
    plot([edges(i),edges(i)],[0,meanN],'color',[1,0.5,0.5],'linewidth',7.5)
    plot([edges(i),edges(i)],[meanN-semN,meanN+semN],'color',[0.2,0.2,0.2],'linewidth',0.5);
end

%
% h = findobj(gca,'Type','patch');
% h.FaceColor = [0.5 0.5 0.5];
% h.EdgeColor = 'w';
xlim([-25,max(xbins)])
% set(gca,'XTick', 0:0.2:1);
xlabel('avr within-clus dist')
ylabel('count')
ylim([0,60])

author: xiuyechen

repo: FishExplorer


text:
function prior = invgammaPriorParamInit(prior, params)

% INVGAMMAPRIORPARAMINIT Inverse gamma prior model's parameter initialisation.
% FORMAT
% DESC initialises the parameters of the inverse gamma prior with some
% default parameters.
% ARG prior : prior structure to be initialised.
% RETURN prior : prior structure with initial values in place.
%
% SEEALSO : priorCreate, gammaPriorParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006
%
% MODIFICATIONS: Andreas C. Damianou, 2013
%
% PRIOR

if nargin < 2
    prior.a = 1e-6;
    prior.b = 1e-6;
else
    prior.a = params(1);
    prior.b = params(2);
end

prior.transforms.index = [1 2];
prior.transforms.type =  optimiDefaultConstraint('positive');
prior.nParams = 2;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [crack_img] = seg2cracks(seg)
%
% [crack_img] = seg2cracks(seg)
%
%  Converts a segmentation image into a crack-coded image.
%

% Bits:
%       |
%       | 1
% 4 ----+---- 2
%       |
%       | 3

% Values:
% tic
UP    = 1;
RIGHT = 2;
DOWN  = 3;
LEFT  = 4;

JUNCTION = 5;

% dx = uint8(seg ~= image_right(seg));
% dy = uint8(seg ~= image_down(seg) );
% crack_img = dx + bitshift(dy,LEFT-1) + ...
%     bitshift(image_right(dy),RIGHT-1) + bitshift(image_down(dx),DOWN-1);

% Removed dependency on image_down and image_right, to make this more
% easily packaged:
dx = uint8(seg ~= seg(:,[2:end end]));
dy = uint8(seg ~= seg([2:end end],:));
crack_img = dx + bitshift(dy,LEFT-1) + ...
    bitshift(dy(:,[2:end end]),RIGHT-1) + bitshift(dx([2:end end],:),DOWN-1);

% Find interior junctions:
junction_map = (crack_img==11 | crack_img==7 | crack_img==14 | ...
    crack_img==13 | crack_img==15);

% Find the junctions along the borders:
junction_map([1 end],:) = junction_map([1 end],:) | bitget(crack_img([1 end],:),UP);
junction_map(:,[1 end]) = junction_map(:,[1 end]) | bitget(crack_img(:,[1 end]),LEFT);

% set the junction bit for all these junctions in the crack_img
crack_img(junction_map) = bitset(crack_img(junction_map), JUNCTION);

% % Also set all borders of the image to be seen as junctions (all bits set
% % == a value of 15) so that we will know to stop in the fragment chaining
% % process later
% crack_img(:,[1 end]) = 15;
% crack_img([1 end],:) = 15;

% fprintf('New method: %.3f seconds\n', toc);
%
% %% Old Method %%
% tic
% [nrows, ncols ] = size(seg);
% crack_img2 = zeros(nrows, ncols, 'uint8');
% % Bits:
% UP    = 1;
% RIGHT = 2;
% DOWN  = 3;
% LEFT  = 4;
%
% index = find(seg ~= image_right(seg));
% crack_img2(index) = bitset(crack_img2(index), UP);
%
% index = find(image_right(seg) ~= image_downright(seg));
% crack_img2(index) = bitset(crack_img2(index), RIGHT);
%
% index = find(image_down(seg) ~= image_downright(seg));
% crack_img2(index) = bitset(crack_img2(index), DOWN);
%
% index = find(seg ~= image_down(seg));
% crack_img2(index) = bitset(crack_img2(index), LEFT);
%
% fprintf('Old method: %.3f seconds\n', toc);
%
% if(all(crack_img(:)==crack_img2(:)))
%     disp('Results agree.')
% end
%

author: Cloud-CV

repo: object-proposals


text:
function [mean_norm, stdd] = spp_feature_stats(imdb, roidb, layer, spp_model)
% [mean_norm, stdd] = spp_feature_stats(imdb, roidb, layer, spp_model)
%
% Adapted from spp code written by Ross Girshick
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

imdbs_name = cell2mat(cellfun(@(x) x.name, imdb, 'UniformOutput', false));
conf = spp_config('sub_dir', fullfile(spp_model.cache_name, imdbs_name));
save_file = sprintf('%s/feature_stats_layer_%d.mat', ...
                    conf.cache_dir, layer);

t_start = tic();
try
  ld = load(save_file);
  mean_norm = ld.mean_norm;
  stdd = ld.stdd;
  clear ld;
catch
  % fix the random seed for repeatability
  prev_rng = seed_rand();

  image_idx_in_imdb = cell2mat(cellfun(@(x) 1:length(x.image_ids), imdb, 'UniformOutput', false));
  image_imdb_id = cell2mat(cellfun(@(x, y) ones(1, length(x.image_ids))*y, imdb, num2cell(1:length(imdb)), 'UniformOutput', false));

  num_images = min(length(image_idx_in_imdb), 200);
  boxes_per_image = 200;

  valid_idx = randperm(length(image_idx_in_imdb), num_images);

  ns = [];
  for i = 1:length(valid_idx)
    image_idx = valid_idx(i);
    tic_toc_print('feature stats: %d/%d\n', i, length(valid_idx));

    imdb_idx = image_imdb_id(image_idx);
    d = roidb{imdb_idx}.rois(image_idx_in_imdb(image_idx));
    d.feat = spp_load_cached_poolX_features(spp_model.spp_pooler, spp_model.feat_cache{imdb_idx}, ...
        imdb{imdb_idx}.name, imdb{imdb_idx}.image_ids{image_idx_in_imdb(image_idx)}, d.boxes);
    if isempty(d.feat)
        continue;
    end
    X = d.feat(:, randperm(size(d.feat,2), min(boxes_per_image, size(d.feat,2))));
    X = spp_poolX_to_fcX(X, layer, spp_model, conf.use_gpu);

    ns = cat(2, ns, sqrt(sum(X.^2, 1)));
  end

  mean_norm = mean(ns);
  stdd = std(ns);
  save(save_file, 'mean_norm', 'stdd');

  % restore previous rng
  rng(prev_rng);
end
fprintf('spp_feature_stats_spm in %f seconds.\n', toc(t_start));

author: ShaoqingRen


text:
% SetValueOfAssignment Sets the value of a variable assignment in a factor.
%
%   F = SetValueOfAssignment(F, A, v) sets the value of a variable assignment,
%   A, in factor F to v. The order of the variables in A are assumed to be the
%   same as the order in F.var.
%
%   F = SetValueOfAssignment(F, A, v, VO) sets the value of a variable
%   assignment, A, in factor F to v. The order of the variables in A are given
%   by the vector VO.
%
%   Note that SetValueOfAssignment *does not modify* the factor F that is
%   passed into the function, but instead returns a modified factor with the
%   new value(s) for the specified assignment(s). This is why we have to
%   reassign F to the result of SetValueOfAssignment in the code snippets
%   shown above.
%
%   See also GetValueOfAssignment.m and SampleFactors.m
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function F = SetValueOfAssignment(F, A, v, VO)

if (nargin == 3),
    indx = AssignmentToIndex(A, F.card);
else
    map = zeros(length(F.var), 1);
    for i = 1:length(F.var),
        map(i) = find(VO == F.var(i));
    end;
    indx = AssignmentToIndex(A(map), F.card);
end;

F.val(indx) = v;

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
% This program obtains the analytical step response of a parallel RLC
% circuit. In addition, graphs of the voltage across the parallel
% combination and branch currents are obtained.
% H. Saadat, Copyright 1999

function PresponseIL
warning('off','MATLAB:dispatcher:InexactMatch')
Rhandle=findobj('Tag','Rtext');
R=eval(get(Rhandle,'String'));
KOhandle=findobj('Tag','KOhm');
KO=get(KOhandle,'Value');
if KO==1
   R=1000*R;
else, end
Lhandle=findobj('Tag','Ltext');
L=eval(get(Lhandle,'String'));
mHhandle=findobj('Tag','mHenry');
ml=get(mHhandle,'Value');
if ml==1
   L=L/1000;
else, end
Chandle=findobj('Tag','Ctext');
C=eval(get(Chandle,'String'));
micFhandle=findobj('Tag','micFarad');
micF=get(micFhandle,'Value');
if micF==1
   C=C/1000000;
else, end
Rhandle=findobj('Tag','IStext');
Is=eval(get(Rhandle,'String'));
Rhandle=findobj('Tag','I0text');
I0=eval(get(Rhandle,'String'));
Rhandle=findobj('Tag','V0text');
V0=eval(get(Rhandle,'String'));
% Step Response of a parallel RLC circuit
if L==inf | C==inf
   plot(0, 0)
   text(-0.8, 0.80, 'Parameters must lie in the following ranges:', 'color', [0 0 0.6275])
   text(-0.8, 0.65, '0 < R \leq inf,   0 < L < inf,   &    0 < C < inf', 'color', [0 0 0.6275])
   text(-0.8, 0.5, 'Enter the correct value and press Solve', 'color', [0 0 0.6275])
   axis([-1  1  -1  1]);return,  else
end
if R==0 | L==0 | C == 0
plot(0, 0)
title('Short-circuit across the parallel circuit', 'color', [ 1 0 0])
text(-0.8, 0.80, 'Parameters must lie in the following ranges:', 'color', [0 0 0.6275])
text(-0.8, 0.65, '0 < R \leq inf,   0 < L < inf,   &    0 < C < inf', 'color', [0 0 0.6275])
text(-0.8, 0.5, 'Enter the correct value and press Solve', 'color', [0 0 0.6275])
axis([-1  1  -1  1]);return,  else
alpha =1/(2*R*C);
end
w02 = 1/(L*C); w0 = sqrt(w02);
if alpha~= 0
   err=(alpha-w0)/alpha;
   if abs(err) <1e-4
   w0=alpha;
   else, end
else, end
dv = Is/C-I0/C-V0/(R*C);
di = V0/L;
if alpha > w0
 % Overdamped response
    s(1) = -alpha + sqrt(alpha^2 - w02);
    s(2) = -alpha - sqrt(alpha^2 - w02);
    A1 = (dv-s(2)*V0)/(s(1)-s(2));
    A2 = (dv-s(1)*V0)/(s(2)-s(1));
    A1i =(di-s(2)*(I0-Is))/(s(1)-s(2));
    A2i =(di-s(1)*(I0-Is))/(s(2)-s(1));
    tf = 6*max(abs(1/s(1)), abs(1/s(2)));
    t=0:tf/100:tf;
    v = A1*exp(s(1)*t)+A2*exp(s(2)*t);
    iL=Is + A1i*exp(s(1)*t)+A2i*exp(s(2)*t);
    iR = v/R;
    iC= Is -iR-iL;
    it=iR+iL+iC;
    plot(t,iL,'erasemode','none','color',[0.9 0 0.8]), grid
    if Is==0
    title(['I_L(t) = (', num2str(A1i), ') e^{(', num2str(s(1)), 't)} + (', num2str(A2i), ') e^{(', num2str(s(2)), 't)}'],'color',[0.9 0 0.8])
    elseif Is~=0
    title(['i_L(t) = ',num2str(Is), ' + (', num2str(A1i), ') e^{(', num2str(s(1)), 't)} + (', num2str(A2i), ') e^{(', num2str(s(2)), 't)}'],'color',[0.9 0 0.8])
    end
    xlabel(['\alpha = ',num2str(alpha), ',  \omega_0 = ',  num2str(w0), '  (O.D.)          t, sec'])
    ylabel('i_L(t), Amps')
    elseif alpha < w0
     if alpha~=0
      % Underdamped response
      tf = 6/alpha;
      t=0:tf/200:tf;
      wd= sqrt(w02 - alpha^2);
      s(1) = -alpha +j*wd;
      s(2) = -alpha -j*wd;
      B1 = V0; B2 = (dv+alpha*B1)/wd;
      B1i = I0 - Is; B2i = (di+alpha*B1i)/wd;
      v = exp(-alpha*t).*(B1*cos(wd*t)+B2*sin(wd*t));
      iL = Is + exp(-alpha*t).*(B1i*cos(wd*t)+B2i*sin(wd*t));
      iR = v/R;
      iC= Is -iR-iL;
      it=iR+iL+iC;
      plot(t,iL,'erasemode','none','color',[.9 0 0.8]), grid
         if Is == 0
         title(['I_L(t) = e^{(-', num2str(alpha), 't)}[(', num2str(B1i), ')cos',num2str(wd),'t + (' , num2str(B2i), ')sin',num2str(wd),'t]'],'color',[0.9 0 0.8])
         elseif Is~= 0
         title(['i_L(t) = ',num2str(Is), ' + e^{(-', num2str(alpha), 't)}[(', num2str(B1i), ')cos',num2str(wd),'t + (' , num2str(B2i), ')sin',num2str(wd),'t]'],'color',[0.9 0 0.8])
         end
      xlabel(['\alpha = ',num2str(alpha), ',  \omega_0 = ',  num2str(w0), '  (U.D.)          t, sec'])
      ylabel('i_L(t), Amps')
      else
      % Undamped Response
      tf = 8*pi/w0;
      t=0:tf/200:tf;
      wd= w0;
      B1 = V0; B2 = (dv)/wd;
      B1i = I0 - Is; B2i = (di)/wd;
      iL = Is + B1i*cos(wd*t)+B2i*sin(wd*t);
      v = B1*cos(wd*t)+B2*sin(wd*t);
      iR = v/R;
      iC= Is -iR-iL;
      it=iR+iL+iC;
      plot(t,iL,'erasemode','none','color',[0.9 0 0.80]), grid
      	if Is == 0
      	title(['I_L(t) = (', num2str(B1i), ')cos',num2str(wd),'t + (' , num2str(B2i), ')sin',num2str(wd),'t'],'color',[0.9 0 0.8])
      	elseif Is ~= 0
      	title(['i_L(t) = ',num2str(Is), ' + (', num2str(B1i), ')cos',num2str(wd),'t + (' , num2str(B2i), ')sin',num2str(wd),'t'],'color',[0.9 0 0.8])
      	end
      xlabel(['\alpha = ',num2str(alpha), ',  \omega_0 = ',  num2str(w0), '  (Undamped)          t, sec'])
      ylabel('i_L(t), Amps')
      end
   elseif alpha ==w0
    % Critically damped response
    s(1)=-alpha;  s(2) = s(1);
    D2 = V0;
    D1 = dv+alpha*D2;
    D2i = I0 - Is; D1i = di+alpha*D2i;
    tf = 8/alpha;
    t=0:tf/100:tf;
    iL=Is + exp(-alpha*t).*(D1i*t+D2i);
    v = exp(-alpha*t).*(D1*t+D2);
    iR = v/R;
    iC= Is -iR-iL;
    it=iR+iL+iC;
    plot(t,iL,'erasemode','none','color',[0.9 0 0.8]), grid
    	if Is == 0
          title(['i_L(t) = e^{(-', num2str(alpha), 't)}[(', num2str(D1i), ')t + (' , num2str(D2i), ')]'],'color',[0.9 0 0.8])
          elseif Is ~= 0
          title(['i_L(t) = ',num2str(Is), ' + e^{(-', num2str(alpha), 't)}[(', num2str(D1i), ')t + (' , num2str(D2i), ')]'],'color',[0.9 0 0.8])
      end
    xlabel(['\alpha = ',num2str(alpha), ',  \omega_0 = ',  num2str(w0), '  (C.D.)          t, sec'])
    ylabel('i_L(t), Amps')
    else, end

author: Sable

repo: mcbench-benchmarks