text:
function checkPropLevelsWarningsErrors(stateLog, celBodyData)
%checkPropLevelsWarningsErrors Summary of this function goes here
%   Detailed explanation goes here
    chunkedStateLog = breakStateLogIntoSoIChunks(stateLog);

    lowFuelOxEventsBodiesWarn = zeros(0,2);
    lowMonopropEventsBodiesWarn = zeros(0,2);
    lowXenonEventsBodiesWarn = zeros(0,2);

    lowFuelOxEventsBodiesAlert = zeros(0,2);
    lowMonopropEventsBodiesAlert = zeros(0,2);
    lowXenonEventsBodiesAlert = zeros(0,2);

    minWarnThresh = 0.1;
    minAlertThresh = 0.01;

    iniState = stateLog(1,:);
    iniFuelOx = iniState(10);
    iniMono = iniState(11);
    iniXenon = iniState(12);

    for(i=1:size(chunkedStateLog,1)) %#ok<*NO4LP>
        for(j=1:size(chunkedStateLog,2))
            subLog = chunkedStateLog{i,j};

            if(isempty(subLog) || (size(subLog,1)==0 || size(subLog,2)==0))
                continue;
            end

            state = subLog(end,:);

            bodyID = state(8);
            %bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);
            eventNum = state(13);

            fuelOx = state(10);
            monoprop = state(11);
            xenon = state(12);

            if(fuelOx / iniFuelOx <= minWarnThresh && fuelOx / iniFuelOx > minAlertThresh)
                lowFuelOxEventsBodiesWarn(end+1,:) = [eventNum, bodyID];
            elseif(fuelOx / iniFuelOx < minAlertThresh)
                lowFuelOxEventsBodiesAlert(end+1,:) = [eventNum, bodyID];
            end

            if(monoprop / iniMono <= minWarnThresh && monoprop / iniMono > minAlertThresh)
                lowMonopropEventsBodiesWarn(end+1,:) = [eventNum, bodyID];
            elseif(monoprop / iniMono < minAlertThresh)
                lowMonopropEventsBodiesAlert(end+1,:) = [eventNum, bodyID];
            end

            if(xenon / iniXenon <= minWarnThresh && xenon / iniXenon > minAlertThresh)
                lowXenonEventsBodiesWarn(end+1,:) = [eventNum, bodyID];
            elseif(xenon / iniXenon < minAlertThresh)
                lowXenonEventsBodiesAlert(end+1,:) = [eventNum, bodyID];
            end
        end
    end

    if(~isempty(lowFuelOxEventsBodiesAlert))
        eventNum = lowFuelOxEventsBodiesAlert(1,1);
        bodyID = lowFuelOxEventsBodiesAlert(1,2);
        addToExecutionErrors(['Low Fuel/Ox Alert (Events: ', makeEventsStr(lowFuelOxEventsBodiesAlert), ')'], -1, -1, celBodyData);
    end

    if(~isempty(lowFuelOxEventsBodiesWarn))
        eventNum = lowFuelOxEventsBodiesWarn(1,1);
        bodyID = lowFuelOxEventsBodiesWarn(1,2);
        C = setdiff(lowFuelOxEventsBodiesWarn,lowFuelOxEventsBodiesAlert,'rows');
        addToExecutionWarnings(['Low Fuel/Ox Warning (Events: ', makeEventsStr(C), ')'], -1, -1, celBodyData);
    end

    if(~isempty(lowMonopropEventsBodiesAlert))
        eventNum = lowMonopropEventsBodiesAlert(1,1);
        bodyID = lowMonopropEventsBodiesAlert(1,2);
        addToExecutionErrors(['Low Monoprop Alert (Events: ', makeEventsStr(lowMonopropEventsBodiesAlert), ')'], -1, -1, celBodyData);
    end

    if(~isempty(lowMonopropEventsBodiesWarn))
        eventNum = lowMonopropEventsBodiesWarn(1,1);
        bodyID = lowMonopropEventsBodiesWarn(1,2);
        C = setdiff(lowMonopropEventsBodiesWarn,lowMonopropEventsBodiesAlert,'rows');
        addToExecutionWarnings(['Low Monoprop Warning (Events: ', makeEventsStr(C), ')'], -1, -1, celBodyData);
    end

    if(~isempty(lowXenonEventsBodiesAlert))
        eventNum = lowXenonEventsBodiesAlert(1,1);
        bodyID = lowXenonEventsBodiesAlert(1,2);
        addToExecutionErrors(['Low Xenon Alert (Events: ', makeEventsStr(lowXenonEventsBodiesAlert), ')'], -1, -1, celBodyData);
    end

    if(~isempty(lowXenonEventsBodiesWarn))
        eventNum = lowXenonEventsBodiesWarn(1,1);
        bodyID = lowXenonEventsBodiesWarn(1,2);
        C = setdiff(lowXenonEventsBodiesWarn,lowXenonEventsBodiesAlert,'rows');
        addToExecutionWarnings(['Low Xenon Warning (Events: ', makeEventsStr(C), ')'], -1, -1, celBodyData);
    end
end

function str = makeEventsStr(eventsBodies)
    eventNums = unique(eventsBodies(:,1));
    str = '';
    for(i=1:length(eventNums))
        if(i==length(eventNums))
            endChar = '';
        else
            endChar = ', ';
        end

        str = [str, num2str(eventNums(i)),endChar];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function x = twos2dec(t)

% TWOS2DEC Convert binary string two's complement to decimal integer.
%
% Usage: X = TWOS2DEC(T)
%
% Converts the two's complement representation as a string given by T to
% decimal integers X. T can be a character array or cell string. Input
% multiple numbers into T along the rows, and the output will be as a
% column vector. Similarly to BIN2DEC, leading spaces in T are sign
% extended (so treated as either 0 or 1 depending on the first non space
% character), and embedded and trailing spaces are removed.
%
% Example:
%     >> twos2dec(['010 111'; ' 010111'; '101011 '; ' 10   1'])
%
%     ans =
%
%         23
%         23
%        -21
%         -3
%
% Inputs:
%   -T: string two's complement numbers to convert to decimal integers.
%
% Outputs:
%   -X: decimal representation of T.
%
% See also: DEC2TWOS, FIX2DEC, DEC2FIX, BIN2DEC, HEX2DEC, BASE2DEC.

error(nargchk(1, 1, nargin));
if iscellstr(t)
    t = char(t);
end

% Convert to numbers.
x = bin2dec(t);

% Get the number of bits as the number of 0's and 1's in each row.
nbits = sum(t == '0' | t == '1', 2);

% Case for negative numbers.
xneg = log2(x) >= nbits - 1;
% xneg = bitshift(x, -(nbits - 1)) == 1;
if any(xneg)
    x(xneg) = -( bitcmp(x(xneg), nbits(xneg)) + 1 );
end

author: Sable

repo: mcbench-benchmarks


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

%%
clear; close all; clc;

%% Syntax
% |varargout=patchEdges(varargin);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% Snake AI: Dynamic Hamiltonian Cycle Repair (with some strategic stuff)
% by Brian Haidet for AlphaPhoenix
% published 4/11/2020
% CC non-commercial, attribution
% See readme for details

%there are two methods that this file uses to reconstruct the hamiltonian
%cycle. first, it checks to see if the turn the snake WANTS to take splits
%the hamcycle into two closed loops. if so, it looks for a location where
%both loops have a flat edge to each other and splices the loops together.

%if that fails, there's a special method for a "loop of two" where a single
%line of two nodes was isolated and needs to be re-fused to tha path.

%if that fails, then there's a more expensive algorythm that attempts to
%draw a brand new hamiltionian cycle by expanding the current path into
%unclaimed regions. It's not vry smart, but it does ocasionally succeed
%(but normally leaves little one-off missing nodes everywhere and doesn't
%complete.)

%if none of these methods solve for a new hamiltionian path, the snake
%follows the old hamiltionian path.

%I'm very sure that somewhere in this file there are some missing minus
%signs or some flipped indeces to matrices or vectors because i sometimes
%see the snake fail to repair the path in a way that I believe it should
%know how, but as you can see, there are an awful lot of IFs and (-)s in
%this file, and I lacked the patience to find every glitch when working on
%this initially. Apparently whatever sometimes fails simply causes this
%file to return no hamiltonian path instead of a wrong or illegal
%hamiltonian path, and that was good enough at the time. if you find
%anything, god you must be bored, but if you do, let me know!!  :)

hypsnake=[nextstepinds,snake];
newcycle=[]; % populated if sucessful

% map cycle #1 (the one with the snake)
oldcutoffpointer_pathind=find(hampathinds==nextstepinds);
cycle1=hampathinds(oldcutoffpointer_pathind:end);
%map cycle #2 (without the snake)
cycle2=hampathinds(1:oldcutoffpointer_pathind-1);
if mod(length(cycle2),2)==1
    disp("it's broken")
end
%does cycle2 close? (AND isn't a 2-node line!)
if (sum(abs(coords(cycle2(1))-coords(cycle2(end))))==1) && length(cycle2)>2 %yes it closes
    cycle1fill=zeros(l);
    cycle1fill(cycle1)=1;
    hypfield=zeros(l);
    hypfield(hypsnake)=1;
    hashperim=10*(cycle1fill-hypfield)+~cycle1fill;
    [hitx,hity]=find(conv2(hashperim,[1,1;1,1])==22);
    for c=1:length(hitx)
        %find out if the adjacent parts of the paths are "flat" ie. connectable
        %for each (manually)
        %find candidate coords
        candidatesite=[hitx(c),hity(c)]+[-1,-1;0,-1;-1,0;0,0];
        candidatesiteinds=inds(candidatesite);
        candidatechains=[0,0,0,0];%candidates belong to cycle1 (0) or cycle2 (1)
        %where does each point?
        candidatepointers=[0,0,0,0];
        cyloc=[0;0;0;0];
        for p=1:4
            if cycle1fill(inds(candidatesite(p,:)))==1
                cyloc(p)=find(cycle1==candidatesiteinds(p,:));
                candidatepointers(p)=cycle1(mod(cyloc(p),length(cycle1))+1);
            else
                candidatechains(p)=1;
                cyloc(p)=find(cycle2==candidatesiteinds(p,:));
                candidatepointers(p)=cycle2(mod(cyloc(p),length(cycle2))+1);
            end
        end
        matchers=intersect(candidatesiteinds,candidatepointers);
        if length(matchers)==2%they match!

            if sum((matchers(1)==candidatesiteinds)&candidatechains') % if matchers(2) belongs to cycle 1
                cy2ind=sum((matchers(1)==candidatesiteinds).*cyloc);
                cy1ind=sum((matchers(2)==candidatesiteinds).*cyloc);
            else
                cy1ind=sum((matchers(1)==candidatesiteinds).*cyloc);
                cy2ind=sum((matchers(2)==candidatesiteinds).*cyloc);
            end
            newcycle=[cycle1(mod((1:cy1ind-1)-1,length(cycle1))+1);cycle2(mod(cy2ind-1:cy2ind+length(cycle2)-2,length(cycle2))+1);cycle1(cy1ind:end)];
            break;
        end
    end
else%no it does not close
    cycle1save=cycle1;
    cycle1=[repmat([-1],l,1);cycle1(1:end-length(snake));repmat([-1],l,1)];
    if length(cycle2) <=900 %try to fix the thing
        makingprogress=false;
        madeprogress=false;
        spliced=false;
        pi=1;%segment being spliced
        while true
            spliced=false;
            if abs(cycle2(pi)-cycle2(pi+1))==1 %splicing vertical segment
                pc1=find(cycle1==cycle2(pi)-l);%look to the left of pi
                if length(pc1)==1% if that spot IS on chain 1
                    if cycle1(pc1+1)==cycle2(pi+1)-l
                        newcycle=[cycle1(1:pc1) ; cycle2(pi) ; cycle2(pi+1) ; cycle1(pc1+1:end)];
                        madeprogress=true;
                        spliced=true;
                    elseif cycle1(pc1-1)==cycle2(pi+1)-l
                        newcycle=[cycle1(1:pc1-1) ; cycle2(pi+1) ; cycle2(pi) ; cycle1(pc1:end)];
                        madeprogress=true;
                        spliced=true;
                    end
                end
                if spliced==false %look to the right if the left isnt on chain 1
                    pc1=find(cycle1==cycle2(pi)+l);%look to the right of pi
                    if length(pc1)==1 % if that spot IS on chain 1
                        if cycle1(pc1+1)==cycle2(pi+1)+l
                            newcycle=[cycle1(1:pc1) ; cycle2(pi) ; cycle2(pi+1) ; cycle1(pc1+1:end)];
                            madeprogress=true;
                            spliced=true;
                        elseif cycle1(pc1-1)==cycle2(pi+1)+l
                            newcycle=[cycle1(1:pc1-1) ; cycle2(pi+1) ; cycle2(pi) ; cycle1(pc1:end)];
                            madeprogress=true;
                            spliced=true;
                        end
                    end
                end
            else %splicing horizontal segment
                if ~(mod(cycle2(pi),l)==1) %if on top row, don't look up
                    pc1=find(cycle1==cycle2(pi)-1);%look to the up of pi
                    if length(pc1)==1 % if that spot IS on chain 1
                        if cycle1(pc1+1)==cycle2(pi+1)-1
                            newcycle=[cycle1(1:pc1) ; cycle2(pi) ; cycle2(pi+1) ; cycle1(pc1+1:end)];
                            madeprogress=true;
                            spliced=true;
                        elseif cycle1(pc1-1)==cycle2(pi+1)-1
                            newcycle=[cycle1(1:pc1-1) ; cycle2(pi+1) ; cycle2(pi) ; cycle1(pc1:end)];
                            madeprogress=true;
                            spliced=true;
                        end
                    end
                end
                if spliced==false%look to the down if the up isnt on chain 1
                    if ~(mod(cycle2(pi),l)==0) %if on bottom row, don't look down
                        pc1=find(cycle1==cycle2(pi)+1);%look to the down of pi
                        if length(pc1)==1 % if that spot IS on chain 1
                            if cycle1(pc1+1)==cycle2(pi+1)+1
                                newcycle=[cycle1(1:pc1) ; cycle2(pi) ; cycle2(pi+1) ; cycle1(pc1+1:end)];
                                madeprogress=true;
                                spliced=true;
                            elseif cycle1(pc1-1)==cycle2(pi+1)+1
                                newcycle=[cycle1(1:pc1-1) ; cycle2(pi+1) ; cycle2(pi) ; cycle1(pc1:end)];
                                madeprogress=true;
                                spliced=true;
                            end
                        end
                    end
                end
            end
            if spliced==true
                cycle2(pi:pi+1)=[];
                cycle1=newcycle;
            else
                pi=pi+2;
            end
            if isempty(cycle2)
%                 disp('stage1')
                break;
            end
            if pi>length(cycle2)
                if madeprogress
                    pi=1;
                    madeprogress=false;
%                     %%plotting for debug
%                     figure(3)
%                     plot(snakecoords(:,2),snakecoords(:,1),'Color',[0 1 0],'LineWidth',15)
%                     xlim([0.5,l+.5]);
%                     ylim([0.5,l+.5]);
%                     axis square
%                     set(gca,'Ydir','reverse')
%                     set(gca,'Color','k')
%                     hold on
%                     plot(optimalPath(:,2),optimalPath(:,1),'Color',[1 0 0],'LineWidth',5)
%                     ncc2=coords(cycle1);
%                     plot(ncc2(:,2),ncc2(:,1),'Color',[0 .3 0],'LineWidth',3)
%                     hold off
                else
%                     disp('stage1 - FAIL')
                    break;
                end
            end

        end
        if ~isempty(cycle2)
            makingprogress=false;
            madeprogress=false;
            spliced=false;
            pi=l+1;%segment being spliced
            while true
                spliced=false;
                if abs(cycle1(pi)-cycle1(pi+1))==1 %splicing vertical segment
                    pc2a=find(cycle2==cycle1(pi)-l);%look to the left of pi
                    pc2b=find(cycle2==cycle1(pi+1)-l);%look to the left of pi+1
                    if (length(pc2a)==1)&&(length(pc2b)==1)% if those spot ARE on chain 2
                        newcycle=[cycle1(1:pi) ; cycle2(pc2a) ; cycle2(pc2b) ; cycle1(pi+1:end)];
                        madeprogress=true;
                        spliced=true;
                    end
                    if spliced==false %look to the right if the left isnt on chain 1
                        pc2a=find(cycle2==cycle1(pi)+l);%look to the right of pi
                        pc2b=find(cycle2==cycle1(pi+1)+l);%look to the right of pi+1
                        if (length(pc2a)==1)&&(length(pc2b)==1)% if those spot ARE on chain 2
                            newcycle=[cycle1(1:pi) ; cycle2(pc2a) ; cycle2(pc2b) ; cycle1(pi+1:end)];
                            madeprogress=true;
                            spliced=true;
                        end
                    end
                else %splicing horizontal segment
                    if ~(mod(cycle1(pi),l)==1) %if on top row, don't look up
                        pc2a=find(cycle2==cycle1(pi)-1);%look to the up of pi
                        pc2b=find(cycle2==cycle1(pi+1)-1);%look to the up of pi+1
                        if (length(pc2a)==1)&&(length(pc2b)==1)% if those spot ARE on chain 2
                            newcycle=[cycle1(1:pi) ; cycle2(pc2a) ; cycle2(pc2b) ; cycle1(pi+1:end)];
                            madeprogress=true;
                            spliced=true;
                        end
                    end
                    if spliced==false%look to the down if the up isnt on chain 1
                        if ~(mod(cycle1(pi),l)==0) %if on bottom row, don't look down
                            pc2a=find(cycle2==cycle1(pi)+1);%look to the down of pi
                            pc2b=find(cycle2==cycle1(pi+1)+1);%look to the down of pi+1
                            if (length(pc2a)==1)&&(length(pc2b)==1)% if those spot ARE on chain 2
                                newcycle=[cycle1(1:pi) ; cycle2(pc2a) ; cycle2(pc2b) ; cycle1(pi+1:end)];
                                madeprogress=true;
                                spliced=true;
                            end
                        end
                    end
                end
                if spliced==true
                    cycle2([pc2a,pc2b])=[];
                    cycle1=newcycle;
                else
                    pi=pi+1;
                end
                if isempty(cycle2)
%                     disp('stage2')
                    break;
                end
                if pi>length(cycle1)-2*l
                    if madeprogress
                        pi=l+1;
                        madeprogress=false;
%                         %%plotting for debug
%                         figure(3)
%                         plot(snakecoords(:,2),snakecoords(:,1),'Color',[0 1 0],'LineWidth',15)
%                         xlim([0.5,l+.5]);
%                         ylim([0.5,l+.5]);
%                         axis square
%                         set(gca,'Ydir','reverse')
%                         set(gca,'Color','k')
%                         hold on
%                         plot(optimalPath(:,2),optimalPath(:,1),'Color',[1 0 0],'LineWidth',5)
%                         ncc2=coords(cycle1);
%                         plot(ncc2(:,2),ncc2(:,1),'Color',[0 .3 .7],'LineWidth',3)
%                         hold off
                    else
%                         disp('stage2 - FAIL')
                        break;
                    end
                end

            end
        end
        if madeprogress
            newcycle=[newcycle(l+1:end-l); cycle1save(end-length(snake)+1:end)];
        else
            newcycle=[];
        end
    end
end

author: BrianHaidet

repo: AlphaPhoenix


text:
function Vout = disp_reaction(V,species)

% DISP_REACTION  Displays reactions for a given stoichiometric matrix.
%
% SYNTAX
%
% disp_reaction(V,species)
%
%   Displays the reactions corresponding given by the stoichiometric matrix
%   V. The columns of V correspond to reactions, the rows of V correspond
%   to species. species is a cell array of strings with labels for the
%   chemical species. Note that species is not parsed, so that labels can
%   need not be valid chemical formulas.
%
% Vout = disp_reaction(V,species)
% Vout = disp_reaction(V)
%
%   When called with an output argument, disp_reaction attempts to produce
%   integer coefficients for the stoichiometric matrix. The goal is to
%   express the coefficients efficiently. No reactions are displayed.
%

% AUTHOR
%
%    Jeff Kantor
%    December 19, 2010

    assert(nargin > 0, 'disp_reaction:input', ['No input. Expects a ', ...
                        'stoichiometric matrix.']);
    assert(nargin < 3, 'atomdisp_reactionic:input', 'Unexpected  inputs.');

    % Get size of the stoichiometric matrix. N is number of species, K is
    % number of reactions

    [N,K] = size(V);

    % Return if there are no species or no reactions to display.

    if ~isnumeric(V) || (N==0) || (K==0)
        disp('No reactions to display.');
        return
    end

    % Create generic species labels if no species are provided. Otherwise
    % check the list of species for obvious errors.

    if nargin < 2
        species = arrayfun(@(n)sprintf('Species_%d',n),1:N,'Uni',false);
    else
        assert(iscellstr(species), 'disp_reaction:input', ...
            'List of species must be a cell array.');
        assert(length(species(:)) == N, 'disp_reaction:input', ...
            'Number of species must equal rows of stoichiometry matrix.');
    end

    species = species(:);

    % The first step is to see if there is a better format for the
    % stoichiometric coefficients. We have 3 options for formatting the
    % stoichiometric coefficients: integer, rational, or floating point. We
    % form all three formats then choose a format for each reaction that
    % minimizes the display length.

    % Compute rational number approximation for the stoichiometric
    % coefficients.

    [num,den] = rat(V);

    % Find least common multiples of the denominators in each reaction. The
    % create output stoichiometric matrix with integer coefficients.

    Vout = zeros(N,K);
    for k = 1:K
        Vout(:,k) = num(:,k)*lcms(den(:,k))./den(:,k);
    end

    % Arrays to hold string representations of the stoichiometric
    % coefficients.

    Vs = cell(N,K);
    lens = zeros(1,K);

    Vf{N,K} = cell(N,K);
    lenf = zeros(1,K);

    Vr{N,K} = cell(N,K);
    lenr = zeros(1,K);

    % Smallest meaningful stoichiometric coefficient

    TOL = 1e-6;

    for k = 1:K
        for n = 1:N
            if abs(V(n,k)) >= TOL

                % Integer Coefficients

                Vs{n,k} = strtrim(sprintf('%d',abs(Vout(n,k))));

                if strcmp(Vs{n,k},'1')
                    Vs{n,k} = '';
                else
                    lens(k) = lens(k) + 1 + length(Vs{n,k});
                end

                % Rational Coefficients

                if (den(n,k) == 1) && (num(n,k) == 1);
                    Vr{n,k} = '';
                elseif (den(n,k) == 1)
                    Vr{n,k} = strtrim(sprintf('%d',abs(num(n,k))));
                    lenr(k) = lenr(k) + 1 + length(Vr{n,k});
                else
                    Vr{n,k} = strtrim(sprintf('%d/%d',abs(num(n,k)),den(n,k)));
                    lenr(k) = lenr(k) + 1 + length(Vr{n,k});
                end

                % Floating Point Coefficients

                Vf{n,k} = strtrim(sprintf('%10.5g',abs(V(n,k))));
                if strcmp(Vf{n,k},'1')
                    Vf{n,k} = '';
                else
                    lenf(k) = lenf(k) + 1 + length(Vf{n,k});
                end

            end
        end

        % Pick the best representation

        if (lenf(k) < lenr(k)) && (lenr(k) < lens(k))
            Vs(:,k) = Vf(:,k);
            lens(k) = lenf(k);
            Vout(:,k) = V(:,k);

        elseif (lenr(k) < lens(k))
            Vs(:,k) = Vr(:,k);
            lens(k) = lenr(k);
            Vout(:,k) = n(:,k)./den(:,k);

        end

    end

    % If there is no output then display the reactions

    if nargout < 1
        fprintf('\n');
        for k = 1:K

            pos = 1;

            % Print Reaction LHS -- Reactants

            n = find(V(:,k) < -TOL);
            pos = fprintrxn(strtrim(strcat(Vs(n,k),{' '},species(n))),pos);

            % Print separator

            if pos > 40
                fprintf('\n  <=> ');
                pos = 7;
            else
                fprintf(' <=> ');
                pos = pos + 5;
            end

            % Print Reaction RHS -- Products

            n = find(V(:,k) > TOL);
            fprintrxn(strtrim(strcat(Vs(n,k),{' '},species(n))),pos);

            fprintf('\n');

        end
        fprintf('\n');
    end

end

function q = lcms(v)

% LCMS  Find the least common multiple of a set of numbers.

    v = v(:);
    q = v(1);
    for k = 2:length(v)
        q = lcm(q,v(k));
    end
end

function pos = fprintrxn(terms,pos)

% FPRINTRXN  Helper function for displaying reactions.
%
% SYNTAX
%
% pos = fprintrxn(terms,pos)
%
%   Prints the reaction terms in the cell array terms separated by ' + '
%   starting at position pos. Breaks into multple lines if necessary.
%   Returns the position of the next character to be printed.

    % Print first term

    fprintf('%s',terms{1});
    pos = pos + length(terms{1});

    % Print remaining terms

    for m = 2:length(terms)

        % See where line would end

        tlen = 3 + length(terms{m});

        if (pos + tlen) <= 70

            fprintf(' + %s',terms{m});
            pos = pos + tlen;

        else % need a new line

            fprintf('\n   + %s',terms{m});
            pos = tlen + 2;

        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function asa266_test03 ( )

%*****************************************************************************80
%
%% TEST03 tests DIGAMMA, R8_PSI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 August 2010
%
%  Author:
%
%    John Burkardt
%
  ntest = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03\n' );
  fprintf ( 1, '  digamma(X) = d ( Log ( Gamma ( X ) ) ) / dX.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  DIGAMMA and\n' );
  fprintf ( 1, '  R8_PSI compute the digamma function:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X  DIGAMMA   R8_PSI\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest

    x = i / ntest;

    fprintf ( 1, '  %12f  %12f  %12f\n', x, digamma ( x ), r8_psi ( x ) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function scale = p39_scale ( neqn )

%*****************************************************************************80
%
%% P39_SCALE returns scale factors for problem p39.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real SCALE(NEQN), the scaling factors.
%
  scale(1:neqn,1) = 1.0;

  return
end


text:
function [cmd, outFile] = ctrScript(p,sName,outFile, outFileType)
%Create unix script to run ConTrack fiber generation
%
%   [cmd, outFile] = ctrScript(p,[sName],outFile, [outFileType])
%
% This saves the bash script used to run ConTrack.  The script is saved in
% the dt6 file directory.  The other files used by this script are located
% relative to this directory in the fibers\contrack directory.
%
% The usage for the Contrack fiber generation program is:
%
% USAGE:
%    ./contrack_gen.glxa64  -i <string> -p <string> [-t <float>] [-A] [-v]
%                           [--] [--version] [-h]
% Where:
%    -i <string>,  --info <string>
%      (required)  Information file, e.g. met_params.txt.
%    -p <string>,  --pdb <string>
%      (required)  Pathway database output file, e.g. paths.Bfloat.
%    -t <float>,  --time <float>
%      Pathway generation time limit.(EXPERIMENTAL--DON'T USE!!)
%    -A,  --Axial
%      Axial pathway generation.(EXPERIMENTAL--DON'T USE!!)
%    -v,  --vox
%      Voxelwise pathway generation.(EXPERIMENTAL--DON'T USE!!)
%    --,  --ignore_rest
%      Ignores the rest of the labeled arguments following this flag.
%    --version
%      Displays version information and exits.
%    -h,  --help
%      Displays usage information and exits.
%
%    ConTrack's pathway generation algorithm.
%
% Example:
%   sName = 'myScript.sh';
%   outFile = 'myTestFibers';
%   outFileType = 'pdb';
%   [cmd, outFile] = ctrScript(p,sName,outFile,outFileType)
%
% Author: AJS, BW
% 2008.09.17 DY & MP Edited the CMD to remove nohup and logFile to deal with
% errors running the command with Cygwin. These are not essential to
% running conTrack.
% 2008.12.04 DY&MP added the system command (chmod) to change the
% permissions of the .sh file (sName).
% 2009.03.10 GS & DY: add option to output pdb-format fiber files rather
% than Bfloat. added optional third input argument OUTFILETYPE, which
% defaults to Bfloat.
% 2009.06.11 DY: removed the '&' on the last line of the .sh file, which
% allows control over which processes are run in the background or not at
% the batch scripting level

if notDefined('p'), error('Parameters needed'); end

if notDefined('sName'),
    sName = ['ctrScript_',p.timeStamp,'.sh'];
    sName = fullfile(p.localSubjDir,'fibers','conTrack',sName);
end

if notDefined('outFileType')
    outFileType = '.pdb';  % should be .pdb
elseif outFileType(1) ~= '.'
    outFileType = ['.' outFileType];
end

if notDefined('outFile'),
    [tmp,roi1] = fileparts(p.roi1File);
    [tmp,roi2] = fileparts(p.roi2File);
    outFile = [roi1,'_',roi2,'_',p.timeStamp,outFileType];
end

% We create a relative path to the sampler file.  It is stored in the
% ..\fibers\contrack directory
[tmp,samplerFile,ext] = fileparts(p.samplerFile);
samplerFile = [samplerFile,ext];

% Get the log file text name.
[tmp,logFile,ext] = fileparts(p.logFile);
logFile = [logFile,ext];

% Create script for fiber generation
fid = fopen(sName,'wt');

% We write out an executiable bash script.  We write the script in the
% ..\fibers\contrack directory and we assume the user executes the script
% from inside the dt6 directory.  The output fiber file will be written to
% the ..\fibers\contrack directory.  The parameters file is also in the
% ..\fibers\contrack directory.
if ismac
    contrackBin = which('contrack_gen.maci64');
else % assume linux
    contrackBin = which('contrack_gen.glxa64');
end

fprintf(fid,'#!/bin/bash\n');
% cmd = sprintf('CMD=''nohup %s -i %s -p %s >%s;''',contrackBin,samplerFile,outFile,logFile);
cmd = sprintf('CMD=''%s -i %s -p %s''',contrackBin,samplerFile,outFile);
fprintf(fid,'%s \n',cmd);
fprintf(fid,'echo $CMD\n');
fprintf(fid,'$CMD\n');
fclose(fid);

% Edit permissions of the .sh file (sName) so that it can be executed.
[status,result] = system(['chmod 775 ' sName]);
if status ~= 0
    disp(['chmod failure in ctrScript Line 88: Permissions need to be edited manually for ' sName]);
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [divTr, divTe]= sample_KFold(label, folds, varargin)
%SAMPLE_KFOLD - Sampling function: random divisions (by default stratified)
%
%Synopsis:
%  [PARTR, PARTE]= sample_KFold(LABEL, FOLDS, <OPT>)
%
%Arguments:
% LABEL  - class label of size [nClasses x nSamples].
% FOLDS  - DOUBLE nFolds: number of folds into which the samples are
%          divided. Or FOLDS can be [nShifts nFolds] in which case
%          all partitions will also be generated in shifted versions.
% OPT    - property/value list of optinal parameters:
%   'stratified' [BOOL] stratified sampling (true, default)
%          or completely random sampling (false).
%
%Returns:
% DIVTR   - Partitions of the training set
%           DIVTR{n}: cell array holding the training sets folds for
%           shuffle #n, more specificially
%           DIVTR{n}{m} holds the indices of the training set of the m-th
%           fold of shuffle #n
% DIVTE   - analogue to DIVTR, for the test sets

% 2014-02 Martijn Schreuder

props = {'Stratified'      true          'BOOL|DOUBLE[1]'
        };

if nargin==0,
  divTr= props;
  return;
end

opt= opt_proplistToStruct(varargin{:});
[opt,~] = opt_setDefaults(opt, props, 1);

misc_checkType(label, 'DOUBLE[- -]');
misc_checkType(folds, 'DOUBLE|DOUBLE[2]');

nSamples = sum(label,2);

if length(folds)==1
  folds= [1 folds];
end

% check that the number of folds is smaller than the smallest class
if any(folds(2) > nSamples),
    error('The number of folds is larger than the number of samples in the smallest class');
end

%divTr= {cell(1,folds(2))};
%divTe= {cell(1,folds(2))};
for nn= 1:folds(1)
  clear idx;
  % prepare indices
  if opt.Stratified
      for cl = 1:length(nSamples)
          clid = find(label(cl,:));
          idx{cl} = clid(randperm(nSamples(cl)));
      end
  else
      idx{1} = randperm(sum(nSamples));
  end

  % make divisions
  for cl = 1:length(idx)
      div{cl}= round(linspace(0, nSamples(cl), folds(2)+1));
  end

  % sample
  for kk= 1:folds(2)
    divTe{nn}{kk} = [];
    for cl = 1:length(div),
      divTe{nn}{kk}= sort([divTe{nn}{kk} idx{cl}(div{cl}(kk)+1:div{cl}(kk+1))]);
    end
    divTr{nn}{kk}= setdiff(1:sum(nSamples), divTe{nn}{kk});

    % check that all classes are inhabited
    if ~all(sum(label(:,divTr{nn}{kk}))),
      error('empty classes in training set');
    end
  end
end

author: bbci


text:
%--- help for polyshape/translate ---
%
%  TRANSLATE Translate a polyshape
%
%  PG = TRANSLATE(pshape, V) translates a polyshape according to a
%  two-element row vector V. The first element of V is the translation
%  distance in the x direction, and the second element is the translation
%  distance in the y direction. Positive values in V translate right and up,
%  and negative values translate left and down. When pshape is an array of
%  polyshapes, each element of pshape is translated according to V.
%
%  PG = TRANSLATE(pshape, x, y) specifies the x and y translation distances
%  as separate arguments.
%
%  See also scale, rotate, polybuffer, polyshape
%
%    Reference page in Doc Center
%       doc polyshape/translate
%
%    Other functions named translate
%
%       dsge/translate
%

author: jmaih


text:
function inspect_menu_viewpoint

% Some of the FT_PLOT_XXX functions that return a 3D object support a
% right-mouse-click context menu with which you can select
% top/bottom/left/right/front/back. This functionality requires that the object being
% plotted has a known coordinate system.

% note that the objects don't really go together, so it does not fit 100%

figure

%%

elec = ft_read_sens('GSN-HydroCel-128.sfp');
elec = ft_convert_units(elec, 'mm');
elec.coordsys = 'ras';

ft_plot_sens(elec, 'label', 'label'); axis off

%%

mesh = ft_read_headshape('cortex_8196.surf.gii');
mesh.coordsys = 'mni';

ft_plot_mesh(mesh, 'facealpha', 0.5, 'unit', 'mm')
lighting gouraud

%%
% dipole pointing to anterior and superior

ft_plot_dipole([40 0 50], [0 1 1], 'coordsys', 'ras', 'unit', 'mm', 'axes', true)

%%

ft_plot_topo3d(elec.elecpos, elec.elecpos(:,3), 'coordsys', 'ras', 'unit', 'mm', 'axes', true)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [trend,cyclic] = beveridgenelson(y,parameters,constant,p,q)
% Beveridge-Nelson decomposition for a trending time series.
%
% USAGE:
%   [TREND,CYCLIC] = beveridgenelson(y,parameters,constant,p,q)
%
% INPUTS:
%   Y          - T by 1 column of trending data data
%   PARAMETERS - 1 by CONSTANT + length(P) + length(Q) vector of parameters estimated on diff(Y)
%   CONSTANT   - Scalar variable: 1 the the fit model included a constant
%   P          - Non-negative integer vector representing the AR orders of the fit model
%   Q          - Non-negative integer vector representing the MA orders of the fit model
%
% OUTPUTS:
%   TREND      - T by 1 vector of containing the trend.  First max(P)+1 observations are the same as Y
%   CYCLIC     - T by 1 vector of containing the cyclic component.  First max(P)+1 observations are
%                  the same as Y
%
% COMMENTS:
%   Uses the exact BN decomposition as presented in
%
%   Paul Newbold, "Precise and efficient computation of the Beveridge-Nelson decomposition of
%     economic time series", Journal of Monetary Economics, Volume 26, Issue 3, December 1990, Pages
%     453-457.
%
%   This decomposition is based on a demeaned ARMA model for the short-run dynamics.  If CONSTANT = 1,
%   this function will demeand the data using the model implied long-run mean (the constant
%   parameter divided by the sum of the AR coefficients).
%
% EXAMPLES:
%   BN decomposition for Real US GDP using an AR(1) for the SR component
%       load GDP
%       parameters = armaxfilter(diff(lnGDP),1,1);
%       [trend,cyclic] = beveridgenelson(lnGDP,parameters,1,1);
%   BN decomposition for Real US GDP using an ARMA(1,4) for the SR component
%       parameters = armaxfilter(diff(lnGDP),1,1,1:4);
%       [trend,cyclic] = beveridgenelson(lnGDP,parameters,1,1,1:4);
%
% See also ARMAXFILTER

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 11/1/2009

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 3
        p=[];
        q=[];
    case 4
        q=[];
    case 5
        % Nothing
    otherwise
        error('3 to 5 inputs required.');
end
if size(y,2)>size(y,1)
    y=y';
end
if size(y,2)~=1
    error('Y must be a T by 1 vector.');
end
if ~ismember(constant,[0 1])
    error('CONSTANT must be 0 or 1.')
end
np = length(p);
nq = length(q);
if np>0
    maxp = max(p);
else
    maxp = 0;
end
if nq>0
    maxq = max(q);
else
    maxq = 0;
end
if length(parameters)~=(constant+np+nq)
    error('PARAMETERS must have CONSTANT + length(P) + length(Q) parameters.');
end
% Remove the constant, if any
if constant
    constant = parameters(1);
    parameters = parameters(2:length(parameters));
    if np>0
        longRunMean = constant/(1-sum(parameters(1:np)));
    else
        longRunMean = constant;
    end
end
deltaY = diff(y) - longRunMean;
if maxp>0
    A = zeros(maxp);
    A(1,p) = parameters(1:np);
    for i=2:maxp
        A(i,i-1) = 1;
    end
    Aeigs = abs(eig(A));
    if  max(Aeigs)>=1
        error('The model for the short-run dynamics is not stationary. Stationarity is required for the BN decomposition to be well defined.');
    end
else
    A = 0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

T = size(deltaY,1);
% Construct the errors
m = max(maxp,maxq);
zerosToPad = max(maxq-maxp,0);
deltaYAugmentedForMA = [zeros(zerosToPad,1);deltaY];
errors = armaxerrors(parameters,p,q,0,deltaYAugmentedForMA,[],m,ones(size(deltaYAugmentedForMA)));
errors = errors(zerosToPad+1:length(errors));

% Initialize the trend and cyclic component
trend = zeros(T+1,1);
trend(1:maxp+1) = y(1:maxp+1);
cyclic = zeros(T+1,1);
deltaYForecast = zeros(T+m,1);
errorsForecast = zeros(1,T+m);
% Initialize the coefficient needed
e = [1 zeros(1,maxp-1)];
longRunScale = (eye(maxp) - A)^(-1)*A;
longRunScale = e*longRunScale;
for t = maxp + 1 : T
    % Construct maxq forecasts, then use companion matrix to compute sum
    deltaYForecast(:) = 0;
    errorsForecast(:) = 0;
    deltaYForecast(1:t) = deltaY(1:t);
    errorsForecast(1:t) = errors(1:t);
    for h = 1 : maxq
        for j=1:np
            deltaYForecast(t+h) = deltaYForecast(t+h) + parameters(j) * deltaYForecast(t+h-p(j));
        end
        for j=1:nq
            if (t+h-q(j))>0
                deltaYForecast(t+h) = deltaYForecast(t+h) + parameters(np+j) * errorsForecast(t+h-q(j));
            end
        end
    end
    % Select the correct forecasts
    cumForecast = sum(deltaYForecast(t+1:t+maxq));
    if np>0
        deltaYHat=deltaYForecast(t+maxq-maxp+1:t+maxq);
        cumForecast = cumForecast + longRunScale*deltaYHat;
    end
    trend(t+1) = y(t+1) + cumForecast;
    cyclic(t+1) = -cumForecast;
end

author: bashtage

repo: mfe-toolbox


text:
function val=convertSpeedUnits(val,unitOrigNum,unitOrigDenom,unitDestNum,unitDestDenom)
%%CONVERTSPEEDUNITS Convert values of speed from one set of units to
%                   another. Speed is given as a ratio of a distance unit
%                   to a time unit. Note that 1 knot (kt)=1 nml/h. Thus, a
%                   conversion from knots to meters per second would have
%                   unitOrigNum='nml';unitOrigDenom='h';unitDestNum='m';
%                   unitDestDenom='s';
%
%INPUTS: val The matrix or vector of values that are to be converted.
% unitOrigNum,unitOrigDenom
% unitDestNum,unitDestDenom Four character strings indicating the units of
%            val and the units into which it is to be converted.
%            unitOrigNum and unitDestNum are the units of distance and can
%            take the values listed in convertTimeUnits. unitOrigDenom and
%            unitDestDenom are the length units and can take the values
%            listed in convertLengthUnits.
%
%OUTPUTS: val The values converted into the desired coordinate system.
%
%This just uses the functions convertLengthUnits and convertTimeUnits to
%get the approprimate multiplication factors to convert the length and time
%components.
%
%May 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

lengthCoeff=convertLengthUnits(1,unitOrigNum,unitDestNum);
%Inverse, because time is in the denominator.
timeCoeff=1/convertTimeUnits(1,unitOrigDenom,unitDestDenom);

val=val*lengthCoeff*timeCoeff;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function env_load_dependencies(dependency_dir,autocompile)
% Execute dependency loaders in the given directory tree.
% env_load_dependencies(Directory)
%
% This function recognizes 3 types of dependency loader scripts in the dependencies (sub-)folders.
%  1. Scripts called env_exec.m. Each such script is executed by env_load_dependencies.
%
%  2. Scripts called env_add.m. Each such function is executed, and if the 'ans' variable is true
%     after executing the script (it is set to true before running the script), and no exception
%     had occured, the containing directory is added using addpath.
%
%  3. Scripts called env_compile.m. The presence of such a file indicates to the loader function
%     that any mex or Java source files in the given directory should be automatically compiled.
%     Additional compiler inputs (such as defines, libraries and include directories), if any, may
%     be declared as Name=value; statements in the file itself; These will be passed as named arguments
%     to hlp_trycompile (examples can be found throughout the existing dependencies).
%     A special parameter is Skip, which can be set to false on platforms on which the compilation
%     shall be skipped (perhaps because it is known to fail).
%
% In:
%   Directory : directory tree with dependencies to be loaded; by default, searches for a directory
%               named 'dependencies', relative to the current path or the path of this file.
%
%   Compile : whether to act on env_compile.m files (default: true)
%
% Example:
%   % load the dependencies in the C:\mydependencies folder.
%   env_load_dependencies('C:\mydependencies')
%
% Notes:
%   If you intend to compile your toolbox using the MATLAB compiler, you cannot have empty .m files.
%   Therefore, your env_add etc. files must contain at least a space character to satisfy the compiler.
%
% See also:
%   env_startup
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-03-06

if ~exist('dependency_dir','var')
    % other possible directories where the dependencies could be located:
    possible_dirs = {pwd,fileparts(mfilename('fullpath'))};
    while ~exist('dependency_dir','var')
        % for each of the possible directories...
        for d=1:length(possible_dirs)
            % check whether it contains a 'dependencies' directory
            curdir = possible_dirs{d};
            if exist([curdir filesep 'dependencies'],'dir')
                dependency_dir = [curdir filesep 'dependencies'];
                break;
            end
            % otherwise peel off the last directory
            fsp = find(curdir==filesep,1,'last');
            if ~isempty(fsp)
                possible_dirs{d} = curdir(1:fsp-1); end
        end
        if cellfun('isempty',possible_dirs)
            error('Did not find a directory called ''dependencies''. Please specify it.'); end
    end
else
    if ~ischar(dependency_dir)
        error('The given dependency directory must be a string.'); end
    if ~exist(dependency_dir,'dir')
        error('The given dependency directory (%s) does not exist.',dependency_dir); end
end

if ~exist('autocompile','var')
    autocompile = true; end

if autocompile
    % check if compilation is suppressed on this platform
    try
        fname = [hlp_homedir filesep '.bcilab' filesep 'bcilab-dont-compile'];
        fid = fopen(fname,'r');
        platforms = fgetl(fid);
        fclose(fid);
        if strfind(platforms,[hlp_hostname ';'])
            autocompile = false; end
    catch
    end
end

load_dependencies(genpath(dependency_dir),autocompile);

% now, make sure that the build-* directories are included after their ancestor directories...
if ~isdeployed
    % get current path order
    added_paths = strsplit(path,pathsep)';
    added_paths = sort(added_paths); %#ok<TRSRT>
    added_paths = added_paths(cellfun(@(s)strncmp(dependency_dir,s,length(dependency_dir)),added_paths));

    super_dirs = {}; % this is a cell array of super directories with build dirs that may have to be moved
    sub_dirs = [];   % this is an array of indices (into added_paths) of sub-directories that need
                     % to be moved in front of the respective super dirs (indexed the same)
    % we start from the end so that we don't have to change indices
    for p=length(added_paths):-1:1
        try
            cur = added_paths{p};
            k = strfind(cur,[filesep 'build-']);
            if ~isempty(k)
                % if this is a build directory, add to the potential-move list
                super_dirs{end+1} = cur(1:k(end)-1); %#ok<AGROW>
                sub_dirs(end+1) = p; %#ok<AGROW>
            end
            % check if this is one of the super dirs
            match = strcmp(cur,super_dirs);
            if any(match)
                % move all the respective sub-dirs to right before this position
                indices = sub_dirs(match);
                added_paths = vertcat(added_paths(1:p-1),added_paths(indices),added_paths(setdiff(p:end,indices)));
                super_dirs = super_dirs(~match);
                sub_dirs = sub_dirs(~match);
            end
        catch
            warning('bcilab:dependency_glitch','Glitch during mex/build path reordering.');
            super_dirs = {};
            sub_dirs = [];
        end
    end

    % rebuild path
    if ~isempty(added_paths)
        rmpath(added_paths{:});
        addpath(added_paths{:});
    end
end

% from the given path list, add all paths that contain an active env_add.m
% file (active: either empty or runs without exceptions and ans is not being
% set to false), and execute all env_exec.m files
function load_dependencies(pathlist,autocompile)
compile_problems = {};
paths = sort(strsplit(pathlist,pathsep));
for px = paths(end:-1:1)
    p = px{1};
    if exist([p filesep 'env_add.m'],'file')
        % env_add file: check if "active"
        [dummy,isactive] = evalc('runscript([p filesep ''env_add.m''])');  %#ok<ASGLU>
        % immediately add all java .jar packages
        dirinfos = dir([p filesep '*.jar']);
        for m={dirinfos.name}
            if isdeployed
                warning off MATLAB:javaclasspath:jarAlreadySpecified; end
            evalc('javaaddpath([p filesep m{1}]);');
        end
    else
        isactive = false;
    end
    % env_exec file: run it
    if exist([p filesep 'env_exec.m'],'file')
        runscript([p filesep 'env_exec.m']); end
    % env_compile file: invoke compilation, if necessary
    if autocompile && exist([p filesep 'env_compile.m'],'file')
        settings = get_settings([p filesep 'env_compile.m']);
        % the Skip field determines whether the current platform should be skipped
        if isfield(settings,'Skip')
            docompile = ~settings.Skip;
        else
            docompile = true;
        end
        if ~isequal(settings,false) && docompile
            if ~hlp_trycompile('Directory',p, 'Style','eager', settings);
                compile_problems{end+1} = p; end %#ok<AGROW>
        end
    end
    if isactive
        % finally add to MATLAB & Java path
        if ~isdeployed
            addpath(p); end
        % also add java files here (note: this may clear global variables)
        if ~isempty(dir([p filesep '*.class']))
            if isdeployed
                warning off MATLAB:javaclasspath:jarAlreadySpecified; end
            javaaddpath(p);
        end
    end
end

% handle compilation errors gracefully
if ~isempty(compile_problems) && ~isdeployed
    disp('There were compiling problems for the dependencies in:');
    for k=1:length(compile_problems)
        disp(['  ' compile_problems{k}]); end
    fprintf('\n');
    disp('This means that the affected features of BCILAB will be disabled.');
    if ispc
        disp('These issues are usually resolved by installing a (free) Microsoft Visual Studio Express version that is between five years to a half year older than your MATLAB release.');
        disp('If you install a new compiler, you might have to type "mex -setup" in the MATLAB command line, before MATLAB recognizes it.');
    elseif ismac
        disp('These issues are usually resolved by installing the Xcode application that comes with your OS (with all install options to include gcc compiler support).');
    elseif isunix
        disp('These issues are usually resolved by a installing a version of the gcc compiler package that is between three years to a half year older than your MATLAB release.');
    end
    yn = lower(input('Should BCILAB try to recompile in the future (or skip recompilation if not)? [y/n]','s'));
    if any(strcmp(yn,{'y','n','yes','no'}))
        if yn(1) == 'n'
            if ~exist([hlp_homedir filesep '.bcilab'],'dir')
                if ~mkdir(hlp_homedir,'.bcilab');
                    disp('Cannot create directory .bcilab in your home folder.'); end
            end
            fname = [hlp_homedir filesep '.bcilab' filesep 'bcilab-dont-compile'];
            disp(['BCILAB will add/update a file called ' fname ' which contains the names of computers on which you chose to disable compilation.']);
            fid = fopen(fname,'a');
            if fid ~= -1
                try
                    fwrite(fid,[hlp_hostname ';']);
                    fclose(fid);
                catch
                    try fclose(fid); catch,end
                    disp('There were permission problems while accessing this file. Did not update it.');
                end
            else
                disp('The file name appears to be invalid on your platform; did not update it.');
            end
        end
    else
        disp('Please type ''y'' or ''n'' next time.');
    end
    fprintf('\n');
end

% run the given script (in its own scope...)
% and return the result of the last line
function ans = runscript(filename__)
try
    ans = 1; %#ok<NOANS>
    fprintf('Loading %s...\n',fileparts(filename__));
    run_script(filename__); % Note: This function is equivalent to run(), but works in deployed mode,
                            %       too. If the code fails here, you probably don't have this
                            %       function in your path. You may retrieve it (it's in the 'misc'
                            %       directory of BCILAB), or replace it here by run().
catch e
    disp(['Error running dependency loader ' filename__ '; reason: ' e.message]);
    disp('The dependency will likely not be fully operational.');
    ans = 0; %#ok<NOANS>
end

% run the given settings script (in its own scope...)
% and return the collected defined variables (or false)
function settings = get_settings(filename__)
try
    evalc('run_script(filename__)');
    infos = whos();
    settings = struct();
    for n = {infos.name}
        settings.(n{1}) = eval(n{1}); end
catch
	disp(['The settings script ' filename__ ' contains an error.']);
    settings = false;
end

% Split a string according to some delimiter(s). Not as fast as hlp_split (and doesn't fuse
% delimiters), but doesn't need bsxfun().
function strings = strsplit(string, splitter)
ix = strfind(string, splitter);
strings = cell(1,numel(ix)+1);
ix = [0 ix numel(string)+1];
for k = 2 : numel(ix)
    strings{k-1} = string(ix(k-1)+1:ix(k)-1); end
strings = strings(~cellfun('isempty',strings));

repo: REST


text:
function bit = i4_bit_hi1 ( n )

%*****************************************************************************80
%
%% I4_BIT_HI1 returns the position of the high 1 bit base 2 in an integer.
%
%  Example:
%
%       N    Binary     BIT
%    ----    --------  ----
%       0           0     0
%       1           1     1
%       2          10     2
%       3          11     2
%       4         100     3
%       5         101     3
%       6         110     3
%       7         111     3
%       8        1000     4
%       9        1001     4
%      10        1010     4
%      11        1011     4
%      12        1100     4
%      13        1101     4
%      14        1110     4
%      15        1111     4
%      16       10000     5
%      17       10001     5
%    1023  1111111111    10
%    1024 10000000000    11
%    1025 10000000001    11
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the integer to be measured.
%    N should be nonnegative.  If N is nonpositive, the value will always be 0.
%
%    Output, integer BIT, the number of bits base 2.
%
  i = floor ( n );
  bit = 0;

  while ( 1 )

    if ( i <= 0 )
      break;
    end

    bit = bit + 1;
    i = floor ( i / 2 );

  end

  return
end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [ data ] = readOffsetRep( ndata, labelset )
%READOFFSETREP read the offset representations and  recover the boxes with
%absolute positions
%
% INPUT: ndata - the ndata is a struct with
%                 "kids" (hierarchy)
%                 "leafreps" (layer+label of boxes)
%                 "relposreps" (relative positions between child
%                 and parent)
%        labelset - the category list for this dataset
%
% OUTPUT: data - the data is a struct with
%                 "kids" (hierarchy) the same with input
%                 "obblist" (boxes with absolute positions)
%                 "labellist" (labels of the boxes)
%                 "layers" (layers of the boxes)

kids = ndata.kids;
leafreps = ndata.leafreps;
relposreps = ndata.relposreps;
leafnum = size(leafreps,2);
nodenum = length(kids);

%% extract the labels
classnum = length(labelset);
labelvec = leafreps(end-classnum+1:end,:);
labellist = {};
for i = 1:leafnum
    label = labelvec(:,i);
    [~,I] = max(label);
    label = labelset{I};
    labellist{i} = label;
end

%% extract the layers
if(size(leafreps,1)-classnum>0)
    layers = leafreps(1,:);
end

%% compute the absolute positions of the boxes
obblist = zeros(12,nodenum);
% obblist(:,end) = [0.0322;0;-0.0885;1;0;0;0;1;0;0.49;0;1.0670]; % the root box
obblist(:,end) = [0;0;0;1;0;0;0;1;0;1;0;1];
% [ psize ] = getRootSize( ndata );
% obblist(10:12,end) = psize;
lastmergerep = relposreps{end};
rootrp = lastmergerep(:,1);
obblist(10,end) = rootrp(2);
obblist(12,end) = rootrp(3);
for i = nodenum:-1:1
    pobb = obblist(:,i);
    k = kids{i};
    rplist = relposreps{i};
    if(length(k)>1)
        for j = 2:length(k)
            rp = rplist(:,j-1);
            cobb = genOBBfrom2dOffset( pobb, rp );
            obblist(:,k(j)) = cobb;
        end
    end
end

% form the data
data.kids = kids;
data.labellist = labellist;
data.obblist = obblist;%(:,1:leafnum);
if(exist('layers','var'))
    data.layers = layers;
end

end

repo: GRAINS


text:
 function X = dtft1(x, omega, n_shift)
%function X = dtft1(x, omega, n_shift)
%| Compute DTFT of 1D signal x at frequency locations wx
%| In
%|	x	[N,ncol]	signal values
%|	omega	[M,1]		frequency locations
%|	n_shift	[1,1]		use [0:(N-1)]-n_shift indices
%| Out
%|	X	[M,ncol]	DTFT values
%|
%| This needs enough memory to store M * N size matrices
%| Copyright 2000-1-9, Jeff Fessler, University of Michigan

if nargin == 1 && streq(x, 'test'), dtft1_test, return, end
if nargin < 2, ir_usage(), end

N = size(x,1);
if ~isvar('n_shift') || isempty(n_shift), n_shift = 0; end
nn = [0:(N-1)] - n_shift;

X = exp(-1i*omega*nn) * x; % compute 1D DTFT

function dtft1_test()
N = 4;
x = [[1:N]', [1 1 2 2]'];	% two test signal
omega = 2*pi*[0:(N-1)]'/N;	% test with uniform frequency locations
Xd = dtft1(x, omega);
Xf = fft(x);
%disp([Xd Xf])
printm('max %% difference = %g', max_percent_diff(Xf,Xd))

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function X_rec = recoverData(Z, U, K)
%RECOVERDATA Recovers an approximation of the original data when using the
%projected data
%   X_rec = RECOVERDATA(Z, U, K) recovers an approximation the
%   original data that has been reduced to K dimensions. It returns the
%   approximate reconstruction in X_rec.
%

% You need to return the following variables correctly.
X_rec = zeros(size(Z, 1), size(U, 1));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the approximation of the data by projecting back
%               onto the original space using the top K eigenvectors in U.
%
%               For the i-th example Z(i,:), the (approximate)
%               recovered data for dimension j is given as follows:
%                    v = Z(i, :)';
%                    recovered_j = v' * U(j, 1:K)';
%
%               Notice that U(j, 1:K) is a row vector.
%

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

end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function Test_quadraticFunction

% Define a linear space
m = 5;
n = 2;
manifold = euclideanfactory(m,n);

% random symmetric square matrix
Q = multisym(randn(m));

% Create the problem structure.
problem.M = manifold;

% Define the problem cost function, Euclidean gradient and Hessian.
problem.cost = @cost;
function f = cost(X)
  f = 0.5*trace(X'*Q*X);
end

problem.egrad = @egrad;
function G = egrad(X)
  G = Q*X;
end

problem.ehess = @ehess;
function H = ehess(X, Xdot)
  H = Q*Xdot;
end

% Check gradient and Hessian correctness
% This is a quadratic function in a linear manifold,
% so the quadratic approximation used in `checkhessian` should be exact
% Other similar cases are those where the unknowns belong to a tangent
% space (also a linear manifold).
checkgradient( problem ); pause
checkhessian(  problem ); pause

end

author: NicolasBoumal

repo: manopt


text:
%% RATE OF CONVERGENCE OF ADAPTIVE FINITE ELEMENT METHOD USING WG ELEMENT
%
% This example is to show the rate of convergence of the lowest order weak
% Galerkin element approximation of the second order elliptic equation.
%
% Reference
%
% L. Chen, J. Wang and X. Ye. A posteriori error estimates for Weak
% Galerkin finite element methods for second order elliptic problems.
% Journal of Scientific Computing, 59(2), 496-511, 2014.

%% Lshape problem
% mesh
[node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
[node,elem] = delmesh(node,elem,'x>0 & y<0');
bdFlag = setboundary3(node,elem,'Dirichlet');
mesh = struct('node',node,'elem',elem,'bdFlag',bdFlag);
% pde
pde = Lshapedata3;
% option
format shorte
option.L0 = 1;
option.maxIt = 50;
option.printlevel = 1;
option.elemType = 'WG';
option.plotflag = 1;
option.maxN = 2e4;
% AFEM
err = afemPoisson3(mesh,pde,option);
% plot
figure;
showrate2(err.N,err.H1,10,'k-*','|| Du-Du_h||',err.N,err.eta,10,'k-+','eta');
% latexerrtable(err.N,[err.H1 err.eta])

%% Jump coefficients
% The diffusion coefficent is piecewise constant with large jump.
[node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
bdFlag = setboundary3(node,elem,'Dirichlet','(x==1) | (x==-1)');
mesh = struct('node',node,'elem',elem,'bdFlag',bdFlag);
% pde
pde = jumpmgdata1;
global epsilon
epsilon = 1e-4;
% option
option.L0 = 3;
option.maxIt = 50;
option.printlevel = 1;
option.elemType = 'WG';
option.plotflag = 1;
option.maxN = 1e4;
option.viewcut = '~(x<=0 & y>=0 & z>=0)';
option.viewangle = [59,20];
% AFEM
err = afemPoisson3(mesh,pde,option);
% plot
figure;
showrate2(err.N,err.H1,10,'k-*','||Du-Du_h||',err.N,err.eta,10,'k-+','eta');
% Note no exact solution to compute the error. Code the energy later on.

repo: ifem


text:
classdef Vector3D < AbstractTensor ...
                    & FirstOrderDescriptor ...
                    & Elasticity3dDescriptor

    methods (Access = public)

        function normalize(obj)
            d = obj.getValue();
            d = d/norm(d);
            obj.setValue(d);
        end

    end

    methods (Access = protected)

        function loadOrderVariable(obj)
            obj.order = 'first';
        end

        function loadTensorSize(obj)
            obj.tensorSize = [3 1];
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function scores = calc_scores(input_dir,GT_dir,shave_width,verbose)

addpath(genpath(fullfile(pwd,'utils')));

%% Loading model
load modelparameters.mat
blocksizerow    = 96;
blocksizecol    = 96;
blockrowoverlap = 0;
blockcoloverlap = 0;

%% Reading file list
file_list = dir([input_dir,'/*.png']);
im_num = length(file_list);

%% Calculating scores
scores = struct([]);

for ii=1:im_num
    if verbose
        fprintf(['\nCalculating scores for image ',num2str(ii),' / ',num2str(im_num)]);
    end

    % Reading and converting images
    input_image_path = fullfile(input_dir,file_list(ii).name);
    input_image = convert_shave_image(imread(input_image_path),shave_width);
    GD_image_path = fullfile(GT_dir,file_list(ii).name);
    GD_image = convert_shave_image(imread(GD_image_path),shave_width);

    % Calculating scores
    scores(ii).name = file_list(ii).name;
    scores(ii).MSE = immse(input_image,GD_image);
    scores(ii).Ma = quality_predict(input_image);
    scores(ii).NIQE = computequality(input_image,blocksizerow,blocksizecol,...
        blockrowoverlap,blockcoloverlap,mu_prisparam,cov_prisparam);
end

end

author: roimehrez


text:
%  Figure 4.9      Feedback Control of Dynamic Systems, 6e
%                   Franklin, Powell, Emami
% Figure 4.9   PID control of motor speed
% function [np, dp]= pid(b,J,K,L,R,kp,ki,kd)
% function to compute the equations of a d.c. motor with inductance.
% and compute the response under control
clf;
K=.0670; L1=0.1; J1=0.0113; R=0.45; b=0.0280;
kp=3; ki= 15; kd=0.3;
np=K;
dp=[L1*J1 R*J1+b*L1 R*b+K*K];
dclp=[L1*J1 R*J1+b*L1 R*b+K*K+K*kp];
nclp=K*kp;
nclpw=[L1 R];
dclpi=[L1*J1 R*J1+b*L1 R*b+K*K+K*kp K*ki];
nclpi=[K*kp K*ki];
nclpiw=[L1 R 0];
dclpid=[  L1*J1 R*J1+b*L1+K*kd R*b+K*K+K*kp K*ki];
nclpid=[K*kd K*kp K*ki];
nclpidw=[L1 R 0];
sysp=tf(nclp,dclp);
syspw=tf(nclpw,dclp);
syspi=tf(nclpi,dclpi);
syspiw=tf(nclpiw,dclpi);
syspid=tf(nclpid,dclpid);
syspidw=tf(nclpidw,dclpid);
figure(1)
t=0:.01:6;
[y1w,t]=step(syspw,t);
[y2w,t]=step(syspiw,t);
[y3w,t]=step(syspidw,t);
plot(t,y1w,t,y2w,t,y3w);
xlabel('Time (msec)');
ylabel('Amplitude');
title('Fig. 4.9(a) Response of P,PI, and PID control to a disturbance step')
gtext('P')
gtext('PI')
gtext('PID')
nicegrid;

figure(2)
[y1,t]=step(sysp,t);
[y2,t]=step(syspi,t);
[y3,t]=step(syspid,t);
plot(t,y1,t,y2,t,y3);
xlabel('Time (msec)');
ylabel('Amplitude');
title('Fig. 4.9 (b) Response of P,PI, and PID control to a reference step')
gtext('P')
gtext('PI')
gtext('PID')
nicegrid;

author: Sable

repo: mcbench-benchmarks


text:
function [W,CovW,X,CovX] = vbgppcamv_full(Y,D,inW,inX,covfuncW,logthetaW, ...
                                          covfuncX, logthetaX,varargin)
% [W,CovW,X,CovX] = vbgppcamv(Y,D,inW,inX,covfuncW,logthetaW,covfuncX, ...
% logthetaX,varargin)
%
% Variational Bayesian (VB) Gaussian process (GP) principal component
% analysis (PCA) with handling of missing values (MV).
%
% This version does NOT factorize with respect to the components.
%
% At the moment, pseudo inputs are not yet supported.
%

[M,N] = size(Y);

opts = struct( ...
    'init',          [],...
    'rotate', true, ...
    'autosavetime', 0,...
    'autosavefile', 'vbgppcamv_autosave',...
    'testset', [], ...
    'pseudodensityx', 1, ...
    'pseudodensityw', 1, ...
    'updatepseudox', true, ...
    'updatepseudow', true, ...
    'initpseudow', [], ...
    'initpseudox', [], ...
    'reconstruct', true, ...
    'loglikelihood', true, ...
    'updatehyper', 1, ...
    'maxsearchx', 10, ...
    'maxsearchw', 10, ...
    'checkgradx', false, ...
    'checkgradw', false, ...
    'maxiter', 100);

[ opts, errmsg, wrnmsg ] = argschk( opts, varargin{:} );
if ~isempty(errmsg), error( errmsg ), end
if ~isempty(wrnmsg), warning( wrnmsg ), end

Obs = ~isnan(Y);
NMobs = sum(Obs(:));

inWorig = inW;
inXorig = inX;

% Initialize posterior values
[W,varW,X,varX,tau] = initialize(opts.init, D, inW,inX,covfuncW, ...
                                 logthetaW,covfuncX, logthetaX);

% Remove empty rows/columns
rowrm = (colsum(Obs) == 0);
colrm = (rowsum(Obs) == 0);
if any(rowrm)
  disp('Removing empty rows');
  inW(:,rowrm) = [];
  Y(rowrm,:) = [];
  Obs(rowrm,:) = [];
  Morig = M;
  M = rows(Y);
  W(rowrm,:) = [];
  varW(rowrm,:) = [];
end
if any(colrm)
  disp('Removing empty columns');
  inX(:,colrm) = [];
  Y(:,colrm) = [];
  Obs(:,colrm) = [];
  Norig = N;
  N = cols(Y);
  X(:,colrm) = [];
  varX(:,colrm) = [];
end

log2pi = log(2*pi);

%%
%% Initialize posterior parameters

% $$$ % Number of pseudo inputs for each component
% $$$ opts.pseudodensityx = opts.pseudodensityx(:) .* ones(D,1);
% $$$ opts.pseudodensityw = opts.pseudodensityw(:) .* ones(D,1);
% $$$ Np = ceil(opts.pseudodensityx.*cols(inX)); % no. of pseudos for X components
% $$$ Mp = ceil(opts.pseudodensityw.*cols(inW)); % no. of pseudos for W components
% $$$
% $$$ dimX = rows(inX);   % dimensionality of the input space of X
% $$$ dimW = rows(inW);   % dimensionality of the input space of W
% $$$
% $$$ % Initialize pseudo inputs
% $$$ pseudoX = cell(D,1);
% $$$ pseudoW = cell(D,1);
% $$$ usepseudoX = 1==zeros(D,1);
% $$$ usepseudoW = 1==zeros(D,1);
% $$$ for d=1:D
% $$$   if Np(d)==N
% $$$     pseudoX{d} = inX; % full, don't use pseudo inputs
% $$$     usepseudoX(d) = false;
% $$$   else
% $$$     permN = randperm(N);
% $$$     pseudoX{d} = inX(:,permN(1:Np(d)));
% $$$     usepseudoX(d) = true;
% $$$   end
% $$$   if Mp(d)==M
% $$$     pseudoW{d} = inW; % full, don't use pseudo inputs
% $$$     usepseudoW(d) = false;
% $$$   else
% $$$     permM = randperm(M);
% $$$     pseudoW{d} = inW(:,permM(1:Mp(d)));
% $$$     usepseudoW(d) = true;
% $$$   end
% $$$ end
% $$$
% $$$ if ~isempty(opts.initpseudow)
% $$$   pseudoW = opts.initpseudow;
% $$$   usepseudoW(:) = true;
% $$$ end
% $$$ if ~isempty(opts.initpseudox)
% $$$   pseudoX = opts.initpseudox;
% $$$   usepseudoX(:) = true;
% $$$ end
% $$$
% $$$ cholKWpp = cell(D,1);
% $$$ cholKXpp = cell(D,1);

mu = zeros(M,1);
v_mu = zeros(M,1);

logP = -inf;
costcpu = 0;

% $$$ % Wp
% $$$ Wp = cell(D,1);
% $$$ CovWp = cell(D,1);
% $$$
% $$$ % Xp
% $$$ Xp = cell(D,1);
% $$$ CovXp = cell(D,1);

%% tau precision (inverse noise)
atau = 1e-4;
btau = 1e-4;
a_tau = nan; % posterior pdf parameter
b_tau = nan; % posterior pdf parameter
logtau = nan;

lastsave = cputime;

%%
%% VB learning
% $$$ K = myfeval(covfuncW{1}{:}, logthetaW{1}, inW, inW);
% $$$ figure
% $$$ imagesc(K)
% $$$ return

% Index order (component-wise, i.e., block-diagonal K)
indsW = reshape(1:(D*M), [M D]);
indsX = reshape(1:(D*N), [N D])';

CovW = eye(D*M);
CovX = eye(D*N);

% $$$ %% TEMP
% $$$ W(:) = 1;
% $$$ CovW(:) = 0;

% Initialize help variables
WW = zeros(D,D,M);
for m=1:M
  WW(:,:,m) = W(m,:)'*W(m,:) + CovW(indsW(m,:),indsW(m,:));
end
XX = zeros(D,D,N);
for n=1:N
  XX(:,:,n) = X(:,n)*X(:,n)' + CovX(indsX(:,n),indsX(:,n));
end

% Evaluate prior covariance matrices
KW = zeros(D*M);
KX = zeros(D*N);
for d=1:D
  first = (d-1)*M + 1;
  last = first + M - 1;
  inds = first:last;
  covfunc = covfuncW{d};
  if ~iscell(covfunc)
    covfunc = {covfunc};
  end
  KW(inds,inds) = feval(covfunc{:}, logthetaW{d}, inW, inW);

  first = (d-1)*N + 1;
  last = first + N - 1;
  inds = first:last;
  covfunc = covfuncX{d};
  if ~iscell(covfunc)
    covfunc = {covfunc};
  end
  KX(inds,inds) = feval(covfunc{:}, logthetaX{d}, inX, inX);
end
% $$$ KW = regularize(KW, 1e1);
% $$$ LKW = chol(KW, 'lower');
% $$$ KX = regularize(KX, 1e2);
% $$$ LKX = chol(KX, 'lower');

% $$$ figure
% $$$ imagesc(KW)
% $$$ figure
% $$$ imagesc(KX)

reg = 1e-6;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VB EM

for ind=1:opts.maxiter

  oldW = W;

  startitercpu = cputime;

  %% UPDATE VARIABLES

  update = false;
  if length(opts.updatehyper) > 1
    if any(ind==opts.updatehyper)
      update = true;
    end
  else
    if mod(ind, opts.updatehyper) == 0 && opts.updatehyper ~=0
      update = true;
    end
  end
  if update
    maxsearchx = opts.maxsearchx;
    maxsearchw = opts.maxsearchw;
  else
    maxsearchx = 0;
    maxsearchw = 0;
  end

  if ind < 5 && false
    % Scale X to unit variance. I think it should improve the lower bound
    % but it doesn't seem to do so in practice (at least always).. ??
    disp('Using adhoc whitening of second moment')
    %xx = diag(sum(XX,3))/N - mean(X,2).^2;   % variance
    xx = diag(sum(XX,3))/N;   % second moment
    sqrtxx = sqrt(xx);
    X = bsxfun(@rdivide, X, sqrtxx);
    XX = bsxfun(@rdivide, XX, (sqrtxx(:)*sqrtxx(:)'));
    W = bsxfun(@times, W, sqrtxx');
    WW = bsxfun(@times, WW, (sqrtxx(:)*sqrtxx(:)'));
  end
  %XX_norm = mean(XX,3)

  %% Update W
  U = spalloc(D*M,D*M, M*D*D);
  z = zeros(D*M,1);
  for m=1:M
    obs = Obs(m,:);
    U(indsW(m,:),indsW(m,:)) = tau * sum(XX(:,:,obs),3);
    z(indsW(m,:)) = tau * X(:,obs) * Y(m,obs)';
  end
  % Optimise hyperparameters
  if update
    logthetaW = optimise_hyperparameter(logthetaW, covfuncW, inW, U, z, D, M);
  end
  % Update prior covariance matrix
  KW = zeros(D*M);
  for d=1:D
    first = (d-1)*M + 1;
    last = first + M - 1;
    inds = first:last;
    covfunc = covfuncW{d};
    if ~iscell(covfunc)
      covfunc = {covfunc};
    end
    KW(inds,inds) = feval(covfunc{:}, logthetaW{d}, inW, inW);

    KW = KW + reg*eye(size(KW));

  end
  %KW = regularize(KW, 1e1);
  %LKW = chol(KW, 'lower');
  % Evaluate posterior
  if false % choose method
    L = chol(KW + KW*U*KW, 'lower');
    sqrtCovW = solve_tril(L, KW);
    CovW = sqrtCovW' * sqrtCovW;
    W(:) = CovW(indsW(:),indsW(:)) * z(indsW(:));
  else
    % Prefer this!
%    L = safechol(KW + inv(U), 1e-8, 'lower');
    L = chol(KW + inv(U), 'lower');
    tmp = solve_tril(L, KW);
    CovW = KW - tmp' * tmp;
    tmp = linsolve_chol(L, U \ z, 'lower');
    W(:) = KW(indsW(:),indsW(:)) * tmp(indsW(:));
  end
  WW = zeros(D,D,M);
  for m=1:M
    WW(:,:,m) = W(m,:)'*W(m,:) + CovW(indsW(m,:),indsW(m,:));
  end
% $$$   figure
% $$$   imagesc(WW(:,:,1))
% $$$   error('jou')

% $$$   figure
% $$$   imagesc(CovW)
% $$$   error('halo')

  %% Update X
  U = spalloc(D*N,D*N, N*D*D);
  z = zeros(D*N,1);
  for n=1:N
    obs = Obs(:,n);
    U(indsX(:,n),indsX(:,n)) = tau * sum(WW(:,:,obs),3);
    z(indsX(:,n)) = tau * W(obs,:)' * Y(obs,n);
  end
  % Optimise hyperparameters
  if update
    logthetaX = optimise_hyperparameter(logthetaX, covfuncX, inX, U, z, D, N);
  end
  % Update prior covariance matrix
  KX = zeros(D*N);
  for d=1:D
    first = (d-1)*N + 1;
    last = first + N - 1;
    inds = first:last;
    covfunc = covfuncX{d};
    if ~iscell(covfunc)
      covfunc = {covfunc};
    end
    KX(inds,inds) = feval(covfunc{:}, logthetaX{d}, inX, inX);

    KX = KX + reg*eye(size(KX));

  end
  %KX = regularize(KX, 1e2);
  %LKX = chol(KX, 'lower');
  % Evaluate posterior
% $$$   figure
% $$$   imagesc(inv(U))
% $$$   figure
% $$$   imagesc(KX)
% $$$   tmp = chol(U, 'lower');
% $$$   figure
% $$$   imagesc(linsolve_chol(tmp, eye(size(tmp))) == 0);
  if false
    L = chol(KX + KX*U*KX, 'lower');
    sqrtCovX = solve_tril(L, KX);
    CovX = sqrtCovX' * sqrtCovX;
    X(:) = CovX(indsX(:),indsX(:)) * z(indsX(:));
  else
    % Prefer this!
    L = chol(KX + inv(U), 'lower');
%    L = safechol(KX + inv(U), 1e-8, 'lower');
    tmp = solve_tril(L, KX);
    CovX = KX - tmp' * tmp;
    tmp = linsolve_chol(L, U \ z, 'lower');
    X(:) = KX(indsX(:),indsX(:)) * tmp(indsX(:));
  end
  XX = zeros(D,D,N);
  for n=1:N
    XX(:,:,n) = X(:,n)*X(:,n)' + CovX(indsX(:,n),indsX(:,n));
  end

  %XX_norm = mean(XX,3)

%  CovX(1:10,1:10)
%  tmp = chol(KX, 'lower');

% $$$   figure
% $$$   imagesc(CovX)
% $$$   error('halo')

% $$$   rcondX = rcond(CovX)
% $$$   rcondestX = rcond(CovX + 1e-4*normest(CovX)*eye(size(CovX)))
% $$$   tmp = chol(CovX,'lower');

%  [X,CovX,W,CovW] = rotate(X,CovX,indsX,LKX,W,CovW);

% $$$   for d=1:D
% $$$     others = [1:(d-1), (d+1):D];
% $$$     V = spalloc(M,M,M);
% $$$     c = zeros(M,1);
% $$$     for m=1:M
% $$$       obs = Obs(m,:);
% $$$       V(m,m) = tau * (X(d,obs)*X(d,obs)' + sum(varX(d,obs)));
% $$$       c(m) = tau * X(d,obs) * (Y(m,obs) - W(m,others)*X(others,obs))';
% $$$     end
% $$$
% $$$     if usepseudoW(d)
% $$$       pseudos = pseudoW{d};
% $$$     else
% $$$       pseudos = [];
% $$$     end
% $$$     [Wp{d}, CovWp{d}, pseudos, logthetaW{d}, tmp, Kwp, cholKWpp{d}] = ...
% $$$         gplearn(logthetaW{d}, covfuncW{d}, inW, [], V, pseudos, 'vy', c, ...
% $$$                 'maxsearch', maxsearchw, 'checkgrad', opts.checkgradw, ...
% $$$                 'updatepseudo', opts.updatepseudow);
% $$$
% $$$     if ~usepseudoW(d)
% $$$       % full
% $$$       W(:,d) = Wp{d};
% $$$       varW(:,d) = diag(CovWp{d});
% $$$     else
% $$$       % using pseudo inputs
% $$$       pseudoW{d} = pseudos;
% $$$       [W(:,d), varW(:,d)] = gppred(pseudoW{d}, Wp{d}, CovWp{d}, inW, ...
% $$$                                    logthetaW{d}, covfuncW{d}, 'cholkx', ...
% $$$                                    cholKWpp{d}, 'khx', Kwp);
% $$$     end
% $$$     ltW_in_vbgppca = logthetaW{d};
% $$$   end

% $$$   %% Update X
% $$$   for d=1:D
% $$$     others = [1:(d-1), (d+1):D];
% $$$     V = spalloc(N,N,N);
% $$$     c = zeros(N,1);
% $$$     for n=1:N
% $$$       obs = Obs(:,n);
% $$$       V(n,n) = tau * (W(obs,d)'*W(obs,d) + sum(varW(obs,d)));
% $$$       c(n) = tau * W(obs,d)' * (Y(obs,n) - W(obs,others)*X(others,n));
% $$$     end
% $$$     if usepseudoX(d)
% $$$       pseudos = pseudoX{d};
% $$$     else
% $$$       pseudos = [];
% $$$     end
% $$$     [Xp{d}, CovXp{d}, pseudos, logthetaX{d}, tmp, Kxp, cholKXpp{d}] = ...
% $$$         gplearn(logthetaX{d}, covfuncX{d}, inX, [], V, pseudos, 'vy', c, ...
% $$$                 'maxsearch', maxsearchx, 'checkgrad', opts.checkgradx, ...
% $$$                 'updatepseudo', opts.updatepseudox);
% $$$
% $$$     if ~usepseudoX(d)
% $$$       % full
% $$$       X(d,:) = Xp{d};
% $$$       varX(d,:) = diag(CovXp{d});
% $$$     else
% $$$       % pseudo inputs
% $$$       pseudoX{d} = pseudos;
% $$$       [X(d,:), varX(d,:)] = gppred(pseudoX{d}, Xp{d}, CovXp{d}, inX, ...
% $$$                                    logthetaX{d}, covfuncX{d}, 'cholkx', ...
% $$$                                    cholKXpp{d}, 'khx', Kxp);
% $$$     end
% $$$
% $$$     ltX_in_vbgppca = logthetaX{d};
% $$$   end

  disp('Updating tau')

  %% Update tau
  % TODO: You could loop over the longer dimension N or M? Then
  % vectorization would help the most. :)
  err2 = 0;
  for n=1:N
    obs = Obs(:,n);
    ymu = Y(obs,n) - mu(obs);
    err2 = err2 ...
           + ymu'*ymu ...
           + traceprod(sum(WW(:,:,obs),3), XX(:,:,n)) ...
           + sum(v_mu) ...
           - 2*(W(obs,:)*X(:,n))' * ymu;
  end
  a_tau = atau + 0.5*NMobs;
  b_tau = btau + 0.5*err2;
  tau = a_tau / b_tau;
  logtau = psi(a_tau) - log(b_tau);

  % CPU time used for calculations
  itercpu = cputime - startitercpu;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %% CALCULATE LOWER BOUND OF THE LOG-LIKELIHOOD

  disp('Evaluating lower bound.');

  if opts.loglikelihood

    startcostcpu = cputime;
    oldlogP = logP;

    % Cost from Y
    cost_Y = -0.5*tau*err2 - 0.5*NMobs*(log2pi - logtau);

    % WTF!?!?!?!!! If using safechol, the coefficient has a HUUUGE effect
    % on the loglikelihood lower bound!!!

    % Cost from W
    %rcond_KW = rcond(KW)
%    LKW = chol(KW(indsW(:),indsW(:)), 'lower');
    LKW = chol(KW(indsW(:),indsW(:))+reg*eye(size(KW)), 'lower');
%    LKW = safechol(KW(indsW(:),indsW(:)), 1e-15, 'lower');
    cost_W = mvnvbcost(W(:),CovW(indsW(:),indsW(:)), 0,0, LKW, 2*logdettri(LKW));

% $$$     cost_W = 0;
% $$$     for d=1:D
% $$$       cost_W = cost_W + mvnvbcost(Wp{d},CovWp{d}, 0,0, cholKWpp{d}, ...
% $$$                                   2*logdettri(cholKWpp{d}));
% $$$     end

    % Cost from X
%    rcond_KX = rcond(KX)
%    LKX = chol(KX(indsX(:),indsX(:)), 'lower');
    LKX = chol(KX(indsX(:),indsX(:))+reg*eye(size(KX)), 'lower');
%    LKX = safechol(KX(indsX(:),indsX(:)), 1e-15, 'lower');
    cost_X = mvnvbcost(X(:),CovX(indsX(:),indsX(:)), 0,0, LKX, 2*logdettri(LKX));
% $$$     cost_X = 0;
% $$$     for d=1:D
% $$$       cost_X = cost_X + mvnvbcost(Xp{d},CovXp{d}, 0,0, cholKXpp{d}, ...
% $$$                                   2*logdettri(cholKXpp{d}));
% $$$     end

    % Cost from tau
    cost_tau = -gamkl(tau,logtau,a_tau,b_tau,atau,btau);
% $$$
% $$$     oldlogP = logP;
% $$$
% $$$     % Lower bound of loglikelihood
    logP = cost_Y + cost_W + cost_X + cost_tau;% + cost_mu + cost_tau + cost_w;
    costcpu = cputime - startcostcpu;

    if logP < oldlogP
      warning('Cost increased!!');
    end

  end

  disp('Monitoring stuff.')

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %% MONITORING STUFF

  % Change in subspace
  angle = subspace(oldW, W);

  %% Print progress
  fprintf('%d step: loglike=%e, angle=%.2e (itercpu=%.2fs, nllcpu=%.2fs)\n', ...
          ind, logP, angle, itercpu, costcpu);

  %% Save
  loglikelihood = logP;
  if opts.autosavetime > 0 && cputime - lastsave > opts.autosavetime
    fprintf('Saving to a file %s..', opts.autosavefile);
    iteration = ind;
    save(opts.autosavefile, 'Wp', 'CovWp', 'W', 'varW', 'pseudoW', ...
         'logthetaW', 'covfuncW', 'inW', 'Xp', 'CovXp', 'X', 'varX', ...
         'pseudoX', 'logthetaX', 'covfuncX', 'inX', 'tau', 'a_tau', 'b_tau', ...
         'loglikelihood', 'ind');
    fprintf(' done\n');
    lastsave = cputime;
  end

end

% $$$ figure
% $$$ imagesc(CovX == 0)

% Predict whole dataset, also rows and columns that were removed at the
% beginning because of they were empty
if opts.reconstruct
  if any(rowrm)
    % Change variables
    W_old = W;
    KW_old = KW;
    CovW_old = CovW;
    Morig = size(inWorig,2);
    W = zeros(Morig,D);
    KW = zeros(D*Morig);
    KWWold = zeros(D*Morig,D*M);
    CovW = zeros(D*Morig);
    indsWorig = reshape(1:(D*Morig), [Morig D]);
    % New prior covariance matrix
    for d=1:D
      first = (d-1)*Morig + 1;
      last = first + Morig - 1;
      inds = first:last;
      first = (d-1)*Morig + 1;
      last = first + Morig - 1;
      inds_old = first:last;
      covfunc = covfuncW{d};
      if ~iscell(covfunc)
        covfunc = {covfunc};
      end
      KWWold(inds,inds_old) = feval(covfunc{:}, logthetaW{d}, inWorig, inW);
      KW(inds,inds) = feval(covfunc{:}, logthetaW{d}, inWorig, inWorig);
      KW = KW + reg*eye(size(KW));
    end
    % New posterior quantities
    LKW_old = chol(KW_old(indsW(:),indsW(:)), 'lower');
    W(:) = KWWold(indsWorig(:),indsW(:)) * linsolve_chol(LKW_old, W_old(:), ...
                                                      'lower');
    A = linsolve_chol(LKW_old, KWWold(indsWorig(:),indsW(:))', 'lower');
    CovW(indsWorig(:),indsWorig(:)) = KW(indsWorig(:),indsWorig(:)) - ...
        A'*(KW_old(indsW(:),indsW(:))-CovW(indsW(:),indsW(:)))*A;

    M = Morig;
    indsW = indsWorig;
    inW = inWorig;
  end
  if any(colrm)
% $$$     figure
% $$$     imagesc(CovX)

    % Change variables
    X_old = X;
    KX_old = KX;
    CovX_old = CovX;
    Norig = size(inXorig,2);
    X = zeros(D,Norig);
    KX = zeros(D*Norig);
    KXXold = zeros(D*Norig,D*N);
    CovN = zeros(D*Norig);
    indsXorig = reshape(1:(D*Norig), [Norig D])';
    % New prior covariance matrix
    for d=1:D
% $$$       first = (d-1)*Morig + 1;
% $$$       last = first + Morig - 1;
% $$$       inds = first:last;
      first = (d-1)*Norig + 1;
      last = first + Norig - 1;
      inds = first:last;
      first = (d-1)*N + 1;
      last = first + N - 1;
      inds_old = first:last;
      covfunc = covfuncX{d};
      if ~iscell(covfunc)
        covfunc = {covfunc};
      end
      KXXold(inds,inds_old) = feval(covfunc{:}, logthetaX{d}, inXorig, inX);
      KX(inds,inds) = feval(covfunc{:}, logthetaX{d}, inXorig, inXorig);
      KX = KX + reg*eye(size(KX));
    end
    % New posterior quantities
    LKX_old = chol(KX_old(indsX(:),indsX(:)), 'lower');
    X(:) = KXXold(indsXorig(:),indsX(:)) * linsolve_chol(LKX_old, X_old(:), ...
                                                      'lower');
    A = linsolve_chol(LKX_old, KXXold(indsXorig(:),indsX(:))', 'lower');
    CovX(indsXorig(:),indsXorig(:)) = KX(indsXorig(:),indsXorig(:)) - ...
        A'*(KX_old(indsX(:),indsX(:))-CovX(indsX(:),indsX(:)))*A;

% $$$     figure
% $$$     imagesc(CovX)

    N = Norig;
    indsX = indsXorig;
    inX = inXorig;
  end
end
% $$$ if opts.reconstruct
% $$$   if any(rowrm)
% $$$     oldW = W;
% $$$     KW_old = KW;
% $$$     KW = feval(covfunc{:}, logthetaW{d}, inWorig, inWorig);
% $$$     W = zeros(Morig,D);
% $$$     varW = zeros(Morig,D);
% $$$     for d=1:D
% $$$       [W(:,d), varW(:,d)] = gppred(inW, oldW(:,d), CovW(indsW(:,d),indsW(:,d)), ...
% $$$                                    inWorig, logthetaW{d}, covfuncW{d});
% $$$     end
% $$$   end
% $$$   if any(colrm)
% $$$     oldX = X;
% $$$     X = zeros(D,Norig);
% $$$     varX = zeros(D,Norig);
% $$$     for d=1:D
% $$$       [X(d,:), varX(d,:)] = gppred(inX, oldX(d,:)', CovX(indsX(d,:),indsX(d,:)), ...
% $$$                                    inXorig, logthetaX{d}, covfuncX{d});
% $$$     end
% $$$   end
% $$$ end

varX = reshape(diag(CovX), [N D])';
varW = reshape(diag(CovW), [M D]);

if nargout == 1
  Q.W = W;
  Q.CovW = CovW;
% $$$   Q.W = W;
  Q.varW = varW;
% $$$   Q.pseudoW = pseudoW;
  Q.logthetaW = logthetaW;
  Q.covfuncW = covfuncW;
  Q.inW = inWorig;

  Q.X = X;
  Q.CovX = CovX;
% $$$   Q.X = X;
  Q.varX = varX;
% $$$   Q.pseudoX = pseudoX;
  Q.logthetaX = logthetaX;
  Q.covfuncX = covfuncX;
  Q.inX = inXorig;

  Q.tau = tau;
  Q.a_tau = a_tau;
  Q.b_tau = b_tau;

  Q.indsW = indsW;
  Q.indsX = indsX;

  Q.loglikelihood = loglikelihood;

  W = Q;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function logtheta = optimise_hyperparameter(logtheta, covfuncs, input, U, ...
                                            z, D, N)

% Because the cell logtheta must be presented as one long vector, mark
% the indeces for each component. That is, the hyperparameters of d-th
% component are vector( (indeces(d-1)+1):indeces(d) ) = logtheta{d}
vector = [];
indeces = zeros(D,1);
for d=1:D
  vector = [vector; logtheta{d}(:)];
  indeces(d) = length(logtheta{d});
end
indeces = cumsum(indeces);

invU = inv(U);

%hyperparameters_before = exp(vector(:))

% $$$ mycheckgrad(@hyperparameter_upperbound, vector, 1e-9, covfuncs, input, ...
% $$$             invU, z, D, N, indeces);

vector = minimize(vector, @hyperparameter_upperbound, 5, covfuncs, input, ...
                  invU, z, D, N, indeces);

%hyperparameters_after = exp(vector(:))

ind = 1;
for d=1:D
  logtheta{d} = vector(ind:indeces(d));
  ind = indeces(d) + 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [f, df] = hyperparameter_upperbound(logtheta_vector, covfuncs, ...
                                             input, invU, z, D, N, ...
                                             logtheta_indeces)

% Projected pseudo-observations
y = invU * z;

% Transform the hyperparameter vector to cell
ind = 1;
for d=1:D
  logtheta{d} = logtheta_vector(ind:logtheta_indeces(d));
  ind = logtheta_indeces(d) + 1;
end

dlogtheta = cell(D,1);

K = zeros(D*N);
for d=1:D
  % Indeces of a component block
  first = (d-1)*N + 1;
  last = first + N - 1;
  inds = first:last;
  % Covariance matrix for the block
  covfunc = covfuncs{d};
  if ~iscell(covfunc)
    covfunc = {covfunc};
  end
  [K(inds,inds), dlogtheta{d}] = feval(covfunc{:}, logtheta{d}, input, input);
end
%K = regularize(K, 1e2);
%LK = chol(K, 'lower');
%K =

%if rcond(K) < 1e-12
%  f = -

L = chol(K + invU, 'lower');

% Loglikelihood lower bound (mnorm_lpdf)
%f = mnorm_lpdf((invU*z)', 0, K + invU);
v = solve_tril(L, (y-0));
f = -0.5*length(L)*log(2*pi) - logdettri(L) - 0.5 * v'*v;

% Gradients
df = [];
% TODO: USE LINSOLVE TO OPTIMISE BACKSUBSTITUTIONS!!!!
%invLL = L' \ (L \ eye(size(L)));
%c = invLL * y;
c = linsolve_chol(L, y, 'lower');
for d=1:D
  % Indeces of a component block
  first = (d-1)*N + 1;
  last = first + N - 1;
  inds = first:last;
  Ld = L(inds,inds);
  invL = linsolve_chol(Ld, eye(size(Ld)), 'lower');
  % Evaluate gradients
  for i=1:length(logtheta{d})
    dK = dlogtheta{d}(:,:,i);
    %sz_c = size(c)
    %sz_dK = size(dK)
    df = [df; 0.5*c(inds)'*dK*c(inds) - 0.5*traceprod(invL,dK)];
  end
end

% Lower bound to upper bound (because using minimize)
f = -f;
df = -df;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X,CovX,W,CovW] = rotate(X,CovX,indsX,LKX,W,CovW)

error('Ei taa oikein toimi, ei kannata kayttaa..')

disp('Finding rotation..')
[D,N] = size(X);
R = orth(randn(D));
R = R*R';

%mycheckgrad(@rotation_cost, R(:), 1e-6, X,CovX,indsX,LKX);

R = minimize(R(:), @rotation_cost, 3, X,CovX,indsX,LKX);

R = reshape(R,[D,D]);

% $$$ % Force positive definiteness
% $$$ [VR,DR] = eig(R);
% $$$ R = VR * exp(DR) * inv(VR);

X = R * X;
kronR = kron(R,eye(N));
CovX = kronR * CovX * kronR';

M = rows(W);
invR = inv(R);
W = W * invR;
kronR = kron(invR,eye(M));
CovW = kronR' * CovW * kronR;
disp('Rotation done.');

R
%error('done')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [f,df] = rotation_cost(R,X,CovX,indsX,LKX)
[D,N] = size(X);
R = reshape(R,[D D]);

% Force positive definiteness
% $$$ [VR,DR] = eig(R);
% $$$ R = VR * exp(DR) * inv(VR);

log_qX = -N * logdet(R);
if ~isreal(log_qX)
  f = inf;
  df = nan*ones(D*D,1);
  return
end

LinvKX = solve_tril(LKX, eye(D*N));
invKX = LinvKX' * LinvKX;
%invKX = solve_triu(LKX', solve_tril(LKX, eye(D*N)));
H = 0;
for n=1:N
  for m=1:N
    in = indsX(:,n);
    im = indsX(:,m);
    H = H + invKX(in,im) * R * (X(:,m)*X(:,n)' + CovX(im,in));
  end
end
error('Check that traceprod')
log_pX = -0.5 * traceprod(H,R);

f = -(log_pX - log_qX)
df = N*inv(R') - H;
df = -df(:);

% Test:
X = R * X;
kronR = kron(R,eye(N));
CovX = kronR * CovX * kronR';
cost_X = mvnvbcost(X(:),CovX(indsX(:),indsX(:)), 0,0, LKX, 2*logdettri(LKX))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ind = row2ind(M,N,m)
ind = (0:(N-1)) * M + m;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ind = col2ind(M,N,n)
ind = (1:M) + (n-1)*M;

%%%%%%%%%%%%%%%%%%%%%
function Y = gpinv(X)
[Q,R] = qr(X);
Y = inv(R) * Q';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [W,varW,X,varX,tau] = initialize(init,D,inW,inX,covfuncW, ...
                                          logthetaW,covfuncX, logthetaX)

M = cols(inW);
N = cols(inX);

% W
if isstruct(init) && isfield(init, 'W') && ~isempty(init.W)
  disp('Initialize W to given value');
  W = init.W;
else
%  W = randn(M,D);
  disp('Initialize W to zero');%by sampling');
  W = zeros(M,D);
% $$$   for d=1:D
% $$$     W(:,d) = gprnd(inW, logthetaW{d}, covfuncW{d});
% $$$   end
end

% varW
if isstruct(init) && isfield(init, 'varW') && ~isempty(init.varW)
  varW = init.varW;
else
  varW = 1*ones(M,D);
end

% X
if isstruct(init) && isfield(init, 'X') && ~isempty(init.X)
  disp('Initialize X to given value');
  X = init.X;
else
  %X = randn(D,N);
  disp('Initialize X by sampling');
  X = zeros(D,N);
  for d=1:D
    X(d,:) = gprnd(inX, logthetaX{d}, covfuncX{d});
  end
end

% varX
if isstruct(init) && isfield(init, 'varX') && ~isempty(init.varX)
  varX = init.varX;
else
  varX = 1 * ones(D,N);
end

%% tau precision (inverse noise)
if isstruct(init) && isfield(init, 'tau') && ~isempty(init.tau)
  tau = init.tau;
else
  tau = 1e3;
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function cdf = burr_cdf ( x, a, b, c, d )

%*****************************************************************************80
%
%% BURR_CDF evaluates the Burr CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Input, real A, B, C, D, the parameters of the PDF.
%    0 < B,
%    0 < C.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x <= a )

    cdf = 0.0;

  else

    cdf = 1.0 / ( 1.0 + ( b / ( x - a ) )^c  )^d;

  end

  return
end


text:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.com
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.

function CstDefaultUnits(mws)

%define units (Default selection)
Geometry = 'mm';
Frequency = 'GHz';
Time = 'ns';
TemperatureUnit = 'Kelvin';
Voltage = 'V';
Current = 'A';
Resistance = 'Ohm';
Conductance  = 'S';
Capacitance = 'PikoF';
Inductance = 'NanoH';

Units = invoke(mws,'Units');
invoke(Units,'Geometry',Geometry);
invoke(Units,'Frequency',Frequency);
invoke(Units,'Time',Time);
invoke(Units,'TemperatureUnit',TemperatureUnit);
invoke(Units,'Voltage',Voltage);
invoke(Units,'Current',Current);
invoke(Units,'Resistance',Resistance);
invoke(Units,'Conductance',Conductance);
invoke(Units,'Capacitance',Capacitance);
invoke(Units,'Inductance',Inductance);

end

repo: CST-MATLAB-API


text:
function [ x, seed ] = latin_center ( dim_num, point_num, seed )

%*****************************************************************************80
%
%% LATIN_CENTER returns center points in a Latin square.
%
%  Discussion:
%
%    In each spatial dimension, there will be exactly one
%    point with the coordinate value
%
%      ( 1, 3, 5, ..., 2*point_num-1 ) / ( 2 * point_num )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer dim_num, the spatial dimension.
%
%    Input, integer point_num, the number of points.
%
%    Input, integer SEED, a seed for UNIFORM, the random number generator.
%
%    Output, real x(dim_num,point_num), the points.
%
%    Output, integer SEED, the updated random number seed.
%
  base = 1;

  for i = 1: dim_num

    [ perm, seed ] = perm_uniform ( point_num, base, seed );

    for j = 1: point_num
      x(i,j) = ( 2 * perm(j) - 1 ) / ( 2 * point_num );
    end

  end

  return
end


text:
function [U_exact] = exact_translating_sphere(sphere_translational_velocity,sphere_radius,mesh_eval)
%+========================================================================+
%|                                                                        |
%|           This function uses the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT :                                                            |
%| PROPERTY  :                                                            |
%| LICENCE   :                                                            |
%| CONTACT   :                                                            |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : exact_translating_sphere.m                    |
%|    #    |   VERSION    : 0.10                                          |
%|   _#_   |   AUTHOR(S)  : Luca Berti                                    |
%|  ( # )  |   CREATION   : 25.12.2018                                    |
%|  / 0 \  |   LAST MODIF :                                               |
%| ( === ) |   SYNOPSIS   :                                               |
%|  `---'  |                                                              |
%+========================================================================+
% COMPUTE THE EXACT VELOCITY FIELD DETERMINED BY A TRANSLATING SPHERE IN
% AN INFINITE FLUID ALONG THE Z AXIS
% Reference: Happel-Brenner, "Low Reynolds number hydrodynamics"
% (pag. 119 and formulas 4-17.17 and 4-17.18) transformed into Cartesian
% coordinates (p.507 A-15.29 for the transformation)

% Input -> mesh_eval: mesh containing the evaluation points for the exact
%                     solution
%          sphere_translational_velocity -> translational velocity of the
%                                           vertically moving sphere
%                                           [0,0,W]
% Output -> U_exact:  3xN_eval vector containing the values of the exact
%                     solution at the evaluation point, expressed IN CARTESIAN
%                     COORDINATES

W = sphere_translational_velocity(3);
% Exact solutions taken from the reference - SPHERICAL COORDINATES
u_r = @(x,y,z) -W.*cos(acos(z./sqrt(x.^2+y.^2+z.^2))).*(sphere_radius^3./(2*sqrt(x.^2+y.^2+z.^2).^3) -3*sphere_radius./(2*sqrt(x.^2+y.^2+z.^2)));
u_theta = @(x,y,z) -W.*sin(acos(z./sqrt(x.^2+y.^2+z.^2))).*(sphere_radius^3./(4*sqrt(x.^2+y.^2+z.^2).^3)+3*sphere_radius./(4*sqrt(x.^2+y.^2+z.^2)));
u_phi =@(x,y,z) 0;

X = mesh_eval.vtx;
N_eval = size(X,1);

% Evaluation vector - SPHERICAL COORDINATES
U_exact = zeros(N_eval*3,1);
U_exact(1:N_eval,1) = feval(u_r,X(:,1),X(:,2),X(:,3));
U_exact(N_eval+1:2*N_eval,1) = feval(u_theta,X(:,1),X(:,2),X(:,3));
U_exact(2*N_eval+1:end,1) = feval(u_phi,X(:,1),X(:,2),X(:,3));

% Expression of the evaluation points in spherical coordinates, to extract
% (r,theta,phi) needed to transform the exact velocities from spherical to
% Cartesian coordinates
Y = X;
[X(:,3),X(:,2),X(:,1)] = cart2sph(Y(:,1),Y(:,2),Y(:,3));
% Latitude is needed instead of colatitude (as computed by Matlab)
X(:,2) = pi/2-X(:,2);

temporary = zeros(N_eval*3,1);
% Transformation of the exact velocity vector field from spherical to
% Cartesian coordinates - CARTESIAN COORDINATES
temporary(1:N_eval,1) = U_exact(1:N_eval,1).*sin(X(:,2)).*cos(X(:,3)) + ...
                          U_exact(N_eval+1:N_eval*2,1).*cos(X(:,2)).*cos(X(:,3))- ...
                          U_exact(1+2*N_eval:end,1).*sin(X(:,3));

temporary(1+N_eval:2*N_eval,1) = U_exact(1:N_eval,1).*sin(X(:,2)).*sin(X(:,3)) + ...
                          U_exact(N_eval+1:N_eval*2,1).*cos(X(:,2)).*sin(X(:,3)) + ...
                          U_exact(1+2*N_eval:end,1).*cos(X(:,3));

temporary(1+2*N_eval:end,1) = U_exact(1:N_eval,1).*cos(X(:,2)) - ...
                          U_exact(N_eval+1:N_eval*2,1).*sin(X(:,2));
U_exact = temporary;

end