text:
%Coded by Dr. T. Gu (gu@ohio.edu). Copyright 2009.
%A Pro version is available for up to 99 equations.
%User can import equations and export results
function varargout = MfigExtract(varargin)
% MFIGEXTRACT M-file for MfigExtract.fig
%      MFIGEXTRACT, by itself, creates a new MFIGEXTRACT or raises the existing
%      singleton*.
%
%      H = MFIGEXTRACT returns the handle to a new MFIGEXTRACT or the handle to
%      the existing singleton*.
%
%      MFIGEXTRACT('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MFIGEXTRACT.M with the given input arguments.
%
%      MFIGEXTRACT('Property','Value',...) creates a new MFIGEXTRACT or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before MfigExtract_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to MfigExtract_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

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

% Last Modified by GUIDE v2.5 23-Apr-2009 13:49:27

% Begin initialization code - DO NOT EDIT

thisver=version; %prepare to check version = R14 or above?
if isempty(strfind(thisver,'R13'))~=1 | isempty(strfind(thisver,'R12'))~=1
    disp('This software requires MATLAB R14 or above! Yours is lower. Quit.')
    return
end

gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @MfigExtract_OpeningFcn, ...
                   'gui_OutputFcn',  @MfigExtract_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

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

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

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

% Update handles structure
guidata(hObject, handles);

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

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

% Get default command line output from handles structure
varargout{1} = handles.output;

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

%Let's close figure windows before opening new figure to extract
figs = findobj('Type','figure');
figs = sort(figs);
if length(figs) > 1  %if more than just GUI main
   for n = 1:(length(figs)-1) %skip last figure (GUI main)
       close(figs(n)); pause(0.0000001); %delay needed
   end
end

extractxydata_Callback(hObject, eventdata, handles)

function extractxydata_Callback(hObject, eventdata, handles)
% hObject    handle to extractxydata (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global pathname2 filename2  %for View Excel file button

[filename1 pathname1]=uigetfile('*.fig','Select a figure file to extract x-y data');
if isequal(filename1,0) | isequal(pathname1,0); return; end %User pressed Cancel

[filename2, pathname2] = uiputfile('*.xls', 'Provide an Excel file name to save the x-y data');
if isequal(filename2,0) | isequal(pathname2,0)
    return
end %User pressed Cancel

if isempty(strfind(filename2,'.xls'));
    filename2=[filename2 '.xls'];
end

s = openfig(fullfile(pathname1,filename1)); % open figure and get handle to figure
delete(legend);  %legends add extra data to output. Remove legends

%get line handles
h = findobj(s,'Type','line');   %line is the type of your figure file.
	if isempty(h)
	line{1}='Your figure file does not contain x-y data.';
    line{2}='Did you accidentally load a wrong file?';
	line{3}='Your can view your figures using the File menu.';
	line{4}='Make sure you do not load a wrong file.';
	uiwait(msgbox(line, 'No x-y data in your file','warn'));
	close(s)
    return
	end

x=get(h,'xdata');
y=get(h,'ydata');
close(s) %close figure with legends removed
s=hgload(fullfile(pathname1,filename1));  %show original figure with legends

[n dummy]=size(x);  %n= number of x-y sets in the figure.

fid=fopen(fullfile(pathname2,filename2), 'wt');
fprintf(fid,'%s %s\n', 'x, y data extracted from ', [pathname1, filename1]);
if n==1   %because x{i} becomes invalid symbol below if n=1
fprintf(fid,'%s\t %s\n', '      x1', '      y1');
fprintf(fid,'%g\t %g\n', [x;y]);
elseif n > 1

%x{1}=x1 series, y{1} is y1 series, etc.
%find out maximum x{i} length
maxsize=length(x{1});
for i=2:n
    if length(x{i}) > maxsize
        maxsize=length(x{i});
    end
end

%find out whether all the x series are the same
xseriessame=1;
for i=2:n
    if isequal(x{1},x{i})==0  %0 means not equal
        xseriessame=0;
    end
end

if xseriessame==1;  %write x y1 y2... header
fprintf(fid,'All the x series in the figure are found to be the same.\n');
fprintf(fid,'       x\t');
    for i=1:n
    fprintf(fid,'       y%g\t',i);
    end
fprintf(fid,'\n'); %carriage return for header line

%Write each x,y pair horizontally in the output file
    for j=1:maxsize
        fprintf(fid,'%g\t',x{1}(j));  %x series are the same! Write x{1} for all
        for i=1:n
            fprintf(fid,'%g\t',y{i}(j));
        end
    fprintf(fid,'\n');  %carriage return at end of each row
    end
elseif xseriessame==0

%write x, y column header
for i=1:n
       fprintf(fid,'       x%g\t        y%g\t',i,i);
end
fprintf(fid,'\n'); %carriage return for header line

%Write each x,y pair horizontally in the output file
for j=1:maxsize
        for i=1:n
            if j > length(x{i})   %the x-y pair has ended earlier than others
                fprintf(fid,'\t \t');  %mark current x, y as empty
            elseif j <= length(x{i})
                fprintf(fid,'%g\t %g\t',x{i}(j),y{i}(j));
            end
        end
fprintf(fid,'\n');  %carriage return at end of each row
end

end   %ending elseif xseriessame==0 statement
end   %ending elseif n>1 statement

fclose(fid);
message{1}='Output from last run:';
message{3}=sprintf('File = %s',filename2);
message{5}=sprintf('Directory = %s',pathname2);
set(handles.text4,'String',message);
set(handles.pushbutton2,'Visible','on')

% --------------------------------------------------------------------
function viewfigure_Callback(hObject, eventdata, handles)
% hObject    handle to viewfigure (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename pathname]=uigetfile('*.fig','Select a figure file to extract x-y data');
if isequal(filename,0) | isequal(pathname,0); return; end %User pressed Cancel
s=hgload(fullfile(pathname,filename));

% --------------------------------------------------------------------
function viewExcelfile_Callback(hObject, eventdata, handles)
% hObject    handle to viewExcelfile (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename pathname]=uigetfile('*.xls; *.xlsx','Select a figure file to extract x-y data');
if isequal(filename,0) | isequal(pathname,0); return; end %User pressed Cancel
winopen(fullfile(pathname,filename));

% --------------------------------------------------------------------
function about_Callback(hObject, eventdata, handles)
% hObject    handle to about (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
message{1}='This freeware was coded by:';
message{2}='Dr. Tingyue Gu ( gu@ohio.edu )';
message{3}='Department of Chemical and Biomolecular Eng.';
message{4}='Ohio University, Athens, Ohio 45701, USA';
msgbox(message,'About MfigExtract','none');

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global pathname2 filename2
winopen(fullfile(pathname2,filename2))


text:
function [ output_args ] = setConfig( config )
%setConfig Summary of this function goes here
%   Detailed explanation goes here

end

author: Cloud-CV

repo: object-proposals


text:
function [ goodClusters ] = ClusterSelection( clusters, image, s, w )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
height = length(image(:,1,1));
width = length(image(1,:,1));
discHeight = height/s;
discWidth = width/s;
numClusters = size(clusters,3);
windowClusters = zeros(discHeight,discWidth);
dists = zeros(1,numClusters);
boxedImage = OutlineRegion(image,ones(discHeight,discWidth));

for i=1:discHeight
    for j=1:discWidth
        hist = SimpleHist1D(image((i-1)*s+1:i*s,(j-1)*s+1:j*s,:),w);
        for k=1:numClusters
            dists(k) = Distance1D(hist,clusters(:,:,k));
        end
        [dist,windowClusters(i,j)] = min(dists); %#ok<ASGLU>
    end
end

display(windowClusters);
goodClusters = [];
response = '';

binaryImage = zeros(discHeight,discWidth);
for i=1:length(goodClusters)
    binaryImage = binaryImage + (windowClusters == goodClusters(i));
end
currentImage = DeleteWindowsImage(boxedImage,binaryImage);
imshow(currentImage);

while (not(strcmp(response,'done')))
    [x,y] = ginput(1);
    i = floor(y/s)+1;
    j = floor(x/s)+1;
    if (any(windowClusters(i,j)==goodClusters))
        goodClusters = goodClusters(goodClusters ~= windowClusters(i,j));
    else
        goodClusters = [goodClusters windowClusters(i,j)]; %#ok<AGROW>
        display(goodClusters);
    end
    binaryImage = zeros(discHeight,discWidth);
    for i=1:length(goodClusters)
        binaryImage = binaryImage + (windowClusters == goodClusters(i));
    end
    currentImage = DeleteWindowsImage(boxedImage,binaryImage);
    imshow(currentImage);
    response = input('more, or done?', 's');
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [lb,ub,A,b,infeasible] = remove_bounds_from_Aeqbeq(A,b,lbin,ubin);

infeasible = 0;
lb = lbin;
ub = ubin;
if size(A,1)>0
    cand_rows_lp = find(sum(A~=0,2)==1);
    if ~isempty(cand_rows_lp)
        [ii,jj,kk] = find(A(cand_rows_lp,:));
        s_pos = find(kk>0);
        s_neg = find(kk<=0);
        if ~isempty(s_pos)
            for s = 1:length(s_pos)
                lb(jj(s_pos(s)),1) = full(b(cand_rows_lp(ii(s_pos(s))))./kk(s_pos(s)));
                ub(jj(s_pos(s)),1) = full(b(cand_rows_lp(ii(s_pos(s))))./kk(s_pos(s)));
            end
        end
        if ~isempty(s_neg)
            for s = 1:length(s_neg)
                lb(jj(s_neg(s)),1) = full(b(cand_rows_lp(ii(s_neg(s))))./kk(s_neg(s)));
                ub(jj(s_neg(s)),1) = full(b(cand_rows_lp(ii(s_neg(s))))./kk(s_neg(s)));
            end
        end
    end
    A(cand_rows_lp,:) = [];
    b(cand_rows_lp,:) = [];
end
if any(lb > ubin)
    j = find(lb > ubin);
    ub(j) = ubin(j);
    infeasible = 1;
end
if any(ub < lbin)
    j = find(ub < lbin);
    lb(j) = lbin(j);
    infeasible = 1;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function H=haarmtx(n)

% HAARMTX Compute Haar orthogonal transform matrix.
%
%   H = HAARMTX(N) returns the N-by-N HAAR transform matrix.  H*A
%   is the HAAR transformation of the columns of A and H'*A is the inverse
%   transformation of the columns of A (when A is N-by-N).
%   If A is square, the two-dimensional Haar transformation of A can be computed
%   as H*A*H'. This computation is sometimes faster than using
%   DCT2, especially if you are computing large number of small
%   Haar transformation, because H needs to be determined only once.
%
%   Class Support
%   -------------
%   N is an integer scalar of class double. H is returned
%   as a matrix of class double.
%
%
%   I/O Spec
%   N - input must be double
%   D - output DCT transform matrix is double
%
%   Author : Frederic Chanal (f.j.chanal@student.tue.nl) - 2004

a=1/sqrt(n);
for i=1:n
    H(1,i)=a;
end
for k=1:n-1
    p=fix(log2(k));
    q=k-2^p+1;
    t1=n/2^p;
    sup=fix(q*t1);
    mid=fix(sup-t1/2);
    inft=fix(sup-t1);
    t2=2^(p/2)*a;
    for j=1:inft
        H(k+1,j)=0;
    end
    for j=inft+1:mid
        H(k+1,j)=t2;
    end
    for j=mid+1:sup
        H(k+1,j)=-t2;
    end
    for j=sup+1:n
        H(k+1,j)=0;
    end
end

author: Sable

repo: mcbench-benchmarks


text:
% The COBRAToolbox: testMgPipe.m
%
% Purpose:
%     - test the mappingInfo function of the Microbiome Toolbox
%
% Authors:
%     - original version: Federico Baldini & Laurent Heirendt, April 2018
%       - Almut Heinken, December 2020: adapted to recent changes in mgPipe
%

global CBTDIR

% define the features required to run the test
requiredToolboxes = { 'distrib_computing_toolbox' };

% test Requirements
prepareTest('requiredToolboxes',requiredToolboxes, 'requiredSolvers', {'ibm_cplex'});

% save the current path
currentDir = pwd;

% initialize the test
cd(fileparts(which(mfilename)));
modPath = [CBTDIR filesep 'test' filesep 'models' filesep 'mat'];

% path to the file with abundance information.
abunFilePath='testData_normCoverageReduced.csv';

% path to the file with dietary information
dietFilePath='AverageEuropeanDiet';

% path to the file with information on individuals (default=none)
indInfoFilePath='';

% to define whether uptake and secretion profiles should be recrutied
computeProfiles = true;

resPath = [pwd filesep 'Results'];

% number of cores dedicated for parallelization
numWorkers = 2;

% name of objective function of organisms
objre={'EX_biomass(e)'};

% if documentations on patient health status is provided (0 not 1 yes)
patStat = 0;

% to enable also rich diet simulations (default=false)
rDiet = 0;

% to enable personalized diet simulayions (default=false)
pDiet = 0;

% set the lower bound on biomass in community models (default=0)
lowerBMBound=0;

% set if diet should be adapted with human-derived metabolites
% (default=true)
adaptDiet=1;

% set if existing results should be overwritten (default=false)
repeatSim=1;

% turn all warnings off
warning('off', 'all')

% test getIndividualSizeName
[sampNames, organisms, exMets] = getIndividualSizeName(abunFilePath,modPath);

assert(length(sampNames) == 4)
assert(length(organisms) == 5)

%% run the complete pipeline
% path to the file with abundance information.
abunFilePath='testData_normCoverageReduced.csv';

% test complete pipeline run
[init, netSecretionFluxes, netUptakeFluxes, Y] = initMgPipe(modPath, abunFilePath, computeProfiles, 'dietFilePath', dietFilePath, 'numWorkers', numWorkers, 'resPath', resPath);

assert(init)

% test detectOutput
mapP = detectOutput(resPath, 'mapInfo.mat');
assert(~isempty(mapP));

%% test each function separately

load(['Results' filesep 'mapInfo.mat']);

% load the models
[sampNames, organisms] = getIndividualSizeName(abunFilePath,modPath);
models = loadUncModels(modPath, organisms, objre);

assert(length(models) == 5);

assert(verifyCobraFunctionError('loadUncModels', 'inputs', {}));
assert(verifyCobraFunctionError('loadUncModels', 'inputs', {modPath, ''}));

warning('off', 'all');
loadUncModels(modPath, organisms);
assert(~isempty(lastwarn()))
warning('on', 'all');

% test if the function throws an error when no arguments are provided
assert(verifyCobraFunctionError('initMgPipe'))

% test with only the path to the models (throws an error that the abundance file is missing)
assert(verifyCobraFunctionError('initMgPipe', 'inputs',{modPath}));

% cleanup
delete simRes.mat
delete rDiet_allFlux.csv

% testing with rich diet
[init, netSecretionFluxes, netUptakeFluxes, Y] = initMgPipe(modPath, abunFilePath, computeProfiles, 'resPath',resPath,'numWorkers',numWorkers, 'rDiet', true);
assert(init)

% cleanup
delete simRes.mat
clear ID fvaCt nsCt presol inFesMat
cd(resPath)
delete *.mat

%% verify that mgPipe also works with a single microbiome model with one strain

abunFilePath='testData_normCoverage_one_model.csv';

[init, netSecretionFluxes, netUptakeFluxes, Y] = initMgPipe(modPath, abunFilePath, computeProfiles, 'resPath',resPath,'numWorkers',numWorkers,'resPath',resPath);

assert(init)

% cleanup
rmdir(resPath,'s')

author: opencobra

repo: cobratoolbox


text:
function updateStackedplot(obj, plotIndex)

    %-------------------------------------------------------------------------%

    %-INITIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);
    lineData = plotData.LineProperties(end:-1:1);

    %-get trace data-%

    sourceTable = plotData.SourceTable;
    displayVariables = plotData.DisplayVariables;
    nTraces = length(plotData.AxesProperties);

    if isempty(sourceTable)
        xData = plotData.XData;
        for t = 1:nTraces
            n = nTraces - t + 1;
            yData{t} = plotData.YData(:, n);
        end

    else
        if istimetable(sourceTable)
            xData = date2NumData(sourceTable.Properties.RowTimes);
        else
            xData = 1:size(sourceTable, 1);
        end

        for t = 1:nTraces
            n = nTraces - t + 1;
            yData{t} = date2NumData(sourceTable.(displayVariables{n}));
        end
    end

    %-------------------------------------------------------------------------%

    %-UPDATE STACKEDPLOT AXIS-%
    updateStackedplotAxis(obj, plotIndex)

    %-------------------------------------------------------------------------%

    %-SET TRACES-%
    traceIndex = plotIndex;

    for t = 1:nTraces

        %-update current trace Index-%
        if t ~= 1
            obj.PlotOptions.nPlots = obj.PlotOptions.nPlots + 1;
            traceIndex = obj.PlotOptions.nPlots;
        end

        %-set current trace-%
        obj.data{traceIndex}.type = 'scatter';
        obj.data{traceIndex}.visible = strcmp(plotData.Visible,'on');
        obj.data{traceIndex}.name = plotData.DisplayLabels{t};
        obj.data{traceIndex}.xaxis = sprintf('x%d', 1);
        obj.data{traceIndex}.yaxis = sprintf('y%d', t);

        %-set current trace data-%
        obj.data{traceIndex}.x = xData;
        obj.data{traceIndex}.y = yData{t};

        %-set current trace marker-%
        switch lineData(t).PlotType
            case 'scatter'
                markerSize = ones(size(xData)) * lineData(t).MarkerSize;
                obj.data{traceIndex}.mode = 'markers';
                obj.data{traceIndex}.marker = extractLineMarker(lineData(t));
                obj.data{traceIndex}.marker.size = markerSize;
                obj.data{traceIndex}.marker.line.width = 1;
            otherwise
                obj.data{traceIndex}.mode = 'lines';
                obj.data{traceIndex}.line = extractLineLine(lineData(t));
        end

    end

    %-------------------------------------------------------------------------%
end

function updateStackedplotAxis(obj, plotIndex)

    %-------------------------------------------------------------------------%

    %-INITIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);

    %-------------------------------------------------------------------------%

    %-SET X-AXIS-%

    [xaxis, xExpoFormat] = getAxis(obj, plotIndex, 'X');
    obj.layout = setfield(obj.layout, 'xaxis1', xaxis{1});

    %-------------------------------------------------------------------------%

    %-SET Y-AXIS-%

    [yaxis, yExpoFormat] = getAxis(obj, plotIndex, 'Y');

    for a = 1:length(yaxis)
        obj.layout = setfield(obj.layout, sprintf('yaxis%d', a), yaxis{a});
    end

    %-------------------------------------------------------------------------%

    %-SET AXES ANOTATIONS-%

    %-trace title-%
    updateTitle(obj, plotData.Title, [1, 3]);

    %-exponent for x-axis-%
    updateExponentFormat(obj, xExpoFormat(1), [1,1], 'X');

    %-exponent for y-axis-%
    for a = 1:length(yExpoFormat)
        updateExponentFormat(obj, yExpoFormat(a), [1, a], 'Y');
    end
end

function [ax, expoFormat] = getAxis(obj, plotIndex, axName)

    %-------------------------------------------------------------------------%

    %-INITIALIZATIONS-%

    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);
    plotData = get(obj.State.Plot(plotIndex).Handle);

    lineFactor = obj.PlotlyDefaults.AxisLineIncreaseFactor;
    axisPos = plotData.Position;
    axisColor = getStringColor(zeros(1,3));
    lineWidth = 1;
    fontSize = plotData.FontSize;
    fontFamily = matlab2plotlyfont(plotData.FontName);;
    tickLen = 5;

    %-------------------------------------------------------------------------%

    %-Parse parameters accorging to axisName (X or Y)

    switch axName
        case {'x', 'X'}
            nAxis = 1;
            nTicks = [5, 12];
            axisLim{nAxis} = plotData.XLimits;
            axisLabel{nAxis} = plotData.XLabel;
            axisDomain{nAxis} = min([axisPos(1) sum(axisPos([1,3]))], 1);
            axisAnchor{nAxis} = 'y1';

        case {'x', 'Y'}
            nAxis = length(plotData.AxesProperties);
            yPos = linspace(axisPos(2), sum(axisPos([2,4])), nAxis+1);
            yOffset = diff(yPos)*0.1; yOffset(1) = 0;

            for a = 1:nAxis
                b = nAxis-a+1;
                axisLim{a} = plotData.AxesProperties(b).YLimits;
                axisLabel{a} = plotData.DisplayLabels{b};
                axisDomain{a} = min([yPos(a)+yOffset(a) yPos(a+1)], 1);
                axisAnchor{a} = 'x1';
            end

            if nAxis < 4
                nTicks = [5, 8];
            else
                nTicks = [3, 5];
            end
    end

    %-------------------------------------------------------------------------%

    %-GET EACH AXIS-%

    for a = 1:nAxis

        %-general-%
        ax{a}.domain = axisDomain{a};
        ax{a}.anchor = axisAnchor{a};
        ax{a}.range = date2NumData(axisLim{a});

        ax{a}.side = 'left';
        ax{a}.mirror = false;
        ax{a}.zeroline = false;

        ax{a}.linecolor = axisColor;
        ax{a}.linewidth = lineWidth;
        ax{a}.exponentformat = obj.PlotlyDefaults.ExponentFormat;

        %-grid-%
        if strcmp(plotData.GridVisible, 'on')
            ax{a}.showgrid = true;
            ax{a}.gridwidth = lineWidth;
            ax{a}.gridcolor = getStringColor(0.15*ones(1,3), 0.15);
        else
            ax{a}.showgrid = false;
        end

        %-ticks-%
        if isnumeric(axisLim{a})
            [tickVals, tickText, expoFormat(a)] = getNumTicks(axisLim{a}, nTicks);

            [tickVals, tickText] = rmTicks(ax{a}.range, tickVals, tickText, 's');
            [tickVals, tickText] = rmTicks(ax{a}.range, tickVals, tickText, 'e');

        elseif isduration(axisLim{a}) || isdatetime(axisLim{a})
            [tickVals, tickText] = getDateTicks(axisLim{a}, nTicks);
            expoFormat(a) = 0;

            [tickVals, tickText] = rmTicks(ax{a}.range, tickVals, tickText, 's');
        end

        ax{a}.showticklabels = true;
        ax{a}.ticks = 'inside';
        ax{a}.ticklen = tickLen;
        ax{a}.tickwidth = lineWidth;

        ax{a}.tickfont.size = fontSize;
        ax{a}.tickfont.family = fontFamily;
        ax{a}.tickfont.color = axisColor;

        ax{a}.tickvals = tickVals;
        ax{a}.ticktext = tickText;

        %-label-%
        if ~isempty(axisLabel{a})
            ax{a}.title = parseString(axisLabel{a});
            ax{a}.titlefont.color = axisColor;
            ax{a}.titlefont.size = fontSize;
            ax{a}.titlefont.family = fontFamily;
        end
    end

    %-------------------------------------------------------------------------%
end

function [tickVals, tickText] = getDateTicks(axisLim, nTicks)

    %-by year-%
    yearLim = year(date2NumData(axisLim));
    isYear = length(unique(yearLim)) > 1;
    refYear = [1, 2, 5];

    if isYear
        yearTick = getTickVals(yearLim, refYear, 1, nTicks);

        for n = 1:length(yearTick)
            tickVals(n) = datenum(datetime(yearTick(n),1,1,'Format','yy'));
            tickText{n} = num2str(yearTick(n));
        end

        return;
    end

    %-by month-%
    monthLim = month(date2NumData(axisLim));
    isMonth = length(unique(monthLim)) > 1;

    if isMonth
        % TODO
        return
    end

    %-by day-%
    dayLim = day(date2NumData(axisLim));
    isDay = length(unique(dayLim)) > 1;
    refDay = [0.5, 1];

    if isDay
        dayTick = getTickVals(dayLim, refDay, 1, nTicks);
        hourTick = 24 * ( dayTick - fix(dayTick) );
        dayTick = fix(dayTick);

        for n = 1:length(dayTick)
            matDate = [yearLim(1),monthLim(1),dayTick(n),hourTick(n),0,0];
            tickVals(n) = datetime(matDate,'Format', 'MMM dd, HH:mm');
            tickText{n} = datestr(tickVals(n), 'mmm dd, HH:MM');
        end

        tickVals = datenum(tickVals);
        return;
    end
end

function [tickVals, tickText, expoFormat] = getNumTicks(axisLim, nTicks)
    refVals = [1, 2, 5];
    refPot = floor(log10(range(axisLim)));
    nonStop = true;

    fixAxisLim = fix(axisLim);

    if ~all(fixAxisLim == 0)
        expoFormat = floor(log10(max(1, range(fixAxisLim))));
    else
        expoFormat = refPot;
    end

    if abs(expoFormat) < 3 expoFormat = 0; end

    tickVals = getTickVals(axisLim, refVals, refPot, nTicks);

    for n = 1:length(tickVals)
        tickText{n} = num2str(tickVals(n)/10^expoFormat);
    end
end

function tickVals = getTickVals(axisLim, refVals, refPot, nTicks)
    done = false;

    for p = refPot:-1:refPot-1
        for v = refVals
            vp = v*10^p;

            startTick = axisLim(1) - mod(axisLim(1), vp);
            if startTick < axisLim, startTick = startTick + vp; end
            endTick = axisLim(2) - mod(axisLim(2), vp);
            nPoints = floor(range([startTick, endTick])/vp) + 1;

            tickVals = startTick:vp:endTick;
            lenTicks = length(tickVals);

            if lenTicks >= nTicks(1) && lenTicks <= nTicks(2)
                done = true;
                break;
            end
        end

        if done, break; end
    end
end

function updateTitle(obj, titleText, xySource)
    xaxis = eval(sprintf('obj.layout.xaxis%d', xySource(1)));
    yaxis = eval(sprintf('obj.layout.yaxis%d', xySource(2)));
    anIndex = 1;
    if ~isempty(titleText), titleText = parseString(titleText); end

    ann.showarrow = false;
    ann.text = sprintf('%s', titleText);
    ann.xref = 'paper';
    ann.yref = 'paper';
    ann.x = mean(xaxis.domain);
    ann.y = yaxis.domain(2);
    ann.xanchor = 'middle';
    ann.yanchor = 'bottom';
    ann.font.size = 1.5*xaxis.tickfont.size;
    ann.font.color = xaxis.tickfont.color;
    ann.font.family = xaxis.tickfont.family;

    obj.layout.annotations{anIndex} = ann;
end

function updateExponentFormat(obj, expoFormat, xySource, axName)

    axName = lower(axName);
    xaxis = eval(sprintf('obj.layout.xaxis%d', xySource(1)));
    yaxis = eval(sprintf('obj.layout.yaxis%d', xySource(2)));
    anIndex = obj.PlotlyDefaults.anIndex + 1;

    if expoFormat ~= 0
        exponentText = sprintf('\\times10^%d', expoFormat);
        exponentText = parseString(exponentText, 'tex');
        % exponentText = ['
' exponentText '
'];

        ann.showarrow = false;
        ann.text = exponentText;
        ann.xref = 'paper';
        ann.yref = 'paper';
        ann.font.size = eval(sprintf('%saxis.tickfont.size', axName));
        ann.font.color = eval(sprintf('%saxis.tickfont.color', axName));
        ann.font.family = eval(sprintf('%saxis.tickfont.family', axName));
        ann.xanchor = 'left';

        switch axName
            case 'x'
                ann.yanchor = 'bottom';
                ann.x = xaxis.domain(2);
                ann.y = yaxis.domain(1);

            case 'y'
                ann.yanchor = 'bottom';
                ann.x = xaxis.domain(1);
                ann.y = yaxis.domain(2);
        end

        obj.layout.annotations{anIndex} = ann;
        obj.PlotlyDefaults.anIndex = anIndex;
    end
end

function [tickVals, tickText] = rmTicks(axisLim, tickVals, tickText, rmCase)
    rangeLim = range(axisLim);

    switch rmCase
        case 's'
            if abs(axisLim(1)-tickVals(1)) < rangeLim * 0.01
                tickVals = tickVals(2:end);
                tickText = tickText(2:end);
            end

        case 'e'
            if abs(axisLim(end)-tickVals(end)) < rangeLim * 0.01
                tickVals = tickVals(1:end-1);
                tickText = tickText(1:end-1);
            end

    end
end

author: plotly


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function [LiDARTag, AprilTag] = get4Corners(opt, mat_file_path, bagfile, target_len, pc_iter, num_scan)
    pc = loadPointCloud(mat_file_path, bagfile);
    X = getPayload(pc, pc_iter, num_scan);

    % cost
    opt = optimizeCost(opt, X, target_len);
    target_lidar = [0 -target_len/2 -target_len/2 1;
                    0 -target_len/2  target_len/2 1;
                    0  target_len/2  target_len/2 1;
                    0  target_len/2 -target_len/2 1]';
    LiDARTag.corners = inv(opt.H_opt) * target_lidar;
    LiDARTag.corners = sortrows(LiDARTag.corners', 3, 'descend')';
    LiDARTag.four_corners_line = point3DToLineForDrawing(LiDARTag.corners);
    LiDARTag.points = X;

%             HCamera = getAprilTagPose(app, app.image_num);
    AprilTag.corners = getAprilTagCorners(bagfile);
%             AprilTag.undistorted_corners = getApriltagUndistortedCorners(app, AprilTag.corners);
    AprilTag.four_corners_line = [];
end

author: UMich-BipedLab


text:
% Defines a stimulus
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 11/30/2009
classdef Stimuli < edu.stanford.covert.cell.kb.KnowledgeBaseObject
    properties
        compartments    = edu.stanford.covert.cell.kb.Compartment.empty(0, 0);
        values          = [];
        initialTimes    = [];
        finalTimes      = [];

        regulatedProteinMonomers = edu.stanford.covert.cell.kb.ProteinMonomer.empty(0, 0);
        regulatedProteinComplexs = edu.stanford.covert.cell.kb.ProteinComplex.empty(0, 0);
    end

    methods
        function this = Stimuli(knowledgeBase, wid, wholeCellModelID, name,...
                comments, crossReferences)

            if nargin == 0; return; end;

            this = edu.stanford.covert.cell.kb.Stimuli.empty(size(wid, 1), 0);
            if size(wid,1) == 0; return; end;

            this(size(wid, 1), 1) = edu.stanford.covert.cell.kb.Stimuli;
            for i = 1:size(wid, 1)
                this(i, 1).idx = i;
                this(i, 1).knowledgeBase = knowledgeBase;
                this(i, 1).wid = wid(i);
                this(i, 1).wholeCellModelID = wholeCellModelID{i};
                this(i, 1).name = name{i};
                if exist('comments', 'var') && ~isempty(comments); this(i,1).comments = comments{i}; end;
                if exist('crossReferences', 'var')
                    if size(crossReferences, 1) > 1
                        this(i, 1).crossReferences = crossReferences(i);
                    else
                        this(i, 1).crossReferences = struct;
                        fields = fieldnames(crossReferences);
                        for j = 1:size(fields,1)
                            vals = crossReferences.(fields{j});
                            this(i,1).crossReferences.(fields{j}) = vals(i);
                        end
                    end
                end
            end
        end

        function serializeLinks(this)
            for i = 1:numel(this)
                this(i).compartments             = this.serializeLinksHelper(this(i).compartments);
                this(i).values                   = this.serializeLinksHelper(this(i).values);
                this(i).initialTimes             = this.serializeLinksHelper(this(i).initialTimes);
                this(i).finalTimes               = this.serializeLinksHelper(this(i).finalTimes);

                this(i).regulatedProteinMonomers = this.serializeLinksHelper(this(i).regulatedProteinMonomers);
                this(i).regulatedProteinComplexs = this.serializeLinksHelper(this(i).regulatedProteinComplexs);

                serializeLinks@edu.stanford.covert.cell.kb.KnowledgeBaseObject(this(i));
            end
        end

        function deserializeLinks(this, kb)
            for i = 1:numel(this)
                this(i).compartments             = this.deserializeLinksHelper(this(i).compartments, kb.compartments);
                this(i).values                   = this.deserializeLinksHelper(this(i).values);
                this(i).initialTimes             = this.deserializeLinksHelper(this(i).initialTimes);
                this(i).finalTimes               = this.deserializeLinksHelper(this(i).finalTimes);

                this(i).regulatedProteinMonomers = this.deserializeLinksHelper(this(i).regulatedProteinMonomers, kb.proteinMonomers);
                this(i).regulatedProteinComplexs = this.deserializeLinksHelper(this(i).regulatedProteinComplexs, kb.proteinComplexs);

                deserializeLinks@edu.stanford.covert.cell.kb.KnowledgeBaseObject(this(i), kb);
            end
        end

        function deleteLinks(this)
            for i = 1:numel(this)
                this(i).compartments             = [];
                this(i).values                   = [];
                this(i).initialTimes             = [];
                this(i).finalTimes               = [];

                this(i).regulatedProteinMonomers = [];
                this(i).regulatedProteinComplexs = [];

                deleteLinks@edu.stanford.covert.cell.kb.KnowledgeBaseObject(this(i));
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [ a, l, r ] = r8row_part_quick_a ( m, n, a )

%*****************************************************************************80
%
%% R8ROW_PART_QUICK_A reorders the columns of an R8ROW.
%
%  Discussion:
%
%    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
%    each of length N.
%
%    The routine reorders the rows of A.  Using A(1,1:N) as a
%    key, all entries of A that are less than or equal to the key will
%    precede the key, which precedes all entries that are greater than the key.
%
%  Example:
%
%    Input:
%
%      M = 8, N = 2
%      A = ( 2 4
%            8 8
%            6 2
%            0 2
%           10 6
%           10 0
%            0 6
%            5 8 )
%
%    Output:
%
%      L = 2, R = 4
%
%      A = ( 0 2    LEFT
%            0 6
%            ----
%            2 4    KEY
%            ----
%            8 8    RIGHT
%            6 2
%           10 6
%           10 0
%            5 8 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the row dimension of A.
%
%    Input, integer N, the column dimension of A.
%
%    Input, real A(M,N).  On input, the array to be checked.
%
%    Output, real A(M,N), the reordered array.
%
%    Output, integer L, R, the indices of A that define the three
%    segments.  Let KEY = the input value of A(1:M,1).  Then
%    I <= L                 A(I,1:N) < KEY;
%         L < I < R         A(I,1:N) = KEY;
%                 R <= I    KEY < A(I,1:N).
%
  if ( m < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8ROWL_PART_QUICK_A - Fatal error!\n' );
    fprintf ( 1, '  M < 1.\n' );
    return
  end

  if ( m == 1 )
    l = 0;
    r = 2;
    return
  end

  key(1:n) = a(1,1:n);
  k = 1;
%
%  The elements of unknown size have indices between L+1 and R-1.
%
  l = 1;
  r = m + 1;

  for i = 2 : m

    if ( r8vec_gt ( n, a(l+1,1:n), key(1:n) ) )
      r = r - 1;
      vec = a(r,1:n);
      a(r,1:n) = a(l+1,1:n);
      a(l+1,1:n) = vec;
    elseif ( r8vec_eq ( n, a(l+1,1:n), key(1:n) ) )
      k = k + 1;
      vec = a(k,1:n);
      a(k,1:n) = a(l+1,1:n);
      a(l+1,1:n) = vec;
      l = l + 1;
    elseif ( r8vec_lt ( n, a(l+1,1:n), key(1:n) ) )
      l = l + 1;
    end

  end
%
%  Shift small elements to the left.
%
  for j = 1 : l - k
    a(j,1:n) = a(j+k,1:n);
  end
%
%  Shift KEY elements to center.
%
  for j = l - k + 1 : l
    a(j,1:n) = key(1:n);
  end
%
%  Update L.
%
  l = l - k;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = BEARapp()
% Launch the appropriate version of the BEAR app
if verLessThan('matlab','9.9')
    eval('app = BEARapp20a');
else
    eval('app = BEARapp21a');
end

if nargout == 1
    varargout = {app};
end

author: european-central-bank

repo: BEAR-toolbox


text:
function halton_test11 ( )

%*****************************************************************************80
%
%% TEST11 tests U2_TO_SPHERE_UNIT_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num2 = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST11\n' );
  fprintf ( 1, '  For the unit sphere in 3 dimensions:\n' );
  fprintf ( 1, '  U2_TO_SPHERE_UNIT_3D samples;\n' );

  dim_num = 2;
  halton_dim_num_set ( dim_num );
  n = 5;
  step = 123456789;
  halton_step_set ( step );
  seed(1:dim_num) = 0;
  halton_seed_set ( seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  DIM_NUM = %12d\n', dim_num );
  fprintf ( 1, '  N =    %12d\n', n );
  fprintf ( 1, '  STEP = %12d\n', step );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A few sample values:\n' );
  fprintf ( 1, '\n' );

  for j = 1 : n
    u = halton ( );
    x = u2_to_sphere_unit_3d ( u );
    fprintf ( 1, '  ' );
    for i = 1 : dim_num2
      fprintf ( 1, '  %6f', x(i) );
    end
    fprintf ( 1, '\n' );
  end

  n = 1000;
  step = 0;
  halton_step_set ( step );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  N =    %12d\n', n );
  fprintf ( 1, '  STEP = %12d\n', step );

  average(1:dim_num2) = 0.0;

  for i = 1 : n
    u = halton ( );
    x = u2_to_sphere_unit_3d ( u );
    average(1:dim_num2) = average(1:dim_num2) + x(1:dim_num2);
  end

  average(1:dim_num2) = average(1:dim_num2) / n;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average the points, which should get a value\n' );
  fprintf ( 1, '  close to zero, and closer as N increases.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Average:        ' );
  for i = 1 : dim_num2
    fprintf ( 1, '  %8e', average(i) );
  end
  fprintf ( 1, '\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now choose a random direction, sample the same\n' );
  fprintf ( 1, '  number of points, and compute the dot product with\n' );
  fprintf ( 1, '  the direction.\n' );
  fprintf ( 1, '  Take the absolute value of each dot product\n' );
  fprintf ( 1, '  and sum and average.\n' );

  for j2 = 1 : 5

    seed_val = get_seed ( );
    step = seed_val + 111 * j2;
    halton_step_set ( step );

    u = halton ( );
    v = u2_to_sphere_unit_3d ( u );

    step = 0;
    halton_step_set ( step );

    dot_average = 0.0;

    for j = 1 : n
      u = halton ( );
      x = u2_to_sphere_unit_3d ( u );
      dot_average = dot_average + abs ( x(1:dim_num2) * v(1:dim_num2)' );
    end

    dot_average = dot_average / n;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Random V:         ' );
    for j = 1 : dim_num2
      fprintf ( 1, '  %8f', v(j) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Average |(XdotV)| %f\n', dot_average );

  end

  return
end


text:
% nbins = 8;
% psize = 10;
% npatches = 8;
% optflowwinsig = 3;
% optflowsig = 2;
% optreliability = 1e-4;
% patchsz = psize*npatches;

params = getParams;
npatches = params.npatches;
psize = params.psize;
nbins = params.nbins;
patchsz = params.patchsz;
scale = params.scale;

optflowwinsig = params.optflowwinsig ;
optflowsig = params.optflowsig ;
optreliability = params.optreliability ;

maxflow = 5;
fly_thres = 90;
flow_thres = params.flow_thres;

makeVideo = false;

% Antennal grooming for andy
expdir = '/home/mayank/Dropbox/ForMayankFromAlice/grooming_GMR_30B01_AE_01_CsChr_RigB_20150903T161828/';
fly = 1;
frames = 5961:5980;

% expdir = '../mated10_20140714T131113';
% % frames = 10900 +(1:200);
% % fly = 1;
% fly = 1;
% % frames = 11043; % large middle leg movement
% frames = 10925; % large middle leg movement

% expdir = '/home/mayank/Work/FlySpaceTime/walkMovies/SS03500';
% fly = 1;
% frames = 10200 + (1:100);
% frames = 735:780;
% fly = 3;
% frames = 1886:1888;
% frames = 1477:1477;

% fly = 2; % This one shows the failure of HS-brightness on rear leg.
% frames = 5503;

[~,expname] = fileparts(expdir);
moviefilestr = 'movie.ufmf';
moviefile = fullfile(expdir,moviefilestr);
trxfilestr = 'trx.mat';
trxfile = fullfile(expdir,trxfilestr);

[readframe,nframes] = get_readframe_fcn(moviefile);
td = load(trxfile);
tracks = td.trx;
dx = []; dy = [];
theta = [];
rdx = [];rdy = [];
im1 = []; im2 = []; im3 = [];
for i = 1:numel(frames),
  curf = readframe(frames(i));
  nextf = readframe(frames(i)+1);

  trackndx = frames(i) - tracks(fly).firstframe + 1;
  locy = round(tracks(fly).y(trackndx));
  locx = round(tracks(fly).x(trackndx));
  curpatch = extractPatch(curf,...
    locy,locx,tracks(fly).theta(trackndx),patchsz);
  locy = round(tracks(fly).y(trackndx+1));
  locx = round(tracks(fly).x(trackndx+1));
  curpatch2 = extractPatch(nextf,...
    locy,locx,tracks(fly).theta(trackndx+1),patchsz);

  curpatch3 = extractPatch(nextf,...
    locy,locx,tracks(fly).theta(trackndx),patchsz);

  im1(:,:,i) = curpatch;
  im2(:,:,i) = curpatch2;
  im3(:,:,i) = curpatch3;
  dx(i) = round(tracks(fly).x(trackndx+1))-round(tracks(fly).x(trackndx));
  dy(i) = round(tracks(fly).y(trackndx+1))-round(tracks(fly).y(trackndx));
  rdx(i) = tracks(fly).x(trackndx+1)-tracks(fly).x(trackndx)-dx(i);
  rdy(i) = tracks(fly).y(trackndx+1)-tracks(fly).y(trackndx)-dy(i);
  dtheta(i) = tracks(fly).theta(trackndx+1)-tracks(fly).theta(trackndx);
  theta(i) = tracks(fly).theta(trackndx);
end

[nr,nc,~] = size(im1);
%
% figure out the histogram bins
tmptheta = (0:179)*pi/180;
res = nan(nbins,numel(tmptheta));
m = single(ones(10,10));
o = single(zeros(10,10));
for i = 1:numel(tmptheta),
%   fprintf('i = %d, theta = %f\n',i,tmptheta(i));
  o(:) = single(tmptheta(i));
  rescurr = gradientHist(m,o,1,nbins,1);
  res(:,i) = rescurr(1,1,:);
end

bincenters = nan(1,nbins);
for i = 1:nbins,
  bincenters(i) = tmptheta(argmax(res(i,:)));
end

% this seems to be what the centers correspond to
bincenters = linspace(0,pi,nbins+1);
bincenters = bincenters(1:nbins);

% mayank divides by 2
bincenters2 = bincenters*2;
dt = mean(diff(bincenters2));
binedges2 = [bincenters2(1)-dt/2,(bincenters2(1:end-1)+bincenters2(2:end))/2,bincenters2(end)+dt/2];

%% make a video

tr = 191;
hfig = 100;
t0 = frames(1);
t1 = frames(end);
colorpos = [1,0,0];
colorneg = [0,.3,1];

maxv2 = 0;
colors = hsv(nbins);
sz = size(im1);
bwimg = zeros(sz(1),sz(2));
ctr = [ceil( (sz(1)+1)/2),ceil( (sz(2)+1)/2)];
bwimg(ctr(1),ctr(2))=1;
dimg = bwdist(bwimg,'euclidean');
[xx,yy]= meshgrid(1:sz(2),1:sz(1));
aimg = atan2(-(yy-ctr(1)),-(xx-ctr(2)));

for t = t0:20:t1,

  im1curr = im1(:,:,t-t0+1);
  im2curr = im2(:,:,t-t0+1);
%   [Vx,Vy,~] = optFlowLk(im1(:,:,i),im2(:,:,i),[],optflowwinsig,optflowsig,optreliability);
%  [Vx,Vy,] = optFlowHorn(im1curr,im2curr,optflowsig);
  uv = estimate_flow_interface(im1curr,im2curr,'hs-brightness',{'max_warping_iters',10});
%   uv = estimate_flow_interface(im1curr,im2curr,'ba-brightness');
  Vx = uv(:,:,1); Vy = uv(:,:,2);

%   Vx = Vx-dx(t-t0+1);
%   Vy = Vy-dy(t-t0+1);

  M = sqrt(Vx.^2 + Vy.^2);
  O = mod(atan2(Vy,Vx)/2,pi);
  O = min(O,pi-1e-6);
  H = gradientHist(single(M),single(O),psize,nbins,1);
  maxv2 = max(maxv2,max(H(:)));
end

%%
if makeVideo
  vid = VideoWriter(sprintf('%s_Fly%d_From%d_To%d_HOF_%s.avi',expname,fly,t0,t1,datestr(now,'yyyymmdd')));
  open(vid);
end

dd_err = dimg/patchsz;
% There is some flow towards the center in certain cases.
% d_err is for that.

for t = t0:t1,

  im1curr = im1(:,:,t-t0+1);
  im2curr = im2(:,:,t-t0+1);

  imsz = size(im1curr);
  pairimg = zeros(2*imsz(1),3*imsz(2),3);
  ttimg = [];
  tt(:,:,1) = im2curr;
  tt(:,:,2) = im1curr;
  tt(:,:,3) = im2curr;
  pairimg(1:imsz(1),:,:) = repmat(tt,[1 3 1]);

%   [Vx,Vy,~] = optFlowLk(im1curr,im2curr,[],3);
%   [Vx,Vy,~] = optFlowLk(im1curr,im2curr,[],optflowwinsig,optflowsig,optreliability);
%   [Vx,Vy,] = optFlowHorn(im1curr,im2curr,optflowsig);
  uv = estimate_flow_interface(im1curr,im2curr,'hs-brightness',...
    {'max_warping_iters',2 });
%   uv = estimate_flow_interface(im1curr,im2curr,'ba-brightness',{'max_warping_iters',2 });
%   uv = estimate_flow_interface(im1curr,im2curr,'classic++');
  uvorig = uv;
  pairimg(imsz(1)+(1:imsz(1)),1:imsz(2),:) = flowToColor(uv,maxflow);

  cdx = dx(t-t0+1);
  cdy = dy(t-t0+1);
  curt = theta(t-t0+1);
  rotd = [cos(curt) sin(curt); -sin(curt) cos(curt)]*[cdx;cdy];
  cdx = -rotd(1); cdy = -rotd(2);

  ctheta = dtheta(t-t0+1);
  rotflowu = dimg.*(cos(aimg+ctheta)-cos(aimg));
  rotflowv = dimg.*(sin(aimg+ctheta)-sin(aimg));

  uv = uvorig;
  dd1 = sqrt( (uv(:,:,1)-cdx-rotflowu).^2 + (uv(:,:,2)-cdy-rotflowv).^2);
  for ndx = 1:2
    tt = uv(:,:,ndx);
    tt(dd1< (dd_err+flow_thres)) = 0;
    uv(:,:,ndx) = tt;
  end
  uvmotion = uv;
  pairimg(imsz(1)+(1:imsz(1)),imsz(2)+(1:imsz(2)),:) = flowToColor(uvmotion,maxflow);

  % Using bkg/fkg.
%   fly_bod = im1curr<fly_thres;
%   fly_flow = zeros(1,2);
%   for ndx = 1:2
%     tt = uv(:,:,ndx);
%     fly_flow(ndx) = median(tt(fly_bod));
%   end

  %using computation.
  crdx = rdx(t-t0+1);  crdy = rdy(t-t0+1);
  rotd = [cos(curt) sin(curt); -sin(curt) cos(curt)]*[crdx;crdy];
  fly_flow = rotd;

  uv = uvorig;
  dd2 = sqrt( (uv(:,:,1)-fly_flow(1)).^2 + (uv(:,:,2)-fly_flow(2)).^2);
  dd = min(dd1,dd2);
  for ndx = 1:2
    tt = uv(:,:,ndx);
    tt( dd<(dd_err+flow_thres)) = 0;
    uv(:,:,ndx) = tt;
  end
  uvall = uv;

  pairimg(imsz(1)+(1:imsz(1)),2*imsz(2)+(1:imsz(2)),:) = flowToColor(uvall,maxflow);

  pairimg = uint8(pairimg);
  pairimg = imresize(pairimg,scale);
  Vx = uvorig(:,:,1); Vy = uvorig(:,:,2);
  M = sqrt(Vx.^2 + Vy.^2);
  O = mod(atan2(Vy,Vx)/2,pi);
  O = min(O,pi-1e-6);
  Horig = gradientHist(single(M),single(O),psize,nbins,1);

  Vx = uvmotion(:,:,1); Vy = uvmotion(:,:,2);
  M = sqrt(Vx.^2 + Vy.^2);
  O = mod(atan2(Vy,Vx)/2,pi);
  O = min(O,pi-1e-6);
  Hmotion = gradientHist(single(M),single(O),psize,nbins,1);

  Vx = uvall(:,:,1); Vy = uvall(:,:,2);
  M = sqrt(Vx.^2 + Vy.^2);
  O = mod(atan2(Vy,Vx)/2,pi);
  O = min(O,pi-1e-6);
  Hall = gradientHist(single(M),single(O),psize,nbins,1);

  if t == t0,
    figure(hfig);
    clf;
    hax = axes('Position',[0,0,1,1]);
    set(hfig,'Units','pixels','Position',get(0,'ScreenSize'));

    him = imshow(pairimg);
    axis image;
    truesize;
    colormap gray;
    hold on;
    axis off;
    htext = text(1,nr-1,sprintf('%.2f s',(t-1)/200),'Color','k','HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',18);

  else
    delete(horig(ishandle(horig)));
    delete(hmotion(ishandle(hmotion)));
    delete(hall(ishandle(hall)));
    set(him,'CData',pairimg);
    set(htext,'String',sprintf('%.2f s',(t-1)/200));
  end

  horig = plotHofFeatures(Horig,hax,maxv2,colors,binedges2,psize*scale,[0 0]);
  hmotion = plotHofFeatures(Hmotion,hax,maxv2,colors,binedges2,psize*scale,[imsz(2) 0]*scale);
  hall = plotHofFeatures(Hall,hax,maxv2,colors,binedges2,psize*scale,[2*imsz(2) 0]*scale);

  drawnow;
  if makeVideo
    fr = getframe(hax);
    writeVideo(vid,fr);
  elseif t0<t1
    pause;
  end
end

if makeVideo,  close(vid); end

author: kristinbranson

repo: JAABA


text:
function vexity = DeriveVexityFromInflection(properties,xL,xU)

% [point (changes to convex = 1/changes to concave = -1) ... ]
if isa(properties.inflection,'function_handle')
    data = properties.inflection(xL,xU);
    if isempty(data)
        vexity = 'none';
        return
    end
    data = [data inf];
else
    data = properties.inflection;
    data = [data inf];
end
convex = [data(2:2:end-1)];
points = data(1:2:end);
vexity = 'none';
for k = 1:length(points)-1
    if xL >= points(k) && xU <= points(k+1)
        if convex(k) == 1
            vexity = 'convex';
            return
        else
            vexity = 'concave';
            return
        end
    end
end


text:
function [val,idx] = of_PCMARE(obs,sim,idx)
% of_MARE Calculates a version of the Mean Absolute Relative Error (MARE)
% of simulated streamflow as a percentage of the MARE of the mean
% observed flow. Ignores time steps with negative flow values.  Adds a
% constant e of 1/100 of mean(obs) to avoid issues with zero flows
% (Pushpalatha et al., 2012).

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% In:
% obs       - time series of observations       [nx1]
% sim       - time series of simulations        [nx1]
% idx       - optional vector of indices to use for calculation, can be
%               logical vector [nx1] or numeric vector [mx1], with m <= n
%
% Out:
% val       - objective function value          [1x1]
% idx       - vector of indices used for calculation

% Pushpalatha, R.; Perrin, C.; le Moine, N. and Andréassian V. (2012). "A
% review of efficiency criteria suitable for evaluating low-flow
% simulations". Journal of Hydrology. 420-421, 171-182.
% doi:10.1016/j.jhydrol.2011.11.055

%% Check inputs and select timesteps
if nargin < 2
    error('Not enugh input arguments')
end

if nargin < 3; idx = []; end
[sim, obs, idx] = check_and_select(sim, obs, idx);

%% Find the constant e
m = mean(obs);
e = m/100;

%% Apply constant and transform flows
obs = obs+e;
sim = sim+e;

%% Calculate metric
MARE = mean(abs((sim-obs)/obs));
MARE_mean = mean(abs((sim-m)/m));

val = MARE/MARE_mean;
end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function d = amsr(im,sigma1,sigma2,sigma3)

if ~exist('sigma1','var'),sigma1 = 15;end
if ~exist('sigma2','var'),sigma2 = 80;end
if ~exist('sigma3','var'),sigma3 = 250;end

im = im2double(im); % double --> im2double
gausKernel1 = fspecial('gaussian',[sigma1 sigma1],5);
gausKernel2 = fspecial('gaussian',[sigma2 sigma2],20);
gausKernel3 = fspecial('gaussian',[sigma3 sigma3],50);
Y = 0.299.*im(:,:,1)+0.587.*im(:,:,2)+0.114.*im(:,:,3);

blur_im1 = (imfilter(Y,gausKernel1,'replicate'));
blur_im2 = (imfilter(Y,gausKernel2,'replicate'));
blur_im3 = (imfilter(Y,gausKernel3,'replicate'));

Y_ssr1 = log(Y) - log(blur_im1);
Y_ssr2 = log(Y) - log(blur_im2);
Y_ssr3 = log(Y) - log(blur_im3);

sr = sort(Y_ssr1(:));
p1 = sr(uint32(0.01*size(im,1)*size(im,2)));
p99 = sr(uint32(0.99*size(im,1)*size(im,2)));
lc1 = Y_ssr1 >= p1 & Y_ssr1 <= p99;
lc2 = Y_ssr1 < p1;
lc3 = Y_ssr1 > p99;
Y_ssr1(lc1) = 255 * (Y_ssr1(lc1)-p1)./(p99-p1);
Y_ssr1(lc2) = 0;
Y_ssr1(lc3) = 255;

sr = sort(Y_ssr2(:));
p1 = sr(uint32(0.01*size(im,1)*size(im,2)));
p99 = sr(uint32(0.99*size(im,1)*size(im,2)));
lc1 = Y_ssr2 >= p1 & Y_ssr2 <= p99;
lc2 = Y_ssr2 < p1;
lc3 = Y_ssr2 > p99;
Y_ssr2(lc1) = 255 * (Y_ssr2(lc1)-p1)./(p99-p1);
Y_ssr2(lc2) = 0;
Y_ssr2(lc3) = 255;

sr = sort(Y_ssr3(:));
p1 = sr(uint32(0.01*size(im,1)*size(im,2)));
p99 = sr(uint32(0.99*size(im,1)*size(im,2)));
lc1 = Y_ssr3 >= p1 & Y_ssr3 <= p99;
lc2 = Y_ssr3 < p1;
lc3 = Y_ssr3 > p99;
Y_ssr3(lc1) = 255 * (Y_ssr3(lc1)-p1)./(p99-p1);
Y_ssr3(lc2) = 0;
Y_ssr3(lc3) = 255;

% Y = 0.299.*im(:,:,1)+0.587.*im(:,:,2)+0.114.*im(:,:,3);
% Y_ssr1 = 0.299.*R_ssr1(:,:,1)+0.587.*R_ssr1(:,:,2)+0.114.*R_ssr1(:,:,3);
% Y_ssr2 = 0.299.*R_ssr2(:,:,1)+0.587.*R_ssr2(:,:,2)+0.114.*R_ssr2(:,:,3);
% Y_ssr3 = 0.299.*R_ssr3(:,:,1)+0.587.*R_ssr3(:,:,2)+0.114.*R_ssr3(:,:,3);
sigma = 32;
mu0 = 32;
mu1 = 96;
mu2 = 160;
mu3 = 224;
P0 = exp(-(Y-mu0).^2./(2*sigma^2));
P1 = exp(-(Y-mu1).^2./(2*sigma^2));
P2 = exp(-(Y-mu2).^2./(2*sigma^2));
P3 = max(P0,exp(-(Y-mu3).^2./(2*sigma^2)));
psum = P0+P1+P2+P3;
omega0 = P0./psum;
omega1 = P1./psum;
omega2 = P2./psum;
omega3 = P3./psum;
Y_amsr = omega0.*Y + omega1.*Y_ssr1 + omega2.*Y_ssr2 + omega3.*Y_ssr3;
ratio = Y_amsr./Y;

% modified: uint8 --> im2uint8
hsv = rgb2hsv(im2uint8(im));
v = hsv(:,:,3);%.*255;
v = 0.5.*(ratio.*(v+Y)+v-Y);
hsv(:,:,3) = v./255;
%d = uint8(hsv2rgb(hsv).*255);
d = double(hsv2rgb(hsv)); % modified
d = min(max(d,0),1); % add

% r = 0.5.*(ratio.*(im(:,:,1)+Y)+im(:,:,1)-Y);
% g = 0.5.*(ratio.*(im(:,:,2)+Y)+im(:,:,2)-Y);
% b = 0.5.*(ratio.*(im(:,:,3)+Y)+im(:,:,3)-Y);
% d = cat(3,r,g,b)./255;
end

author: dawnlh

repo: awesome-low-light-image-enhancement


text:
classdef AGEMOEA < ALGORITHM
% <multi/many> <real/integer/label/binary/permutation> <constrained/none>
% Adaptive geometry estimation-based many-objective evolutionary algorithm

%------------------------------- Reference --------------------------------
% A. Panichella, An adaptive evolutionary algorithm based on non-euclidean
% geometry for many-objective optimization, Proceedings of the Genetic and
% Evolutionary Computation Conference, 2019, 595-603.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Annibale Panichella

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();
            [~,FrontNo,CrowdDis] = EnvironmentalSelection(Population,Problem.N);

            %% Optimization
            while Algorithm.NotTerminated(Population)
              MatingPool = TournamentSelection(2,Problem.N,FrontNo,-CrowdDis);
              Offspring  = OperatorGA(Problem,Population(MatingPool));
              [Population,FrontNo,CrowdDis] = EnvironmentalSelection([Population,Offspring],Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [P,Y] = preprocess_data(Y,p,options)
% PREPROCESS_DATA - Preprocess data for CNMF analysis of image signals
%
%   [P_PARMS,Y] = PREPROCESS_DATA(Y, P_ORDER, OPTIONS)
%
% Inputs:
%    Y - the image data to be examined
%    P_ORDER - The order of the autoregressive system
%    OPTIONS - A structure of options (see help CNMFSetParms)
%
% Outputs:
%    P_PARMS - A structure of parameters extracted from the pre-processed data
%    P_PARMS has the following fields:
%
%       Fieldname:                       | Description
%       --------------------------------------------------------------------------
%       p                                |  P_ORDER as above
%       mis_data                         |  Index locations of missing (NaN) values in Y
%       mis_values                       |  Interpolated values at the missing index locations
%       pixels                           |  Index values of pixels that are not saturated
%       sn                               |  Noise power of each pixel
%       g                                |  Autoregressive model parameters
%
%    Y - A processed version of the data Y with the following changes:
%      (i)   identifying and interpolating missing entries (assumed to have the
%            value NaN). Interpolated entries are passed back to Y.
%      (ii)  identifying saturated pixels
%      (iii) estimating noise level for every pixel
%      (iv)  estimating global discrete time constants (if needed)
%
% This function replaces ARPFIT, present in the previous versions of the code.
%

% Author: Eftychios A. Pnevmatikakis
%           Simons Foundation, 2015

defoptions = CNMFSetParms;

if nargin < 3 || isempty(options); options = defoptions; end
if nargin < 2 || isempty(p); p = 2; end
P.p = p;

if ~isfield(options,'noise_range'); options.noise_range = defoptions.noise_range; end
if ~isfield(options,'noise_method'); options.noise_method = defoptions.noise_method; end
if ~isfield(options,'block_size'); options.block_size = defoptions.block_size; end
if ~isfield(options,'flag_g'); options.flag_g = defoptions.flag_g; end
if ~isfield(options,'lags'); options.lags = defoptions.lags; end
if ~isfield(options,'include_noise'); options.include_noise = defoptions.include_noise; end; include_noise = options.include_noise;
if ~isfield(options,'split_data'); split_data = defoptions.split_data; else split_data = options.split_data; end
if ~isfield(options,'cluster_pixels'); cluster_pixels = defoptions.cluster_pixels; else cluster_pixels = options.cluster_pixels; end
if ~isfield(options,'extract_max'); extract_max = defoptions.extract_max; else extract_max = options.extract_max; end
if ~isfield(options,'max_nlocs'); options.max_nlocs = defoptions.max_nlocs; end
if ~isfield(options,'max_width'); options.max_width = defoptions.max_width; end

%% interpolate missing data

if any(isnan(Y(:)))
    Y_interp = interp_missing_data(Y);      % interpolate missing data
    mis_data = find(Y_interp);
    Y(mis_data) = Y_interp(mis_data);       % introduce interpolated values for initialization
else
    Y_interp = sparse(size(Y));
    mis_data = [];
end
P.mis_values = full(Y_interp(mis_data));
P.mis_entries = mis_data;

%% indentify saturated pixels

P.pixels = find_unsaturatedPixels(Y);                % pixels that do not exhibit saturation

%% estimate noise levels

fprintf('Estimating the noise power for each pixel from a simple PSD estimate...');
[sn,psx] = get_noise_fft(Y,options);
P.sn = sn(:);
fprintf('  done \n');

%% cluster pixels based on PSD
if cluster_pixels
    psdx = sqrt(psx(:,3:end-1));
    X = psdx(:,1:min(size(psdx,2),1500));
    P.psdx = X;
    X = bsxfun(@minus,X,mean(X,2));     % center
    X = bsxfun(@times,X,1./sqrt(mean(X.^2,2)));
    [L,Cx] = kmeans_pp(X',2);
    [~,ind] = min(sum(Cx(max(1,end-49):end,:),1));
    P.active_pixels = (L==ind);
    P.centroids = Cx;

    if (0) % not necessary at the moment
        %[P.W,P.H] = nnmf(psdx,2); %,'h0',H0);
        psdx = psdx(:,1:min(size(psdx,2),600));
        r = sort(rand(1,size(psdx,2)),'descend');
        er = ones(1,length(r))/sqrt(length(r));
        H = [r/norm(r); er];
        W_ = rand(size(psdx,1),2);
        for iter = 1:100
            W = max((H*H')\(H*psdx'),0)';
            %H = max((W'*W)\(W'*psdx),0);
            r = max((W(:,1)'*psdx - (W(:,1)'*W(:,2))*er)/norm(W(:,1))^2,0);
            H = [r/norm(r); er];
            if norm(W-W_,'fro')/norm(W_,'fro') < 1e-3
                break;
            else
                W_ = W;
            end
        end
        disp(iter)
        W = max((H*H')\(H*psdx'),0)';
        P.W = W;
        P.H = H;
    end
end

%% extract maximum activity for each pixel
if extract_max
    [LOCS,Ym] = extract_max_activity(Y,options.max_nlocs,options.max_width);
    P.max_locs = LOCS;
    P.max_data = Ym;
end
%% estimate global time constants

if options.flag_g
    if ndims(Y) == 3
        Y = reshape(Y,size(Y,1)*size(Y,2),size(Y,3));
    end

    ff = options.pixels;
    np = length(ff);

    fprintf('Estimating time constant through autocorrelation function.. \n');
    tt1 = tic;
    mp = max(p);
    lags = options.lags + mp;
    if split_data
        Ycl = mat2cell(Y(ff,:),ones(np,1),size(Y,2));
        XC = cell(np,1);
        parfor j = 1:np
            XC{j} = xcov(Ycl{j},lags,'biased');
        end
        XC = cell2mat(XC);
    else
        XC = zeros(np,2*lags+1);
        for j = 1:np
            XC(j,:) = xcov(Y(ff(j),:),lags,'biased');
        end
    end

    gv = zeros(np*lags,1);
    if ~include_noise
        g =  XC(:,lags:-1:1);
        clear XC;
        lags = lags - p;
    end
    A = zeros(np*lags,p);
    for i = 1:np
        if ~include_noise
            A((i-1)*lags + (1:lags),:) = toeplitz(g(i,p:p+lags-1),g(i,p:-1:1));
        else
            A((i-1)*lags + (1:lags),:) = toeplitz(XC(i,lags+(1:lags)),XC(i,lags+(1:p))) - sn(i)^2*eye(lags,p);
            gv((i-1)*lags + (1:lags)) = XC(i,lags+2:end)';
        end
    end
    if ~include_noise
        gv = g(:,p+1:end)';
    end
    %ph = pinv(A)*gv(:);
    ph = A\gv(:);
    disp(ph);
    fprintf('Done after %2.2f seconds. \n',toc(tt1));
    P.g = ph(:);
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
%function [meanmap,pmap]=randperm(srp_data,compmat,numperm,corp)

function [realcomp]=randperm(srp_data,compmat,corp)

if nargin <2
    error('need to input raw data, comparison matrix,and p level');
end

if nargin <3
    corp=0.05;
end

a=size(srp_data);
if length(a)>2;
    error('for 2d only. use anovacont() or permtest()');
end

numsub=size(srp_data,1);
numcond=size(srp_data,2);

%disp('NOT correcting across the third dimension');

srp_data=srp_data(:,:,:,:);
%disp([(numsub*factorial(numcond)),' of your permutations will be the same as point estimate']);

%%calculate degrees of freedom
df_tot=(numsub*numcond)-1;
df_group=numcond-1;
df_subs=numsub-1;
df_error=df_tot-df_group-df_subs;

%    disp('calculating critical F....');

    smtot=squeeze(sum(srp_data,1));               %sum of each condition
    smtot1=squeeze(sum(srp_data,2));
    smmean=squeeze(mean(srp_data,1));          %mean of each condition
    ss_srp_data=srp_data.^2;                     %matrix of squared entries
    g_srp_data=squeeze(sum(smtot,2));        % this is G, the sum of all entries
    sum_ss_srp_data=sum(ss_srp_data(:)); %sum of squares
    sstot=sum_ss_srp_data-((g_srp_data.^2)/((numcond*numsub)));  %ss total
    ssgroup=((squeeze(sum(smtot.^2)))/(numsub))-((g_srp_data.^2)/(numcond*numsub));   %ss group
    sssubs=((squeeze(sum(smtot1.^2)))/numcond)-((g_srp_data.^2)/(numcond*numsub));   %ss group
    sserror=sstot-ssgroup-sssubs;
    msgroup=ssgroup/df_group;
    mserror=sserror/df_error;  %df error=36: 56(dftot)-2(dfgroup)-18(dfsubjects)
    realF=msgroup./(mserror+0.00001);         %this is the matrix of critical F values for frequency*time
    for c=1:size(compmat,1);
        for w=1:numcond;
        smtot2(w)=smmean(w)*compmat(c,w);
        end;
        comp1=squeeze(sum(smtot2));
        estvar1=(mserror./numsub)*(sum(compmat(c,:).^2));
        realcomp=comp1./sqrt(estvar1+0.000001);    %critical F values for frequency*time for each planned comparison
    end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function  vemPoisson3(cube, pde, option, varargin)

maxIt = option.maxIt;

%% Initialize err
errL2 = zeros(maxIt,1);   errH1 = zeros(maxIt,1);
erruIuh = zeros(maxIt,1); errMax = zeros(maxIt,1);
errTime = zeros(maxIt,1); solverTime = zeros(maxIt,1);
assembleTime = zeros(maxIt,1); meshTime = zeros(maxIt,1);
itStep = zeros(maxIt,1);  stopErr = zeros(maxIt,1); flag = zeros(maxIt,1);
N = zeros(maxIt,1);
h = zeros(maxIt,1);
h(1) = option.h0;

% [node,elem, interfaceData] = interfacemesh3(cube, pde.phi, h(1));
% checkinterfacemesh3(node,elem,interfaceData)
% [face, face2elem] = getpolymesh(elem,interfaceData);
% clear elem interfaceData

%         [node, elem] = cubehexmesh(cube,h(k+1));
%         idx = (1:size(elem, 1))';
%         face =  [elem(:,[1,4,3,2]); elem(:,[6, 7, 8, 5]);...
%                  elem(:,[2,3,7,6]); elem(:,[1, 5, 8, 4]);...
%                  elem(:,[1,2,6,5]); elem(:,[4, 8, 7, 3])];
%         face2elem = [idx;idx;idx;idx;idx;idx];
%         clear elem;

[node, elem] = cubemesh(cube, h(1));
elem = fixorder3(node,elem);
idx = (1:size(elem,1))';
face = [elem(:,[2, 3, 4]);elem(:,[1, 4,3]);elem(:,[1,2,4]);elem(:,[1,3,2])];
face(:,4) = 0;
face2elem = [idx;idx;idx;idx];
clear elem

%% Vertual Element Method
for k = 1:maxIt
    N(k) = size(node,1);
    [u, info] = Poisson3VEM(node, face, face2elem, pde, option);
    errMax(k) = max(abs(u - pde.exactu(node)));
    if k < maxIt
          h(k+1) = h(k)/2;

%         [node, elem, interfaceData] = interfacemesh3(cube, pde.phi, h(k+1));
%         h(k+1)
%         checkinterfacemesh3(node,elem,interfaceData)
%         [face, face2elem] = getpolymesh(elem,interfaceData);
%         clear elem interfaceData

%         [node, elem] = cubehexmesh(cube,h(k+1));
%         idx = (1:size(elem, 1))';
%         face =  [elem(:,[1,4,3,2]); elem(:,[6, 7, 8, 5]);...
%                  elem(:,[2,3,7,6]); elem(:,[1, 5, 8, 4]);...
%                  elem(:,[1,2,6,5]); elem(:,[4, 8, 7, 3])];
%         face2elem = [idx;idx;idx;idx;idx;idx];
%         clear elem;

        [node, elem] = cubemesh(cube, h(k+1));
        elem = fixorder3(node,elem);
        idx = (1:size(elem,1))';
        face = [elem(:,[2, 3, 4]);elem(:,[1, 4, 3]);...
                elem(:,[1, 2, 4]);elem(:,[1, 3, 2])];
        face(:,4) = 0;
        face2elem = [idx;idx;idx;idx];
        clear elem
    end
end

set(gcf, 'Units', 'normal');
set(gcf, 'Position', [0.25, 0.25, 0.55, 0.4]);
showrate(1./h(1:k), errMax(1:k), 1, 'k-+', '||u-u_h||');

end

repo: ifem


text:
function progFigure = progressdlg(varargin)
%progressdlg show/update progress dialog
%
% progressdlg creates or updates an existing progress dialog that reports on the
% progress of a task including the estimated time remaining as well as a button to
% cancel the task. It does not use drawnow in order to improve performance when using
% one or several other figures (drawnow refreshes all graphics objects in all figures).
%
% Syntax:
%   progressdlg                 creates default progress dialog
%
%   progressdlg('message');     creates default progress dialog with custom message
%
%   progressdlg(n);             creates default progress dialog with initial value n
%                               or updates existing progress dialog to value n
%
%   progressdlg(n,'message');   creates default progress dialog with message and value n
%                               or updates existing progress dialog to value n and message
%
%   progressdlg('PropertyName',PropertyValue,...);  creates or updates progress dialog with
%                                                   properties PropertyName/PropertyValue
%
%   handle = progressdlg(...);  returns the handle to the progress dialog
%                               handle is empty if user aborted by pressing cancel
%
%  PropertyName  | Description                                | Type           | Default
% ---------------|--------------------------------------------|----------------|------------------
%  Cancel        | Displays a cancel button to abort          | 'on,off'       | 'off'
%  Indeterminate | Sets the ProgressBar to indeterminate      | 'on,off'       | 'off'
%  Min           | Minimum value of ProgressBar               | double         | 0
%  Max           | Minimum value of ProgressBar               | double         | 100
%  Parent        | Parent Figure for Positioning              | figure_handle  | []
%  Position      | Override position of dialog                | [int x, int y] | []
%  ShowTime      | Displays an estimate of the remaining time | 'on/off'       | 'on'
%  Size          | Factor for dialog size (width=400px)       | double         | 1
%  String        | Descriptive text above ProgressBar         | char           | 'Progress'
%  Title         | Title of the dialog                        | char           | 'Progress'
%  Value         | Current Value of ProgressBar               | double         | 0
%  WindowStyle   | Normal or modal window behavior            | 'normal/modal' | 'modal'
%
%
% Notes:
%   1) Does not use drawnow to refresh the progress bar in order to improve performance,
%      especially when using one or several additional figures with many graphics objects.
%   2) Performance is slightly decreased when using the remaining time estimation and/or
%      cancel button
%
% Examples:
%   progressdlg('Calculating magic square');
%   for n=1:100
%       magic(1000);
%       progressdlg(n);
%   end
%
%   progressdlg('Title','Progress Dialog','String','Calculating magic squares','Min',0,'Max',1000,'Value',0);
%   for n=1:1000
%       magic(400);
%       progressdlg(n);
%   end
%
%   progressdlg('Title','Progress Dialog','String','Calculating magic squares','Position',[100 100],'WindowStyle','normal');
%   for n=1:10
%       magic(10000);
%       progressdlg(n*10,['Calculating magic square - Round: ' num2str(n)]);
%   end
%
%   h = progressdlg('Title','Progress','String','Calculating magic square','ShowTime','on','Cancel','on','Size',1.5);
%   for n=1:100
%       magic(4000);
%       h = progressdlg(n);
%       if isempty(h)
%           break;
%       end
%   end
%
% Warning:
%   This code heavily relies on undocumented and unsupported Matlab
%   functionality. It works on Matlab 7+, but use at your own risk!
%
% Bugs and suggestions:
%   Please send to Felix Ruhnow (felix.ruhnow at tu-dresden dot de)
%
% See also: WAITBAR, TIMEBAR, PROGRESSBAR (on the File Exchange)
%
% Adapted from:
% Daniel Claxton's WORKBAR
%
% Inspired by:
% Yair M. Altman's STATUSBAR
%
% Release history:
%    1.0  2013-06-10: initial version
%
% License to use and modify this code is granted freely to all interested, as long as the original author is
% referenced and attributed as such. The original author maintains the right to be solely associated with this work.
%
% Programmed and Copyright by Felix Ruhnow: felix.ruhnow(at)tu-dresden.de

    % keep variables for next call of progressdlg
    persistent p min max cancel progHandles showtime starttime;
    message = 'Progress';
    value = 0;
    text = [];
    est_text = 'Estimated time remaining: ';

    % check input arguments
    if nargin < 1
        % progressdlg; creates default progress dialog
        progHandles = [];

    elseif nargin == 1
        if ischar(varargin{1})
            % progressdlg('close'); closes progress dialog
            if strcmp(varargin{1},'close')
                if isfield(progHandles,'figure') && ishandle(progHandles.figure)
                    close(progHandles.figure);
                end
                progHandles = [];
                progFigure =[];
                return
            else
                % progressdlg('message'); creates default progress dialog with custom message
                message = varargin{1};
                progHandles = [];
            end

        elseif isnumeric(varargin{1})
            % progressdlg(n); creates default progress dialog with initial value n
            %                 or updates existing progress dialog to value n
            value = varargin{1};
        end
        varargin(:) = [];

    elseif nargin == 2 && isnumeric(varargin{1})
        % progressdlg(n,'message');
        if isempty(progHandles)
            % creates default progress dialog with message and value n
            value = varargin{1};
            message = varargin{2};
            varargin(:) = [];
        else
            % updates existing progress dialog to value n and message
            value = varargin{1};
            text = varargin{2};
        end

    elseif isfield(progHandles,'figure') && ishandle(progHandles.figure)
        % progressdlg('PropertyName',PropertyValue,...); updates progress dialog
        parse(p,varargin{1:end});
        value = p.Results.value;
        text = p.Results.string;

    else
        % progressdlg('PropertyName',PropertyValue,...); creates progress dialog
        progHandles = [];
    end

    if isempty(progHandles)
        % if handles are empty create new dialog
        cancel = 0;
        showtime = 0;

        % get initial input arguments
        p = checkInput(value,message);
        parse(p,varargin{1:end});
        prog = p.Results;
        min = prog.min;
        max = prog.max;

        % define height w/o estimated time remaining or cancel button
        height = 90;
        if strcmp(prog.showtime,'on') && ~strcmp(prog.indeterminate,'on')
            showtime = 1;
            height = height + 20;
        end
        if strcmp(prog.cancel,'on')
            cancel = 1;
            height = height + 30;
        end

        % get position of progress dialog (middle of screen, if not otherwise)
        if isempty(prog.position)
            if isempty(prog.parent)
                screen = get(0,'ScreenSize');
                pos = ceil([(screen(3)-400*prog.size)/2 (screen(4)-height*prog.size)/2 400 height]);
            else
                units = get(prog.parent,'Units');
                set(prog.parent,'Units','Pixel');
                pos_parent = get(prog.parent,'Position');
                set(prog.parent,'Units',units);
                pos = ceil([pos_parent(1)+(pos_parent(3)-400*prog.size)/2 pos_parent(2)+(pos_parent(4)-height*prog.size)/2 400 height]);
            end
        else
           pos = ceil([prog.position 400 height]);
        end

        % create dialog dialog
        progFigure = figure('Menubar','none','Numbertitle','off','Resize','off','Tag','ProgressDlg',...
                            'Units','pixels','Visible','off','DockControls','off','IntegerHandle','off',...
                            'Name',prog.title,'Position',pos,'WindowStyle','normal');%prog.windowstyle);

        % get background color of figure to adjust background of its children
        uiColor = get(progFigure,'Color');
        jColor = java.awt.Color(uiColor(1),uiColor(2),uiColor(3));

        if ispc
            barsize = 20;
            buttonsize = 25;
            buttonfont = 11;
        else
            barsize = 40;
            buttonsize = 30;
            buttonfont = 14;
        end

        % create descriptive text from PropertyName 'String', use java for updating dialog
        text1 = javaObjectEDT('javax.swing.JTextArea');
        [jText(1),cText(1)] = javacomponent(text1,[20 height-30 300 20],progFigure);
        set(jText(1),'Background',jColor,'Text',prog.string,...
                     'Font',java.awt.Font('Helvetica',java.awt.Font.PLAIN,14*prog.size));
        set(cText(1),'units','normalized');

        % create text for estimated remaining time if 'ShowTime' is 'on', use java for updating dialog
        if showtime
            text2 = javaObjectEDT('javax.swing.JTextArea');
            [jText(2),cText(2)] = javacomponent(text2,[20 height-90 300 20],progFigure);
            set(jText(2),'Background',jColor,'Text',est_text,...
                         'Font',java.awt.Font('Helvetica',java.awt.Font.PLAIN,14*prog.size));
            set(cText(2),'Units','Normalized');

        end

        % create cancel button if 'Cancel' is 'on', Click on button closes the dialog
        if cancel
            button = uicontrol(progFigure,'Style','pushbutton','Position',[120 10 160 buttonsize],...
                                          'String','cancel','Callback','close','FontSize',buttonfont*prog.size);
            set(button,'Units','Normalized');
        end

        % create ProgressBar with java JProgressBar class
        jProgressBar = javaObjectEDT('javax.swing.JProgressBar');
        jProgressBar.setStringPainted(true);
        set(jProgressBar,'Background',uiColor,'Minimum',prog.min,'Maximum',prog.max,'Value',prog.value);
        if strcmp(prog.indeterminate,'on')
            jProgressBar.setIndeterminate(true);
            showtime = 0;
        end
        [hProgressBar,cProgressBar] = javacomponent(jProgressBar,[20 height-65 360 barsize],progFigure);
        set(cProgressBar,'units','normalized','Tag','ProgressBar');

        % resize progress dialog according to PropertyName 'Size'
        set(progFigure,'Position',[pos(1:2) pos(3:4)*prog.size],'Visible','on');

        % save handles of the created objects
        progHandles.figure = progFigure;
        progHandles.progress = hProgressBar;
        progHandles.text = jText;

        % update dialog (only once!)
        drawnow;

        %save start time for progress
        starttime = clock;
    else
        % if handles are present update progress dialog
        progFigure = progHandles.figure;

        % update ProgressBar
        set(progHandles.progress,'Value',value);

        % update descriptive text if required
        if ~isempty(text)
            set(progHandles.text(1),'Text',text);
        end

        % calculate and update remaining time if required
        if showtime
            if value == min
                starttime=clock;
            else
                runtime = etime(clock,starttime);
                timeleft = runtime*(max-min)/(value-min) - runtime;
                timeleftstr = sec2timestr(timeleft);
                set(progHandles.text(2),'Text',[est_text timeleftstr]);
            end
        end

        % check event queue for callback of the 'cancel' button
        if cancel
            pause(10^-12); %flush event queue
            if ~ishandle(progHandles.figure)
                % if figure was delete via 'close' callback of 'cancel' button
                progHandles = [];
                progFigure = [];
                return;   % return empty handle
            end
        end
    end

    if value==max
        % if 'Maximum' value reached close and reset progress dialog
        close(progHandles.figure);
        progHandles =[];
        progFigure = [];
    end

%% Define input parameters and their default values
function p = checkInput(value,message)
    p = inputParser;
    states = {'on','off'};
    style = {'normal','modal'};
    addParamValue(p,'title','Progress');
    addParamValue(p,'size',1);
    addParamValue(p,'value',value);
    addParamValue(p,'min',0);
    addParamValue(p,'max',100);
    addParamValue(p,'position',[],@isnumeric);
    addParamValue(p,'string',message);
    addParamValue(p,'parent',[],@ishandle);
    addParamValue(p,'cancel','off',@(x) any(validatestring(x,states)))
    addParamValue(p,'showtime','on',@(x) any(validatestring(x,states)))
    addParamValue(p,'indeterminate','off',@(x) any(validatestring(x,states)))
    addParamValue(p,'windowstyle','modal',@(x) any(validatestring(x,style)))

%% Convert seconds to hh:mm:ss
function timestr = sec2timestr(sec)
    h = floor(sec/3600); % Hours
    sec = sec - h*3600;
    m = floor(sec/60); % Minutes
    sec = sec - m*60;
    s = floor(sec); % Seconds
    timestr = sprintf('%02d:%02d:%02d',h,m,s);


text:
%RANDOMFORESTC Breiman's random forest
%
%   W = RANDOMFORESTC(A,L,N)
%   W = A*RANDOMFORESTC([],L,N)
%
% INPUT
%   A       Dateset used for training
%   L       Number of decision trees to be generated (default 50)
%   N       Size of feature subsets to be used (default 1)
%
% OUTPUT
%   W       Resulting, trained feature space classifier
%
% DESCRIPTION
% Train a decision forest on A, using L decision trees, each trained on
% a bootstrapped version of dataset A. Each decison tree is using random
% feature subsets of size N in each node.  When N=0, no feature subsets
% are used.
%
% REFERENCES
% [1] L. Breiman, Random Forests, Machine learning, vol. 45 (1), 5-32, 2001
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, DTC

% Copyright: D.M.J. Tax, D.M.J.Tax@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function out = randomforestc(varargin)

argin = setdefaults(varargin,[],50,1);
if mapping_task(argin,'definition')

  out = define_mapping(argin,'untrained',['RandForest' int2str(argin{2})]);

elseif mapping_task(argin,'training')

  [a,L,featsubset] = deal(argin{:});
	isvaldfile(a,2,2); % at least 2 obj/class, 2 classes
	opt = [];
	[n,dim,opt.K] = getsize(a);
	opt.featsubset = featsubset;
	v = cell(L,1);
	for i=1:L
		[x,z] = gendat(a);
    if exist('decisiontree','file')==3
      v{i} = decisiontree(+x,getnlab(x),opt.K,opt.featsubset);
    else
		 prwarning(2,'No compiled decisiontree found, using the slower Matlab implementation.');
    	v{i} = tree_train(+x,getnlab(x),opt);
    end
  end
  out = trained_classifier(a,v);

elseif mapping_task(argin,'execution')

  [a,w] = deal(argin{1:2});
	v = getdata(w);
	n = size(a,1);  % nr objects
	K = size(w,2);  % nr of classes
	nrv = length(v); % nr of trees
    out = zeros(n,K);
    if exist('decisiontree')==3
      for j=1:nrv
        I = decisiontree(v{j},+a);
        out = out + accumarray([(1:n)' I],ones(n,1),[n K]);
      end
    else
      % the old fashioned slow Matlab code
      for i=1:n
        x = +a(i,:);
        for j=1:nrv
          I = tree_eval(v{j},x);
          out(i,I) = out(i,I)+1;
        end
      end
      out = out./repmat(sum(out,2),1,K);
    end
    out = setdat(a,out,w);

else
  error('Illegal call')
end

return

%    out = tree_eval(w,x)
%
function out = tree_eval(w,x)

n = size(x,1);
out = zeros(n,1);

for i=1:n

	v=w;
	% if the first split is already solving everything (1 obj. per class)
	if isa(v,'double')
		out(i,1) = v;
	end
	while (out(i,1)==0)
		if (x(i,v.bestf)<v.bestt)
			v = v.l;
		else
			v = v.r;
		end
		if isa(v,'double')
			out(i,1) = v;
		end
	end
end

%
%    w = tree_train(x,y,opt)
%
function w = tree_train(x,y,opt)

% how good are we in this node?
err = tree_gini(y,opt.K);
if (err==0)

	w = y(1); % just predict this label

else
	% we split further
	n = size(x,1);

	% optionally, choose only from a subset
	if (opt.featsubset>0)
		fss = randperm(size(x,2));
		fss = fss(1:opt.featsubset);
	else
		fss = 1:size(x,2);
	end

	% check each feature separately:
	besterr = inf; bestf = []; bestt = []; bestj = []; bestI = [];
	for i=fss
		% sort the data along feature i:
		[xi,I] = sort(x(:,i)); yi = y(I);
		% run over all possible splits:
		for j=1:n-1
			% compute the gini
			err = j*tree_gini(yi(1:j),opt.K) + (n-j)*tree_gini(yi(j+1:n),opt.K);
			% and see if it is better than before.
			if (err<besterr)
				besterr = err;
				bestf = i;
				bestj = j;
				bestt = mean(xi(j:j+1));
				bestI = I;
			end
		end
	end

	% store
	w.bestf = bestf;
	w.bestt = bestt;
	%  now find the children:
	w.l = tree_train(x(bestI(1:bestj),:),y(bestI(1:bestj)),opt);
	w.r = tree_train(x(bestI(bestj+1:end),:),y(bestI(bestj+1:end)),opt);
end

function g = tree_gini(y,K)

out = zeros(1,K);
for k=1:K
	out(k) = mean(y==k);
end

g = out*(1-out)';

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% The COBRAToolbox: testOptGene.m
%
% Purpose:
%     - test the optGene function
%
% Authors:
%     - Jacek Wachowiak

% Check requirements

requiredToolboxes = {'gads_toolbox'};  % This is the Global optimization toolbox
solvers = prepareTest('needsLP', true, 'needsMILP', true, 'toolboxes', requiredToolboxes);
% If we have more than one solver per type (LP/MILP), only use those that
% are available for both.
if numel(solvers.LP) > 1 && numel(solvers.MILP) > 1
    commonSolvers = intersect(solvers.LP, solvers.MILP);
    if ~isempty(commonSolvers)  % If there is any such solver
        solvers.MILP = commonSolvers;
        solvers.LP = commonSolvers;
    else  % Use only one solver. otherwise we can get into troubles.
        solvers.LP = solvers.LP(1);
        solvers.MILP = solvers.MILP(1);
    end
end
% save the current path
currentDir = pwd;

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

% test variables
model = getDistributedModel('ecoli_core_model.mat');
targetRxn = model.rxns{39};  % Succinate
fructose_substrateRxn = model.rxns{26};  % Fructose, even though this has no incluence whatsoever.
generxnList = model.rxns(setdiff([1:95], [11, 13, 26, 39]));  % Everything besides the ATP Maintenance, The biomass reaction and the substrate and target reactions.

klt=min([length(solvers.LP),length(solvers.MILP)]);
for k = 1:klt
    changeCobraSolver(solvers.LP{k}, 'LP', 0);
    changeCobraSolver(solvers.MILP{k}, 'MILP', 0);
    fprintf(' -- Running testOptGene using the solver interfaces: LP: %s ; MILP: %s... ', solvers.LP{k}, solvers.MILP{k});
    basicsolution = optimizeCbModel(model);
    % function outputs
    % requires Global Optimization Toolbox
    % Set the rng, for reproducability
    rng(0);
    [x, population, scores, optGeneSol] = optGene(model, targetRxn, fructose_substrateRxn, generxnList, 'StallTimeLimit', 15, 'TimeLimit', 30);
    % Check, that we get the expected solution from a previous run.
    optSols = population((optGeneSol.scores == min(optGeneSol.scores)), :);  % Get the set of optimal Solutions.
    optReacs = optSols(1, :);
    model2 = model;
    model2.lb(ismember(model2.rxns, generxnList(optReacs))) = 0;
    model2.ub(ismember(model2.rxns, generxnList(optReacs))) = 0;
    sol = optimizeCbModel(model2);
    % Lets only assert, that we have some improvement.
    assert(sol.v(39) - basicsolution.v(39) > 0);

end
% close the open windows
close all

% Remove the output, to keep the toolbox updateable.
delete([fileDir filesep 'MILPProblem.mat']);

fprintf('Done.\n');
% change to old directory
cd(currentDir);


text:
% this script represents the evolution of the covariance of an OU process in terms of the dispersion ellipsoid
% see A. Meucci (2009)
% "Review of Statistical Arbitrage, Cointegration, and Multivariate Ornstein-Uhlenbeck"
% available at ssrn.com

% Code by A. Meucci, April 2009
% Most recent version available at www.symmys.com > Teaching > MATLAB

clear; clc; close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input parameters of multivariate OU process
K=1;
J=1;

x0=rand(K+2*J,1);

Mu=rand(K+2*J,1);

A=rand(K+2*J,K+2*J)-.5;
ls=rand(K,1)-.5;
gs=rand(J,1)-.5;
os=rand(J,1)-.5;

S=rand(K+2*J,K+2*J)-.5;

ts=.01*[0:10:100];
NumSimul=10000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% process inputs
Gamma=diag(ls);
for j=1:J
    G=[gs(j)  os(j)
           -os(j) gs(j)];
    Gamma=blkdiag(Gamma,G);
end
Theta=A*Gamma*inv(A);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% one-step exact simulation
Sigma=S*S';
X_0=repmat(x0',NumSimul,1);
[X_t1, MuHat_t1, SigmaHat_t1]=OUstep(X_0,ts(end),Mu,Theta,Sigma);

% multi-step simulation: exact and Euler approximation
X_t=repmat(x0',NumSimul,1);
X_tE=X_t;
for s=1:length(ts)
    Dt=ts(1);
    if s>1
        Dt=ts(s)-ts(s-1);
    end
    [X_t,MuHat_t,SigmaHat_t]=OUstep(X_t,Dt,Mu,Theta,Sigma);
    %[X_tE,MuHat_tE,SigmaHat_tE]=OUstepEuler(X_tE,Dt,Mu,Theta,Sigma);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plots
Pick=[K+2*J-1 K+2*J];

% horizon simulations
hold on
h5=plot(X_t1(:,Pick(1)),X_t1(:,Pick(2)),'.');
set(h5,'color','r','markersize',4)

% horizon location
hold on
h4=plot(MuHat_t1(Pick(1)),MuHat_t1(Pick(2)),'.');
set(h4,'color','k','markersize',5)

% horizon dispersion ellipsoid
hold on
h3=TwoDimEllipsoid(MuHat_t1(Pick),SigmaHat_t1(Pick,Pick),2,0,0);
set(h3,'color','k','linewidth',2);

% starting point
hold on
h2=plot(x0(Pick(1)),x0(Pick(2)),'.');
set(h2,'color','b','markersize',5)

% starting generating dispersion ellipsoid
hold on
h1=TwoDimEllipsoid(x0(Pick),Sigma(Pick,Pick),2,0,0);
set(h1,'color','b','linewidth',2);

legend([h1 h3],'generator','horizon');


text:
% Copyright 2013 The MathWorks, Inc
%% Ground-Based Monostatic Radar - Range and Doppler Estimation
% This example shows how to model a ground-based monostatic pulse radar to
% estimate the range and speed of fluctuating targets.
%% Model
maxrange  = 8000; % Maximum range (m)
range_res = 50;   % Range resolution (m)
pd        = 0.9;  % Probability of detection
pfa       = 1e-6; % Probability of false alarm
nint      = 48;   % Number of pulses to integrate

% The transmitted waveform is a chirp. An S-band pyramidal antenna array is
% used in this case. The antenna is mounted on a vehicle moving at 67 mph
% (30 m/s).
apos = [0; 0; 0];                  % Antenna position
avel = [30*cos(30);30*sin(30);0];  % Antenna velocity
fc   = 3e9;                        % Operating frequency
sAnt = lowProfileArray('FrequencyRange',[2/3*fc 4/3*fc],'ViewArray',false);
[sWav,sTx,sAntPlat,sRad,fs,prf] = setupTx(maxrange,range_res,pd,pfa,nint,sAnt,apos,avel,fc);

% Three targets are set in motion in a free space environment.
tgtRCS = [1.2 1.1 1.05];
tgtpos = [2000 4560 5825; 0 0 0; 0 0 0];
tgtvel = [100 -400 350;0 0 0; 0 0 0];    % m/s
[sTgt,sTgtMotion,sChan] = setupTheater(tgtRCS,tgtpos,tgtvel,fc,fs);

% On the receiver side, matched filter, time-varying gain control and
% non-coherent pulse integration are applied to improve SNR. A
% range-Doppler map is generated and a Neyman-Pearson (NP) decision rule is
% used to achieve the desired Pfa.
nf          = 0;                                   % Noise figure
fast_time   = 0:1/fs:1/prf-1/fs;                   % Fast time grid
range_gates = physconst('LightSpeed')*fast_time/2; % Range gates
pulses      = zeros(numel(fast_time),nint);        % Pre-allocate
intpulses   = zeros(numel(fast_time),1);
[sCol,sRx,sRD,sMFilt,sTVG,threshold] = setupRx(nint,nf,pfa,maxrange,range_gates,sWav,sAnt,fc);

%% Simulate
rsig = zeros(336,3);
ang  = zeros(2,3);
for m = 1:2000
    [s,tx_status] = step(sTx,step(sWav));               % Transmit pulse
    [apos,avel]   = step(sAntPlat,1/prf);               % Move antenna
    for n = 1:3                                         % For each target
        [tpos,tvel]   = step(sTgtMotion{n},1/prf);      % Move target
        [~, ang(:,n)] = rangeangle(tpos,apos);          % Angle between antenna and target
        tsig          = step(sRad,s,ang(:,n));          % Radiate signal
        tsig          = step(sChan{n},tsig,apos,tpos,avel,tvel); % Propagate two ways
        rsig(:,n)     = step(sTgt{n},tsig,true);        % Reflect off target
    end
    rsig = step(sCol,rsig,ang);                         % Collect
    rsig = sum(rsig,2);                                 % Beamform
    nn   = mod(m-1,nint)+1;
    pulses(:,nn) = step(sRx,rsig,~(tx_status>0));       % Receiver pre-amp
    pulses(:,nn) = step(sTVG,pulses(:,nn));             % Time varying gain
    [rdmap,rgrid,sgrid] = step(sRD,pulses,sMFilt.Coefficients); % Range-Doppler estimate
    pulses(:,nn) = step(sMFilt,pulses(:,nn));           % Matched filter
    if nn == nint
        intpulses  = pulsint(pulses,'noncoherent');     % Pulse integration
        [pmax,detect] = findpeaks(intpulses,'MinPeakHeight',sqrt(threshold)); % Detection
        tgtrange   = range_gates(detect-(numel(sMFilt.Coefficients)-1));      % Range estimation
    end
    viewSignals
end