text:
function par_filehost(varargin)
% This function starts a file host on the network.
%
% Notes:
%     The request format accepted by this host on a UDP broadcast port (default 34575) is (without <>'s)
%     <command>\n
%     <argument1>=<value1>
%     <argument2>=<value2>
%     <argument3>=<value3>
%     ...
%
%     Supported commands:
%     getfile -- return a file, either from memory cache or disk (if this host is responsible)
%     filname=<name of file to retrieve or expression that can be evaluated>
%     force=<whether to force-load the data even if it would exceed the allotted capacity>
%     return_address=<IP address or hostname to send the result to>
%     return_port=<port number of return endpoint>
%
%     evict -- evict file from memory cache if present
%     filename=<name of file to evict from memory cache>
%     anyhost=<evict even if not the responsible host>
%
%     terminate -- terminate this file host
%     return_address=<IP address or hostname to send the result to>
%     return_port=<port number of return endpoint>
%
%     getrank -- get the rank of this host among all known file server hosts
%     return_address=<IP address or hostname to send the result to>
%     return_port=<port number of return endpoint>
%
%     allhosts -- get the list of hosts as a cell array
%     return_address=<IP address or hostname to send the result to>
%     return_port=<port number of return endpoint>
%
%     hostname -- get the hostname of this machine
%     return_address=<IP address or hostname to send the result to>
%     return_port=<port number of return endpoint>
%
%     max_capacity -- get the allowed cache capacity of this machine in GB
%     return_address=<IP address or hostname to send the result to>
%     return_port=<port number of return endpoint>
%
%     clear -- clear cache
%     return_address=<IP address or hostname to send the result to>
%     return_port=<port number of return endpoint>
%
% See also:
%   par_fileclient, io_load
%
%                               Christian Kothe, Syntrogi
%                               2015-06-27

% the database (not persistent so we can change the code without erasing it)
global db;      % struct mapping of cache-tag --> data-blob
global lru;     % cell array of recently-used cache tags (sorted by increasing age)
if ~iscell(lru)
    lru = {}; end

fprintf('parsing arguments...');
settings = arg_define(varargin, ...
    arg({'all_hosts','AllHosts'}, {}, [], 'Hostnames where filehosts are running. Used for load-balancing. If this is set to empty, then we are only accepting connections from localhost.'), ...
    arg({'rank_override','RankOverride'}, 0, int32([0 1000]), 'Rank of this instance (1-based). Only used if the current hostname is not in the set AllHosts.'), ...
    arg({'max_capacity','MaxCapacity'}, 90, [1 10 256 100000], 'Maximum memory capacity. In Gigabytes.'), ...
    arg({'accept_pattern','AcceptPattern'}, '*', [], 'File name pattern that is accepted by this host.'), ...
    arg({'per_pattern_capacity','PerPatternCapacity'}, {}, [], 'Per-pattern capacity. This is a cell array of {pattern,capacity,pattern,capacity,...} where pattern are wildcard patterns.','type','expression'), ...
    arg({'evict_policy','EvictPolicy'}, 'lru', {'lru','none'}, 'Eviction policy. What entries will be evicted from memory if capacity is exceeded. If none, no more files will be cached after capacity has been exceeded.'), ...
    arg({'request_port','RequestPort'}, 34575, uint32([1 65535]), 'UDP port on which requests are being received.'), ...
    arg({'timeout_backconnect','TimeoutBackconnect',},1000,[],'Timeout for back-connection. In miliseconds. If it takes longer than this to connect back to the requester we give up.','guru',true), ...
    arg({'timeout_transmit','TimeoutTransmit',},60000,[],'Timeout for transmission. In miliseconds. If we get a hiccup during back-transmission that is longer than this we give up.','guru',true), ...
    arg({'accept_packet_size','AcceptPacketSize'},1400,uint32([1 508 1472 65507]),'Maximum packet size accepted by this host.','guru',true));
fprintf('done.\n');

% maximum amount of memory used inside the Java VM
max_java_memory = 2^26; % 64 MB

fprintf('Importing java packages...');
import java.io.*
import java.net.*
import java.lang.*
fprintf('done.\n');

myhost = hlp_hostname;
fprintf('This is par_filehost on %s.\n',myhost);

% determine which rank we have
matching = strcmp(settings.all_hosts,myhost);
if any(matching)
    myrank = find(matching);
    fprintf('Rank of this host is %i/%i (based on AllHosts).\n',myrank,length(matching));
    fprintf('Allhosts was: %s\n',hlp_tostring(settings.all_hosts));
elseif settings.rank_override
    myrank = settings.rank_override;
    fprintf('Rank of this host is %i (based on RankOverride).\n',myrank);
elseif isempty(settings.all_hosts)
    myrank = -1; % -1 encodes that we're operating in local mode
    fprintf('Neither AllHosts nor RankOverride given: responding only to localhost requests.\n');
else
    fprintf('Error: this host is not in AllHosts and no RankOverride was given. Exiting...\n');
    return;
end

% initialize request socket
fprintf('opening request socket...');
req_socket = DatagramSocket(settings.request_port,InetAddress.getByName('0.0.0.0')); %req_socket = DatagramSocket(settings.request_port,InetAddress.getByName('255.255.255.255'));
if myrank > 0
    req_socket.setBroadcast(true); end
req_socket.setSoTimeout(500);
closer = onCleanup(@()req_socket.close());
fprintf('done.\n');

% get own IP address
myaddr = char(InetAddress.getLocalHost());
myaddr = hlp_split(myaddr,'/');
myaddr = myaddr{2};
fprintf('listening on %s...\n',myaddr);

while true
    try
        % get next request
        fprintf('\n\nwaiting for next packet...');
        packet = DatagramPacket(uint8(zeros(1,settings.accept_packet_size)),settings.accept_packet_size);
        while true
            try
                req_socket.receive(packet);
                break;
            catch e
                if isempty(strfind(e.message,'timed out'))
                    rethrow(e); end
            end
        end
        content = char(packet.getData');
        srcport = packet.getPort();
        srcaddr = char(packet.getAddress());
        srcaddr = srcaddr(2:end);
        if myrank < 0 && ~(strcmp(srcaddr,myaddr) || strcmp(srcaddr,'127.0.0.1'))
            % not from this host while we're in local mode
            continue;
        else
            fprintf('got packet:\n%s\n\n',content);
        end

        % handle the request
        [command,args] = parse_request(content);
        fprintf('handling command %s; %i arguments.\n',command,length(fieldnames(args)));
        reply = [];
        switch command
            % --- primary commands ---
            case 'getfile'
                get_file(args);
            case 'evict'
                evict_file(args);
                % --- control commands ---
            case 'terminate'
                fprintf('  received termination command; exiting...\n');
                send_reply_udp(args,'ReplyContent',[myhost ' terminated']);
                return;
                % ---- diagnostics ---
            case 'getrank'
                reply = hlp_tostring(myrank);
            case 'allhosts'
                reply = hlp_tostring(settings.all_hosts);
            case 'hostname'
                reply = hlp_hostname;
            case 'max_capacity'
                reply = hlp_tostring(max_capacity);
            case 'clear'
                fprintf('  clearing database...');
                db = [];
                lru = {};
                fprintf('done.\n');
                reply = [myhost ' database cleared'];
            otherwise
                fprintf('  unrecognized command.\n');
        end
        % handle misc replies
        if ~isempty(reply)
            send_reply_udp(args,'ReplyContent',reply); end
    catch e
        fprintf('error: %s',hlp_handleerror(e));
        pause(2);
    end
end

    % === command handlers ===

    function get_file(varargin)
        import java.net.*
        import java.io.*
        % retrieve a file either from memory cache or from disk (if necessary)
        opts = hlp_varargin2struct(varargin, ...
            {'filename','FileName'},'', ...
            {'force','Force'},false, ...
            {'return_address','ReturnAddress'},'', ...
            {'return_port','ReturnPort'},'');

        % check if we're responsible, etc.
        [tag,namehash,responsible_host] = db_tag_if_responsible(opts.filename);
        try
            responsible_hostname = settings.all_hosts{responsible_host};
        catch
            responsible_hostname = 'N/A';
        end
        fprintf('Handling get_file request for file "%s" (hash: %s), responsible_host=%i (%s)...\n',opts.filename,namehash,responsible_host,responsible_hostname);
        if isempty(tag)
            fprintf('not responsible; ignoring.\n');
        else
            % got a request that we need to handle
            result = [];

            % check if the file matches our accept pattern
            pattern = ['^',strrep(strrep(settings.accept_pattern,'?','.'),'*','.{0,}'),'$'];
            if isempty(regexp(opts.filename,pattern, 'once'))
                fprintf('file does not match accept pattern; ignoring.\n');
            else
                if isfield(db, tag)
                    % record in cache: return it
                    fprintf('found record in cache.\n');
                    result = db.(tag);
                    lru = [{tag} lru(~strcmp(lru,tag))];
                else
                    % record not in cache
                    is_expression = any(opts.filename=='(') && any(opts.filename==')');
                    if ~is_expression && ~exist(opts.filename,'file')
                        fprintf('file does not exist.\n');
                    else
                        % check if it would fit in our DB, calculate itemsize and optionally result
                        capacity = settings.max_capacity*2^30;
                        if is_expression
                            fprintf('item appears to be a MATLAB expression, evaluating...');
                            % this is actually not a file name but an expression that we can
                            % evaluate
                            try
                                result = eval(opts.filename);
                                fprintf('success.\n');
                            catch e
                                fprintf('could not evaluate expression: %s',hlp_handleerror(e));
                                result = e;
                            end
                            if all(isfield(result,{'head','parts'}))
                                try
                                    fprintf('result is a BCILAB expression, evaluating...');
                                    result = exp_eval(result);
                                    fprintf('success.\n');
                                catch e
                                    fprintf('could not evaluate expression: %s',hlp_handleerror(e));
                                    result = e;
                                end
                            end
                            try
                                t0=tic; fprintf('serializing data structure into bytes...');
                                result = hlp_serialize(result);
                                fprintf('success (%.1f seconds).\n',toc(t0));
                            catch e
                                fprintf('could not serialize data structure: %s',hlp_handleerror(e));
                                result = hlp_serialize(e);
                            end
                            itemsize = getfield(whos('result'),'bytes');
                        else
                            % this is a file, check the size
                            itemsize = getfield(dir(opts.filename),'bytes');
                        end
                        dbsize = getfield(whos('db'),'bytes');
                        fprintf('\ndb parameters: capacity:%.3fGB, itemsize:%.3fGB, dbsize:%.3fGB, records:%f\n',capacity/2^30,itemsize/2^30,dbsize/2^30,length(lru));
                        % determine if this fits in our cache...
                        if (strcmp(settings.evict_policy,'lru') && (itemsize < capacity)) || (strcmp(settings.evict_policy,'none') && (dbsize+itemsize < capacity)) || opts.force
                            if isempty(result)
                                % don't have the result yet
                                t0=tic; fprintf('loading file...');
                                f = fopen(opts.filename,'r');
                                filecloser = onCleanup(@()fclose(f));
                                result = fread(f,inf,'*uint8');
                                fprintf('%.1f seconds (read bandwidth: %.2fMB/s).\n',toc(t0),itemsize/2^20/toc(t0));
                            else
                                % (have the result already, no need to load)
                            end
                            % cache it if applicable
                            if strcmp(settings.evict_policy,'none') && (dbsize+itemsize < capacity)
                                fprintf('committing to db.\n');
                                db.(tag) = result;
                                lru = [{tag} lru(~strcmp(lru,tag))];
                            elseif strcmp(settings.evict_policy,'lru') && (itemsize < capacity)
                                % evict least-recently used records until we're within the capacity limits
                                fprintf('committing to db.\n');
                                db.(tag) = result;
                                lru = [{tag} lru(~strcmp(lru,tag))];
                                if itemsize+dbsize > capacity
                                    fprintf('evicting least-recently used cache records');
                                    initialsize = getfield(whos('db'),'bytes');
                                    while ~isempty(fieldnames(db))
                                        fprintf('.');
                                        db = rmfield(db,lru(end));
                                        lru(end) = [];
                                        if getfield(whos('db'),'bytes') < capacity
                                            break; end
                                    end
                                    fprintf('done. Evicted %.3fGB of data.\n',(initialsize - getfield(whos('db'),'bytes'))/2^30);
                                end
                            end
                        else
                            fprintf('does not fit in allowed cache capacity; not loading.\n');
                        end
                    end
                end
            end

            % send off the result
            if isa(result,'uint8')
                [return_address,return_port] = get_return_endpoint(opts);
                fprintf('preparing to return result to requestor at %s:%i\n',char(return_address.toString()),return_port);
                % send the data back
                fprintf('opening back-connection to transmit data...');
                outconn = Socket();
                destination = InetSocketAddress(return_address, return_port);
                outconn.connect(destination, round(settings.timeout_backconnect));
                sockcloser = onCleanup(@()outconn.close());
                outconn.setSoTimeout(round(settings.timeout_transmit));
                outconn.setSoLinger(true,3);
                fprintf('connected.\nnow sending (%.0fkb)...',length(result)/1024);
                t0 = tic;
                out = DataOutputStream(outconn.getOutputStream());
                out.writeDouble(length(result));
                % we need to efficiently split the data into chunks no larger than max_java_memory
                num_blocks = ceil(length(result)/max_java_memory);
                rest_size = mod(length(result),max_java_memory);
                sizes = [max_java_memory*ones(1,num_blocks-1) rest_size];
                blocks = cell(length(sizes),1);
                [blocks{:}] = chopdeal(result,sizes);
                % send the blocks off
                for i=1:length(blocks)
                    out.write(blocks{i},0,sizes(i)); end
                out.flush();
                fprintf('done (%.1f seconds; %.2fMB/s bandwidth)\nnow closing...',toc(t0),length(result)/2^20/toc(t0));
                fprintf('done.\ntransaction complete.\n');
            end
        end
    end

    function evict_file(varargin)
        % evict a file record from cache
        opts = hlp_varargin2struct(varargin, ...
            {'filename','FileName'},'', ...
            {'anyhost','AnyHost'},false);
        [tag,namehash,responsible_host] = db_tag_if_responsible(opts.filename, opts.anyhost);
        fprintf('Handling evict request for file "%s" (hash: %s), responsible_host=%i (%s)...',opts.filename,namehash,responsible_host,settings.all_hosts{responsible_host});
        if isfield(db,tag)
            try
                db = rmfield(db,tag);
                lru(strcmp(lru,tag)) = [];
                fprintf('evicted record.\n');
            catch e
                fprintf('could not evict record (%s); database may now be inconsistent.\n', e.message);
            end
        elseif ~isempty(tag)
            fprintf('no record in database.\n');
        else
            fprintf('no record in database or not responsible.\n');
        end
    end

    % === helper functions ===

    function send_reply_udp(varargin)
        import java.net.*
        % send a reply message over UDP
        opts = hlp_varargin2struct(varargin, ...
            {'reply_content','ReplyContent'},'');
        [return_address,return_port] = get_return_endpoint(opts);
        fprintf('sending reply packet to %s:%i...',char(return_address.toString()),return_port);
        repl_sock = DatagramSocket();
        repl_packet = DatagramPacket(uint8(opts.reply_content),length(opts.reply_content),return_address,return_port);
        repl_sock.send(repl_packet)
        fprintf('done.\n')
    end

    function [return_address, return_port] = get_return_endpoint(varargin)
        import java.net.*
        % determine the return endpoint based on the global srcaddr/port or override options
        opts = hlp_varargin2struct(varargin, ...
            {'return_address','ReturnAddress'},'', ...
            {'return_port','ReturnPort'},'');
        % allow overriding the return address / port
        fprintf('determining return address/port...');
        if isfield(opts,'return_address') && ~isempty(opts.return_address)
            return_address = InetAddress.getByName(opts.return_address);
        else
            return_address = InetAddress.getByName(srcaddr);
        end
        if isfield(opts,'return_port') && ~isempty(opts.return_port)
            if ischar(opts.return_port)
                return_port = str2num(opts.return_port);
            else
                return_port = opts.return_port;
            end
        else
            return_port = srcport;
        end
        fprintf('done (%s:%i).\n', char(return_address), return_port);
    end

    function [tag,namehash,responsible_host] = db_tag_if_responsible(filename, anyhost)
        % returns a database tag if we're responsible for the given filename or empty otherwise
        % get a hash of the filename
        if nargin < 2
            anyhost = false; end
        namehash = hlp_cryptohash(filename);
        % check if we're responsible for this file
        responsible_host = 1+mod(hex2dec(namehash(1:5))+1,length(settings.all_hosts));
        if anyhost || myrank < 0 || responsible_host == myrank
            tag = ['x_' namehash];
        else
            tag = [];
        end

    end

    function [command,args] = parse_request(content)
        % parses HTTP-style requests (see top of this file for format) into command string and
        % argument struct
        content(content==0) = [];
        % parse a network request packet
        lines = hlp_split(content,sprintf('\n'));
        % first line is the command
        command = strtrim(lower(lines{1}));
        % followed by one argument per line
        args = [];
        for l=2:length(lines)
            line = strtrim(lines{l});
            % strip comments (beginning with ;)
            comment = find(line == ';',1);
            if ~isempty(comment)
                line = line(1:comment-1); end
            if isempty(line)
                continue; end
            % handle assignments
            equals = find(line == '=',1);
            if isempty(equals)
                fprintf('the line "%s" is not a valid argument declaration; ignoring.\n',line); end
            key = strtrim(line(1:equals-1));
            val = strtrim(line(equals+1:end));
            try
                args.(key) = val;
            catch
                fprintf('argument "%s" is not a valid field name; ignoring.\n',key);
            end
        end
    end
end

repo: REST


text:
function undersamplingRINEX(filenameIN, filenameOUT, base, step, delta, wait_dlg)

% SYNTAX:
%   undersamplingRINEX(filenameIN, filenameOUT, base, step, delta, wait_dlg);
%
% INPUT:
%   filenameIN  = input RINEX observation file
%   filenameOUT = output RINEX observation file
%   base  = base timing (e.g. if first epoch should be 13  4  5 11 35  1, then base = 1) [sec]
%   step  = new sampling rate [sec]
%   delta = original sampling rate [sec]
%   wait_dlg = optional handler to waitbar figure (optional)
%
% OUTPUT:
%
%
% DESCRIPTION:
%   Undersamples a RINEX observation file.

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

%default values
if (nargin < 3)
    base  = 0;
end

if (nargin < 4)
    step  = 30;
end

if (nargin < 5)
    delta = 1;
end

% ------------------------------------------
% header copy and paste
% ------------------------------------------

fidIN  = fopen(filenameIN,'rt');
fidOUT = fopen(filenameOUT,'wt');

line = fgets(fidIN);
while (isempty(strfind(line,'END OF HEADER')))
   fprintf(fidOUT,'%s',line);
   line = fgets(fidIN);
   if (~isempty(strfind(line,'INTERVAL')))
       line = sprintf('%10.3f                                                  INTERVAL            \n', step);
   end
end
fprintf(fidOUT,line);

% fclose(fidIN);
% fclose(fidOUT);

% ------------------------------------------
% undersampling
% ------------------------------------------

while (feof(fidIN) == 0)

    % time (seconds) extraction
	sec = roundmod(str2num(line(17:27)),delta);
	%fprintf('%d\n',sec);           % debugging

	if (mod(sec-base,step) == 0)    % if it is a header message
		fprintf(fidOUT,'%s',line);
		line = fgets(fidIN);
		%fprintf('%s',line);        % debugging
		while (feof(fidIN) == 0) & ...
		      ((line(1) ~= ' ') | (line(2) == ' ') | ...
		       (line(3) == ' ') | (line(4) ~= ' '))
			fprintf(fidOUT,'%s',line);
			line = fgets(fidIN);
			%fprintf('%s',line);    % debugging
        end
        if (feof(fidIN) == 1) && ~isempty(line)
            fprintf(fidOUT,'%s',line);
        end
	else                            % if it is a data message
		line = fgets(fidIN);
		%fprintf('%s',line);        % debugging
		while (feof(fidIN) == 0) & ...
		      ((line(1) ~= ' ') | (line(2) == ' ') | ...
		       (line(3) == ' ') | (line(4) ~= ' '))
			line = fgets(fidIN);
			%fprintf('%s',line);    % debugging
		end
	end
end

fclose(fidIN);      % input file closure
fclose(fidOUT);     % output file closure

author: goGPS-Project


text:
classdef TestConvertTo
    %TestConvertTo

    methods (Static)
        function test_1
            A = magic(5);
            B = cv.convertTo(A, 'RType','single', 'Alpha',2, 'Beta',3);
            validateattributes(B, {'single'}, {'size',size(A)});
            assert(isequal(B, single(A*2+3)));
        end

        function test_2
            img = imread(fullfile(mexopencv.root(),'test','cat.jpg'));
            dst = cv.convertTo(img, 'RType','double', 'Alpha',1/255);
            validateattributes(dst, {'double'}, {'size',size(img)});
            assert(norm(dst(:) - double(img(:))/255) < 1e-6);
        end

        function test_error_argnum
            try
                cv.convertTo();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function I = mat2cell_ov(X,xx_s,xx_f,yy_s,yy_f,zz_s,zz_f,overlap,sz)

% converts a matrix into a cell array with overlapping elements
% INPUTS:
% X:            Input matrix
% grid_size:    size of each element without overlap
% overlap:      amount of overlap
% sz:           spatial size of X

% OUTPUT:
% I:            output cell array

% Written by Eftychios A. Pnevmatikakis, Simons Foundation, 2016

I = cell(length(xx_s),length(yy_s),length(zz_s));
nd = length(sz);
if nd == 2; sz(3) = 1; end
for i = 1:length(xx_s)
    for j = 1:length(yy_s)
        for k = 1:length(zz_s)
            extended_grid = [max(xx_s(i)-overlap(1),1),min(xx_f(i)+overlap(1),sz(1)),max(yy_s(j)-overlap(2),1),min(yy_f(j)+overlap(2),sz(2)),max(zz_s(k)-overlap(3),1),min(zz_f(k)+overlap(3),sz(3))];
            if nd == 2
                I{i,j} = X(extended_grid(1):extended_grid(2),extended_grid(3):extended_grid(4),:);
            else
                I{i,j,k} = X(extended_grid(1):extended_grid(2),extended_grid(3):extended_grid(4),extended_grid(5):extended_grid(6),:);
            end
        end
    end
end

author: flatironinstitute

repo: NoRMCorre


text:
%% Copyright (C) 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym gammaln (@var{x})
%% @defmethodx @@sym lgamma (@var{x})
%% Symbolic logarithm of the gamma function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = gammaln(x)
%%   @result{} y = (sym) loggamma(x)
%% y = lgamma(x)
%%   @result{} y = (sym) loggamma(x)
%% @end group
%% @end example
%%
%% @seealso{gammaln, @@sym/gamma, @@sym/psi}
%% @end defmethod

function y = lgamma(x)
  y = gammaln(x);
end

%!test
%! % tested by gammaln
%! assert (isequal (lgamma (sym ('x')), gammaln (sym ('x'))))

repo: octsympy


text:
function [rgb,points3d,points3dMatrix]=read_3d_pts_general(depthInpaint,K,depthInpaintsize,imageName,crop)
    %K is [fx 0 cx; 0 fy cy; 0 0 1];
    %K = frames.K;
    if ~isempty(K)
        cx = K(1,3); cy = K(2,3);
        fx = K(1,1); fy = K(2,2);
    else
        fx = 5.1885790117450188e+02;
        fy = 5.1946961112127485e+02;
        cx = 3.2558244941119034e+02;
        cy = 2.5373616633400465e+02;
    end
    invalid = depthInpaint==0;
    if ~isempty(imageName)
        im = imread(imageName);
        rgb = im2double(im);
    else
        rgb =double(cat(3,zeros(depthInpaintsize(1),depthInpaintsize(2)),...
                    ones(depthInpaintsize(1),depthInpaintsize(2)),...
                    zeros(depthInpaintsize(1),depthInpaintsize(2))));
    end
    rgb = reshape(rgb, [], 3);
    %3D points
    [x,y] = meshgrid(crop(2)-1+(1:depthInpaintsize(2)), crop(1)-1+(1:depthInpaintsize(1)));
    x3 = (x-cx).*depthInpaint*1/fx;
    y3 = (y-cy).*depthInpaint*1/fy;
    z3 = depthInpaint;
    points3dMatrix =cat(3,x3,z3,-y3);
    points3dMatrix(cat(3,invalid,invalid,invalid))=NaN;
    points3d = [x3(:) z3(:) -y3(:)];
    points3d(invalid(:),:) =NaN;
end

author: thusiyuan


text:
% std_precomp() - Precompute measures (ERP, spectrum, ERSP, ITC) for channels in a
%                 study. If channels are interpolated before computing the measures,
%                 the updated EEG datasets are also saved to disk. Called by
%                 pop_precomp(). Follow with pop_plotstudy(). See Example below.
% Usage:
% >> [STUDY ALLEEG customRes] = std_precomp(STUDY, ALLEEG, chanorcomp, 'key', 'val', ...);
%
% Required inputs:
%   STUDY        - an EEGLAB STUDY set of loaded EEG structures
%   ALLEEG       - ALLEEG vector of one or more loaded EEG dataset structures
%   chanorcomp   - ['components'|'channels'| or channel cell array] The string
%                  'components' forces the program to precompute all selected
%                  measures for components. The string 'channels' forces the
%                  program to compute all measures for all channels.
%                  A channel cell array containing channel labels will precompute
%                  the selected measures. Note that the name of the channel is
%                  not case-sensitive.
% Optional inputs:
%  'design'   - [integer] use specific study index design to compute measure.
%  'cell'     - [integer] compute measure only for a give data file.
%  'erp'      - ['on'|'off'] pre-compute ERPs for each dataset.
%  'spec'     - ['on'|'off'] pre-compute spectrum for each dataset.
%               Use 'specparams' to set spectrum parameters.
%  'ersp'     - ['on'|'off'] pre-compute ERSP for each dataset.
%               Use 'erspparams' to set time/frequency parameters.
%  'itc'      - ['on'|'off'] pre-compute ITC for each dataset.
%               Use 'erspparams' to set time/frequency parameters.
%  'scalp'    - ['on'|'off'] pre-compute scalp maps for components.
%  'allcomps' - ['on'|'off'] compute ERSP/ITC for all components ('off'
%               only use pre-selected components in the pop_study interface).
%  'erpparams'   - [cell array] Parameters for the std_erp function. See
%                  std_erp for more information.
%  'specparams'  - [cell array] Parameters for the std_spec function. See
%                  std_spec for more information.
%  'erspparams'  - [cell array] Optional arguments for the std_ersp function.
%  'erpimparams' - [cell array] Optional argument for std_erpimage. See
%                  std_erpimage for the list of arguments.
%  'recompute'   - ['on'|'off'] force recomputing ERP file even if it is
%                  already on disk.
%  'rmicacomps'  - ['on'|'off'|'processica'] remove ICA components pre-selected in
%                  each dataset (EEGLAB menu item, "Tools > Reject data using ICA
%                  > Reject components by map). This option is ignored when
%                  precomputing measures for ICA clusters. Default is 'off'.
%                  'processica' forces to process ICA components instead of
%                  removing them.
%  'rmclust'     - [integer array] remove selected ICA component clusters.
%                  For example, ICA component clusters containing
%                  artifacts. This option is ignored when precomputing
%                  measures for ICA clusters.
%  'savetrials'  - ['on'|'off'] save single-trials ERSP. Requires a lot of disk
%                  space (dataset space on disk times 10) but allow for refined
%                  single-trial statistics.
%  'customfunc'  - [function_handle] execute a specific function on each
%                  EEGLAB dataset of the selected STUDY design. The fist
%                  argument to the function is an EEGLAB dataset. Example is
%                  @(EEG)mean(EEG.data,3)
%                  This will compute the ERP for the STUDY design. EEG is the
%                  EEGLAB dataset corresponding to each cell design. It
%                  corresponds to a dataset computed dynamically based on
%                  the design selection. If 'rmclust', 'rmicacomps' or 'interp'
%                  are being used, the channel data is affected
%                  accordingly. Anonymous and non-anonymous functions may be
%                  used. The output is returned in CustomRes or saved on
%                  disk. The output of the custom function may be an numerical
%                  array or a structure.
%  'customparams' - [cell array] Parameters for the custom function above.
%  'customfileext' - [string] file extension for saving custom data. Use
%                    function to read custom data. If left empty, the
%                    result is returned in the customRes output. Note that
%                    if the custom function does not return a structure,
%                    the data is automatically saved in a variable named
%                    'data'.
%  'customclusters' - [integer array] load only specific clusters. This is
%                    used with SIFT. chanorcomp 3rd input must be 'components'.
%
% Outputs:
%   ALLEEG       - the input ALLEEG vector of EEG dataset structures, modified
%                  by adding preprocessing data as pointers to Matlab files that
%                  hold the pre-clustering component measures.
%   STUDY        - the input STUDY set with pre-clustering data added,
%                  for use by pop_clust()
%   customRes    - cell array of custom results (one cell for each pair of
%                  independent variables as defined in the STUDY design).
%                  If a custom file extension is specified, this variable
%                  is empty as the function assumes that the result is too
%                  large to hold in memory.
%
% Example:
%   >> [STUDY ALLEEG customRes] = std_precomp(STUDY, ALLEEG, { 'cz' 'oz' }, 'interp', ...
%               'on', 'erp', 'on', 'spec', 'on', 'ersp', 'on', 'erspparams', ...
%               { 'cycles' [ 3 0.5 ], 'alpha', 0.01, 'padratio' 1 });
%
%           % This prepares, channels 'cz' and 'oz' in the STUDY datasets.
%           % If a data channel is missing in one dataset, it will be
%           % interpolated (see eeg_interp()). The ERP, spectrum, ERSP, and
%           % ITC for each dataset is then computed.
%
% Example of custom call:
%   The function below computes the ERP of the EEG data for each channel and plots it.
%   >> [STUDY ALLEEG customres] = std_precomp(STUDY, ALLEEG, 'channels', 'customfunc', @(EEG,varargin)(mean(EEG.data,3)'));
%   >> std_plotcurve([1:size(customres{1},1)], customres, 'chanlocs', eeg_mergelocs(ALLEEG.chanlocs)); % plot data
%
%   The function below uses a data file to store the information then read
%   the data and eventyally plot it
%   >> [STUDY ALLEEG customres] = std_precomp(STUDY, ALLEEG, 'channels', 'customfunc', @(EEG,varargin)(mean(EEG.data,3)), 'customfileext', 'tmperp');
%   >> erpdata = std_readcustom(STUDY, ALLEEG, 'tmperp');
%   >> std_plotcurve([1:size(erpdata{1})], erpdata, 'chanlocs', eeg_mergelocs(ALLEEG.chanlocs)); % plot data
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, 2006-

% Copyright (C) Arnaud Delorme, SCCN, INC, UCSD, 2006, arno@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [ STUDY, ALLEEG customRes ] = std_precomp(STUDY, ALLEEG, chanlist, varargin)

    if nargin < 2
        help std_precomp;
        return;
    end;

    if nargin == 2
        chanlist = 'channels'; % default to clustering the whole STUDY
    end
    customRes = [];
    Ncond = length(STUDY.condition);
    if Ncond == 0
        Ncond = 1;
    end

    g = finputcheck(varargin, { 'erp'         'string'  { 'on','off' }     'off';
                                'interp'      'string'  { 'on','off' }     'off';
                                'ersp'        'string'  { 'on','off' }     'off';
                                'recompute'   'string'  { 'on','off' }     'off';
                                'spec'        'string'  { 'on','off' }     'off';
                                'erpim'       'string'  { 'on','off' }     'off';
                                'scalp'       'string'  { 'on','off' }     'off';
                                'allcomps'    'string'  { 'on','off' }     'off';
                                'itc'         'string'  { 'on','off' }     'off';
                                'savetrials'  'string'  { 'on','off' }     'off';
                                'rmicacomps'  'string'  { 'on','off','processica' }     'off';
                                'cell'        'integer' []                 [];
                                'design'      'integer' []                 STUDY.currentdesign;
                                'rmclust'     'integer' []                 [];
                                'rmbase'      'integer' []                 []; % deprecated, for backward compatibility purposes, not documented
                                'specparams'        'cell'    {}                 {};
                                'erpparams'         'cell'    {}                 {};
                                'customfunc'  {'function_handle' 'integer' } { { } {} }     [];
                                'customparams'      'cell'    {}                 {};
                                'customfileext'     'string'  []                 '';
                                'customclusters'    'integer' []                 [];
                                'erpimparams'       'cell'    {}                 {};
                                'erspparams'        'cell'    {}                 {}}, 'std_precomp');
    if isstr(g), error(g); end;
    if ~isempty(g.rmbase), g.erpparams = { g.erpparams{:} 'rmbase' g.rmbase }; end;

    % union of all channel structures
    % -------------------------------
    computewhat = 'channels';
    if isstr(chanlist)
        if strcmpi(chanlist, 'channels')
            chanlist = [];
        else % components
            computewhat = 'components';
            if strcmpi(g.allcomps, 'on')
                chanlist = {};
                for index = 1:length(STUDY.datasetinfo)
                    chanlist = { chanlist{:} [1:size(ALLEEG(STUDY.datasetinfo(index).index).icaweights,1)] };
                end;
            else
                chanlist = { STUDY.datasetinfo.comps };
            end;
        end;
    end;
    if isempty(chanlist)
        alllocs = eeg_mergelocs(ALLEEG.chanlocs);
        chanlist = { alllocs.labels };
    elseif ~isnumeric(chanlist{1})
        alllocs = eeg_mergelocs(ALLEEG.chanlocs);
        [tmp c1 c2] = intersect_bc( lower({ alllocs.labels }), lower(chanlist));
        [tmp c2] = sort(c2);
        alllocs = alllocs(c1(c2));
    end;

    % test if interp and reconstruct channel list
    % -------------------------------------------
    if strcmpi(computewhat, 'channels')
        if strcmpi(g.interp, 'on')
            STUDY.changrp = [];
            STUDY = std_changroup(STUDY, ALLEEG, chanlist, 'interp');
            g.interplocs = alllocs;
        else
            STUDY.changrp = [];
            STUDY = std_changroup(STUDY, ALLEEG, chanlist);
            g.interplocs = struct([]);
        end;
    end;

    % components or channels
    % ----------------------
    if strcmpi(computewhat, 'channels')
         curstruct = STUDY.changrp;
    else curstruct = STUDY.cluster;
    end;

    % compute custom measure
    % ----------------------
    if ~isempty(g.customfunc)
        nc = max(length(STUDY.design(g.design).variable(1).value),1);
        ng = max(length(STUDY.design(g.design).variable(2).value),1);
        allinds = curstruct(1).allinds; % same for all channels and components (see std_selectdesign)
        setinds = curstruct(1).setinds; % same for all channels and components (see std_selectdesign)
        if ~isempty(g.customclusters)
            allinds = curstruct(g.customclusters).allinds; % same for all channels and components (see std_selectdesign)
            setinds = curstruct(g.customclusters).setinds; % same for all channels and components (see std_selectdesign)
        end;

        for cInd = 1:nc
            for gInd = 1:ng
                if ~isempty(setinds{cInd,gInd})
                    desset = STUDY.design(g.design).cell(setinds{cInd,gInd}(:));
                    for iDes = 1:length(desset)
                        if strcmpi(computewhat, 'channels')
                             [tmpchanlist opts] = getchansandopts(STUDY, ALLEEG, chanlist, desset(iDes).dataset, g);
                             TMPEEG = std_getdataset(STUDY, ALLEEG, 'design', g.design, 'cell', setinds{cInd,gInd}(iDes), opts{:}); % trial indices included in cell selection
                        else TMPEEG = std_getdataset(STUDY, ALLEEG, 'design', g.design, 'cell', setinds{cInd,gInd}(iDes), 'cluster', g.customclusters);
                        end;
                        addopts = { 'savetrials', g.savetrials, 'recompute', g.recompute }; % not currently used

                        tmpData = feval(g.customfunc, TMPEEG, g.customparams{:});
                        if isempty(g.customfileext)
                            resTmp(:,:,iDes) = tmpData;
                        else
                            fileName = [ desset(iDes).filebase '.' g.customfileext ];
                            clear data;
                            data.data = tmpData;
                            data.datafile   = computeFullFileName( { ALLEEG(desset(iDes).dataset).filepath }, { ALLEEG(desset(iDes).dataset).filename });
                            data.datatrials = desset(iDes).trials;
                            data.datatype = upper(g.customfileext);
                            if ~isempty(g.customparams) data.parameters = g.customparams; end;
                            std_savedat(fileName, data);
                        end;
                    end;
                    if isempty(g.customfileext)
                        customRes{cInd,gInd} = resTmp;
                    end;
                    clear resTmp;
                end;
            end;
        end;
    end;

    % compute ERPs
    % ------------
    if strcmpi(g.erp, 'on')
        % check dataset consistency
        % -------------------------
        allPnts = [ALLEEG([STUDY.design(g.design).cell.dataset]).pnts];
        if iscell(allPnts), allPnts = [ allPnts{:} ]; end;
        if length(unique(allPnts)) > 1
            error([ 'Cannot compute ERPs because datasets' 10 'do not have the same number of data points' ])
        end;

        for index = 1:length(STUDY.design(g.design).cell)
            if ~isempty(g.cell)
                 desset = STUDY.design(g.design).cell(g.cell);
            else desset = STUDY.design(g.design).cell(index);
            end;
            addopts = { 'savetrials', g.savetrials, 'recompute', g.recompute, 'fileout', desset.filebase, 'trialindices', desset.trials };
            if strcmpi(computewhat, 'channels')
                [tmpchanlist opts] = getchansandopts(STUDY, ALLEEG, chanlist, desset.dataset, g);
                std_erp(ALLEEG(desset.dataset), 'channels', tmpchanlist, opts{:}, addopts{:}, g.erpparams{:});
            else
                if length(desset.dataset)>1 && ~isequal(chanlist{desset.dataset})
                    error(['ICA decompositions must be identical if' 10 'several datasets are concatenated to build' 10 'the design, abording' ]);
                end;
                std_erp(ALLEEG(desset.dataset), 'components', chanlist{desset.dataset(1)}, addopts{:}, g.erpparams{:});
            end;
            if ~isempty(g.cell), break; end;
        end;
        if isfield(curstruct, 'erpdata')
            curstruct = rmfield(curstruct, 'erpdata');
            curstruct = rmfield(curstruct, 'erptimes');
        end;
    end;

    % compute spectrum
    % ----------------
    if strcmpi(g.spec, 'on')
        % check dataset consistency
        % -------------------------
        for index = 1:length(STUDY.design(g.design).cell)
            if ~isempty(g.cell)
                 desset = STUDY.design(g.design).cell(g.cell);
            else desset = STUDY.design(g.design).cell(index);
            end;
            addopts = { 'savetrials', g.savetrials, 'recompute', g.recompute, 'fileout', desset.filebase, 'trialindices', desset.trials };
            if strcmpi(computewhat, 'channels')
                [tmpchanlist opts] = getchansandopts(STUDY, ALLEEG, chanlist, desset.dataset, g);
                std_spec(ALLEEG(desset.dataset), 'channels', tmpchanlist, opts{:}, addopts{:}, g.specparams{:});
            else
                if length(desset.dataset)>1 && ~isequal(chanlist{desset.dataset})
                    error(['ICA decompositions must be identical if' 10 'several datasets are concatenated to build' 10 'the design, abording' ]);
                end;
                std_spec(ALLEEG(desset.dataset), 'components', chanlist{desset.dataset(1)}, addopts{:}, g.specparams{:});
            end;
            if ~isempty(g.cell), break; end;
        end;
        if isfield(curstruct, 'specdata')
            curstruct = rmfield(curstruct, 'specdata');
            curstruct = rmfield(curstruct, 'specfreqs');
        end;
    end;

    % compute spectrum
    % ----------------
    if strcmpi(g.erpim, 'on')
        % check dataset consistency
        % -------------------------
        allPnts = [ALLEEG([STUDY.design(g.design).cell.dataset]).pnts];
        if iscell(allPnts), allPnts = [ allPnts{:} ]; end;
        if length(unique(allPnts)) > 1
            error([ 'Cannot compute ERSPs/ITCs because datasets' 10 'do not have the same number of data points' ])
        end;

        % check consistency with parameters on disk
        % -----------------------------------------
        guimode = 'guion';
        tmpparams = {};
        tmpparams = g.erpimparams;
        if strcmpi(g.recompute, 'off')
            for index = 1:length(STUDY.design(g.design).cell)
                desset = STUDY.design(g.design).cell(index);
                if strcmpi(computewhat, 'channels')
                     filename = [ desset.filebase '.daterpim'];
                else filename = [ desset.filebase '.icaerpim'];
                end;
                [guimode, g.erpimparams] = std_filecheck(filename, g.erpimparams, guimode, { 'fileout' 'recompute', 'channels', 'components', 'trialindices'});
                if strcmpi(guimode, 'cancel'), return; end;
            end;
            if strcmpi(guimode, 'usedisk') || strcmpi(guimode, 'same'), g.recompute = 'off';
            else                                                        g.recompute = 'on';
            end;
            if ~isempty(g.erpimparams) && isstruct(g.erpimparams)
                tmpparams      = fieldnames(g.erpimparams); tmpparams = tmpparams';
                tmpparams(2,:) = struct2cell(g.erpimparams);
            end;
        end;

        % set parameters in ERPimage parameters
        % -------------------------------------
        STUDY = pop_erpimparams(STUDY, tmpparams{:}); % a little trashy as the function pop_erpimparams does not check the fields

        % compute ERPimages
        % -----------------
        for index = 1:length(STUDY.design(g.design).cell)
            if ~isempty(g.cell)
                 desset = STUDY.design(g.design).cell(g.cell);
            else desset = STUDY.design(g.design).cell(index);
            end;
            addopts = { 'recompute', g.recompute, 'fileout', desset.filebase, 'trialindices', desset.trials };
            if strcmpi(computewhat, 'channels')
                [tmpchanlist opts] = getchansandopts(STUDY, ALLEEG, chanlist, desset.dataset, g);
                std_erpimage(ALLEEG(desset.dataset), 'channels', tmpchanlist, opts{:}, addopts{:}, tmpparams{:});
            else
                if length(desset.dataset)>1 && ~isequal(chanlist{desset.dataset})
                    error(['ICA decompositions must be identical if' 10 'several datasets are concatenated to build' 10 'the design, abording' ]);
                end;
                std_erpimage(ALLEEG(desset.dataset), 'components', chanlist{desset.dataset(1)}, addopts{:}, tmpparams{:});
            end;
            if ~isempty(g.cell), break; end;
        end;
        if isfield(curstruct, 'erpimdata')
            curstruct = rmfield(curstruct, 'erpimdata');
            curstruct = rmfield(curstruct, 'erpimtimes');
            curstruct = rmfield(curstruct, 'erpimtrials');
            curstruct = rmfield(curstruct, 'erpimevents');
        end;
    end;

    % compute component scalp maps
    % ----------------------------
    if strcmpi(g.scalp, 'on') && ~strcmpi(computewhat, 'channels')
        for index = 1:length(STUDY.datasetinfo)

            % find duplicate
            % --------------
            found = [];
            ind1 = STUDY.datasetinfo(index).index;
            inds = strmatch(STUDY.datasetinfo(index).subject, { STUDY.datasetinfo(1:index-1).subject });
            for index2 = 1:length(inds)
                ind2 = STUDY.datasetinfo(inds(index2)).index;
                if isequal(ALLEEG(ind1).icawinv, ALLEEG(ind2).icawinv)
                    found = ind2;
                end;
            end;

            % make link if duplicate
            % ----------------------
            if ~isempty(g.cell)
                desset = STUDY.design(g.design).cell(g.cell);
                [path,tmp] = fileparts(desset.filebase);
            else path = ALLEEG(index).filepath;
            end;

            fprintf('Computing/checking topo file for dataset %d\n', ind1);
            if ~isempty(found)
                clear tmp;
                tmpfile1 = fullfile( path, [ ALLEEG(index).filename(1:end-3) 'icatopo' ]);
                tmp.file = fullfile( ALLEEG(found).filepath, [ ALLEEG(found).filename(1:end-3) 'icatopo' ]);
                std_savedat(tmpfile1, tmp);
            else
                std_topo(ALLEEG(index), chanlist{index}, 'none', 'recompute', g.recompute,'fileout',path);
            end;
        end;
        if isfield(curstruct, 'topo')
            curstruct = rmfield(curstruct, 'topo');
            curstruct = rmfield(curstruct, 'topox');
            curstruct = rmfield(curstruct, 'topoy');
            curstruct = rmfield(curstruct, 'topoall');
            curstruct = rmfield(curstruct, 'topopol');
        end;
    end;

    % compute ERSP and ITC
    % --------------------
    if strcmpi(g.ersp, 'on') || strcmpi(g.itc, 'on')
        % check dataset consistency
        % -------------------------
        allPnts = [ALLEEG([STUDY.design(g.design).cell.dataset]).pnts];
        if iscell(allPnts), allPnts = [ allPnts{:} ]; end;
        if length(unique(allPnts)) > 1
            error([ 'Cannot compute ERSPs/ITCs because datasets' 10 'do not have the same number of data points' ])
        end;

        if strcmpi(g.ersp, 'on') & strcmpi(g.itc, 'on'), type = 'both';
        elseif strcmpi(g.ersp, 'on')                   , type = 'ersp';
        else                                             type = 'itc';
        end;

        % check for existing files
        % ------------------------
        guimode = 'guion';
        [ tmpX tmpt tmpf g.erspparams ] = std_ersp(ALLEEG(1), 'channels', 1, 'type', type, 'recompute', 'on', 'getparams', 'on', 'savetrials', g.savetrials, g.erspparams{:});
        if strcmpi(g.recompute, 'off')
            for index = 1:length(STUDY.design(g.design).cell)
                desset = STUDY.design(g.design).cell(index);
                if strcmpi(computewhat, 'channels')
                     filename = [ desset.filebase '.datersp'];
                else filename = [ desset.filebase '.icaersp'];
                end;
                [guimode, g.erspparams] = std_filecheck(filename, g.erspparams, guimode, { 'plotitc' 'plotersp' 'plotphase' });
                if strcmpi(guimode, 'cancel'), return; end;
            end;
            if strcmpi(guimode, 'usedisk') || strcmpi(guimode, 'same'), g.recompute = 'off';
            else                                                        g.recompute = 'on';
            end;
        end;

        % check for existing files
        % ------------------------
        if isempty(g.erspparams),
            tmpparams = {};
        elseif iscell(g.erspparams),
            tmpparams = g.erspparams;
        else
            tmpparams      = fieldnames(g.erspparams); tmpparams = tmpparams';
            tmpparams(2,:) = struct2cell(g.erspparams);
        end;
        tmpparams = { tmpparams{:} 'recompute' g.recompute };
        for index = 1:length(STUDY.design(g.design).cell)
            if ~isempty(g.cell)
                 desset = STUDY.design(g.design).cell(g.cell);
            else desset = STUDY.design(g.design).cell(index);
            end;
            if strcmpi(computewhat, 'channels')
                [tmpchanlist opts] = getchansandopts(STUDY, ALLEEG, chanlist, desset.dataset, g);
                std_ersp(ALLEEG(desset.dataset), 'channels', tmpchanlist, 'type', type, 'fileout', desset.filebase, 'trialindices', desset.trials, opts{:}, tmpparams{:});
            else
                if length(desset.dataset)>1 && ~isequal(chanlist{desset.dataset})
                    error(['ICA decompositions must be identical if' 10 'several datasets are concatenated to build' 10 'the design, abording' ]);
                end;
                std_ersp(ALLEEG(desset.dataset), 'components', chanlist{desset.dataset(1)}, 'type', type, 'fileout', desset.filebase, 'trialindices', desset.trials, tmpparams{:});
            end;
            if ~isempty(g.cell), break; end;
        end;
        if isfield(curstruct, 'erspdata')
            curstruct = rmfield(curstruct, 'erspdata');
            curstruct = rmfield(curstruct, 'ersptimes');
            curstruct = rmfield(curstruct, 'erspfreqs');
        end;
        if isfield(curstruct, 'itcdata')
            curstruct = rmfield(curstruct, 'itcdata');
            curstruct = rmfield(curstruct, 'itctimes');
            curstruct = rmfield(curstruct, 'itcfreqs');
        end;
    end;

    % components or channels
    % ----------------------
    if strcmpi(computewhat, 'channels')
         STUDY.changrp = curstruct;
    else STUDY.cluster = curstruct;
    end;

    return;

    % find components in cluster for specific dataset
    % -----------------------------------------------
    function rmcomps = getclustcomps(STUDY, rmclust, settmpind);

        rmcomps   = cell(1,length(settmpind));
        for idat = 1:length(settmpind) % scan dataset for which to find component clusters
            for rmi = 1:length(rmclust) % scan clusters
                comps = STUDY.cluster(rmclust(rmi)).comps;
                sets  = STUDY.cluster(rmclust(rmi)).sets;
                indmatch = find(sets(:) == settmpind(idat));
                indmatch = ceil(indmatch/size(sets,1)); % get the column number
                rmcomps{idat} = [rmcomps{idat} comps(indmatch(:)') ];
            end;
            rmcomps{idat} = sort(rmcomps{idat});
        end;

    % make option array and channel list (which depend on interp) for any type of measure
    % ----------------------------------------------------------------------
    function [tmpchanlist, opts] = getchansandopts(STUDY, ALLEEG, chanlist, idat, g);

        opts = { };

        if ~isempty(g.rmclust) || strcmpi(g.rmicacomps, 'on') || strcmpi(g.rmicacomps, 'processica')
            rmcomps = cell(1,length(idat));
            if ~isempty(g.rmclust)
                rmcomps = getclustcomps(STUDY, g.rmclust, idat);
            end;
            if strcmpi(g.rmicacomps, 'on')
                for ind = 1:length(idat)
                    rmcomps{ind} = union_bc(rmcomps{ind}, find(ALLEEG(idat(ind)).reject.gcompreject));
                end;
            elseif strcmpi(g.rmicacomps, 'processica')
                for ind = 1:length(idat)
                    rmcomps{ind} = union_bc(rmcomps{ind}, find(~ALLEEG(idat(ind)).reject.gcompreject));
                end;
            end;
            opts = { opts{:} 'rmcomps' rmcomps };
        end;
        if strcmpi(g.interp, 'on')
            tmpchanlist = chanlist;
            allocs = eeg_mergelocs(ALLEEG.chanlocs);
            [tmp1 tmp2 neworder] = intersect_bc( {allocs.labels}, chanlist);
            [tmp1 ordertmp2] = sort(tmp2);
            neworder = neworder(ordertmp2);
            opts = { opts{:} 'interp' allocs(neworder) };
        else
            newchanlist = [];
            tmpchanlocs = ALLEEG(idat(1)).chanlocs;
            chanlocs = { tmpchanlocs.labels };
            for i=1:length(chanlist)
                newchanlist = [ newchanlist strmatch(chanlist{i}, chanlocs, 'exact') ];
            end;
            tmpchanlocs =  ALLEEG(idat(1)).chanlocs;
            tmpchanlist = { tmpchanlocs(newchanlist).labels };
        end;

    % compute full file names
    % -----------------------
    function res = computeFullFileName(filePaths, fileNames);
        for index = 1:length(fileNames)
            res{index} = fullfile(filePaths{index}, fileNames{index});
        end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%ROTX Rotation about X axis
%
% R = ROTX(THETA) is an SO(3) rotation matrix (3x3) representing a rotation of THETA
% radians about the x-axis.
%
% R = ROTX(THETA, 'deg') as above but THETA is in degrees.
%
% See also ROTY, ROTZ, ANGVEC2R, ROT2.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function R = rotx(t, deg)

    if nargin > 1 && strcmp(deg, 'deg')
        t = t *pi/180;
    end

    ct = cos(t);
    st = sin(t);
    R = [
        1   0    0
        0   ct  -st
        0   st   ct
        ];

author: Allopart

repo: rbpf-gmapping


text:
function feq = imequi(f, iscale, hgram)
%------------------------------------------------------------------------------
% imequi
% Under the assumption that f is of class double, we linearly transform f into
% an intensity image according to the definition of the Image Processing Toolbox
% meaning that the values are contained in the range [0, 1]. Moreover, the
% contrast is enhanced. In case the Image Processing Toolbox is available this
% is performed by histogram equalization, if this toolbox is not available it is
% performed by removing the outliers.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: July 31, 2003.
%  2003 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
if nargin ~=2  && nargin ~=3
  error(' imequi - number of arguments should be either 2 or 3 ');
else
  if exist('histeq','file') ~= 2 && iscale == 3
%   Note: if histeq exists then imshow exists as well (in same toolbox).
    isca = 2;
    disp(' imequi - WARNING histeq does not exist ');
  else
    isca = iscale;
  end
end
%
if isca == 1
  feq = intensim(f);
elseif isca == 2
  showmean=mean(f(:));           showstd=std(f(:));
  showlow=showmean-2.0*showstd;  showhgh=showmean+2.0*showstd;
  supermin=min(min(f));          supermax=max(max(f));
  showlow=max(supermin,showlow); showhgh=min(showhgh,supermax);
  fs = (f < showlow)*showlow + (f > showhgh)*showhgh + ...
       (f >= showlow).*(f <= showhgh).*f;
  feq = intensim(fs);
elseif isca == 3
  if nargin == 3
    feq = histeq(intensim(f), hgram);
  else
    feq = histeq(intensim(f));
  end
else
  error(' imequi - value of argument iscale should be either 1 or 2 or 3 ');
end
%------------------------------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
addpath('bin');
% volsz = [512, 512, 255];
volsz = [100, 100, 100];
% volsz = [250, 250, 200];
spacing = [8, 8, 8];
spacing = [4, 4, 4];
k = ceil(volsz ./ spacing) + 1;
Knots = (rand([k, 3]) - 0.5) * 20;
sKnots = single(Knots);
dKnots = double(Knots);
gsKnots = gpuArray(sKnots);
gdKnots = gpuArray(dKnots);
vol = single(rand(volsz))*0;
vol(20:end-20, 20:end-20, 20:end-20) = 1;
dvol = double(vol);
gvol = gpuArray(vol);
gdvol = gpuArray(dvol);
%%
tic
disp = mex_3d_linear_disp_float(sKnots, single(volsz), single(spacing));
toc

tic
dispg = mex_3d_linear_disp_GPU_float(gsKnots, single(volsz), single(spacing));
toc

% dispg = gather(dispg);
sum(abs(disp(:) - gather(dispg(:))))
max(abs(disp(:) - gather(dispg(:))))
%%

tic
vdef = mex_3d_volume_warp_float(vol, disp, single(0));
toc

tic
vdefg = mex_3d_volume_warp_GPU_float(gvol, dispg, single(0));
toc
max(abs(vdef(:) - vdefg(:)))

%%
display('---');
tic
for i = 1 : 10
disp = mex_3d_linear_disp_float(sKnots, single(volsz), single(spacing));
vdef = mex_3d_volume_warp_float(vol, disp, single(0));
end
toc
%%
display('...');
tic
for i = 1 : 10
dispg = mex_3d_linear_disp_GPU_float(gsKnots, single(volsz), single(spacing));
vdefg = mex_3d_volume_warp_GPU_float(gvol, dispg, single(0));
end
toc

%%
volsz = [200, 255, 257];
spacing = [3, 3, 3];
ksz = ceil(volsz ./ spacing) + 1;

G1 = rand(volsz);
G2 = rand(volsz) * 0;
G2(30:end-30, 30:end-30, 30:end-30) = 1;
G3 = rand(volsz);

G1 = single(G1);
G2 = single(G2);
G3 = single(G3);

gG1 = gpuArray(G1);
gG2 = gpuArray(G2);
gG3 = gpuArray(G3);

tic
[gr1, gr2, gr3] = mex_3d_linear_partial_conv_float(G1, G2, G3, single(ksz), single(spacing));
toc

tic
[ggr1, ggr2, ggr3] = mex_3d_linear_partial_conv_GPU_float(gG1, gG2, gG3, single(ksz), single(spacing));
toc

[max(abs(gr1(:) - ggr1(:))), max(abs(gr2(:) - ggr2(:))), max(abs(gr3(:) - ggr3(:)))]

%%
addpath('..');
% addpath('mex_functions/bin');

tic
for i = 1 : 30
    voldef = linear_disp_and_warp_3d((sKnots), (vol), spacing, 0, 0, 0);
end
toc
%%
tic
for i = 1 :30
voldefg = linear_disp_and_warp_3d((gdKnots), (gdvol), spacing, 0, 0, 0);
end
toc
%%
display('----')
tic
D = mex_3d_linear_disp_GPU_float(gsKnots, single(size(gvol)), single(spacing));
toc
tic
voldef1 = mex_3d_volume_warp_GPU_float(gvol, D, single(0));
toc
% tic
% voldef2 = mex_3d_volume_warp_GPU_float(gvol, D, single(0));
% toc
% tic
% voldef3 = mex_3d_volume_warp_GPU_float(gvol, D, single(0));
% toc

%%
addpath('bin');
mex  -output bin/mex_3d_cubic_disp_double CXX='g++' CXXFLAGS='\$CXXFLAGS -O3 -funroll-loops -fopenmp -I../../mex_helpers -I../../deformation_tools_cpp' LDFLAGS="\$LDFLAGS -fopenmp" mex_3d_cubic_disp.cpp ../cubic_grid.cpp
vols = [270,270,270];
spc = [8,8,8];
k = ceil(vols ./ spc) + 3;
Knots = rand([k, 3]);
Knots = Knots * 0;
Knots(2, 2, 2,1) = 5;
Knots(5, 5, 2,1) = 10;
% Knots(2, 3, 1,1) = -5;
tic
D = mex_3d_cubic_disp_double(Knots, vols, spc);
toc
subplot(121)
imagesc(D(:,:,1,1)); colorbar
subplot(122)
imagesc(D(:,:,2,1)); colorbar
%%
spc = 4;
b = -ones(spc*4, 1);
for i = 1 : 2*spc + 1
    t = (i-1)/spc;
    if i <= spc
        k = t;
        f = 3*k^3 - 5 * k * k + 2;
    else
        k = t - 1;
        f = -k^3 + 2*k^2 - k;
    end
    b(2*spc+i-1) = f/2;
    if i <= 2*spc
        b(2*spc - i+1) = f/2;
    end
end
plot(b(1:end-1), 'r.-');
size(b)
%%
% mex  -output bin/mex_3d_cubic_partial_conv_double CXX='g++' CXXFLAGS='\$CXXFLAGS -O3 -funroll-loops -fopenmp -I../../mex_helpers -I../../deformation_tools_cpp' LDFLAGS="\$LDFLAGS -fopenmp" mex_3d_cubic_partial_conv.cpp ../cubic_grid.cpp
tr = rand(270, 270, 270);
spc = [2,2,2];
spc = [8,8,8];
ksz = ceil(size(tr)./spc) + 3;
tic
[g1, g2, g3] = mex_3d_cubic_partial_conv_double(tr, tr, tr, ksz, spc);
toc
ksz2 = ceil(size(tr)./spc) + 1;
tic
% [gg1, gg2, gg3] = mex_3d_linear_partial_conv_double(tr, tr, tr, ksz2, spc);
toc

%% upsmpl test
vol = zeros(100, 100, 100);
spc = [8,8,8];
k1 = ceil( (size(vol)-1) ./ spc) + 3;
k2 = ceil( (2*size(vol)-1) ./ spc) + 3;
Knots1 = rand([k1, 3]);
Knots1(4:round(end), 4:round(end/2), :, 1) = 13;
Knots2 = refine_cubic_grid_3d(Knots1, spc, spc, [1,1,1], size(vol)*2, size(vol));
%
D1 = mex_3d_cubic_disp_double(Knots1, size(vol), spc);
D2 = mex_3d_cubic_disp_double(Knots2, 2*size(vol), spc);
subplot(221)
imagesc(D1(:,:, 1, 1));
subplot(222)
imagesc(D2(:,:, 1, 1));

subplot(223)
imagesc(Knots1(:,:, round(end/2), 1));
subplot(224)
imagesc(Knots2(:,:, round(end/2), 1));
%%

subplot(121);
hold off;
ti1 = 0:(size(vol, 1)-1);
ti2 = (0:(2*size(vol, 1)-1))/2;
plot(ti1, squeeze(D1(round(end/3),:, 1,1)), 'bo-')
hold on;
plot(ti2, squeeze(D2(round(end/3),:, 1,1)), 'r.-')
subplot(121);
t1 = interp1(ti1, squeeze(D1(33,:, 1,1)), ti2);
plot(ti2, t1 - squeeze(D2(round(end/3),:,1,1)), 'k');
grid on;
subplot(122);
hold off
plot( (0:(k1(1)-1))-spc(1), squeeze(Knots1(15, :, 1,1)), 'bx-');
hold on;
plot((0:k2(1)-1)/2 - spc(1)*7/8, squeeze(Knots2(15, :, 1,1)), 'ro-');
%%
hold off;
plot(1:size(vol, 1), D1(:, 1, 1, 1), 'r.-');
hold on;
Ktmp =Knots1(2:end-1, 2, 2, 1);
plot(1 + (0:size(Ktmp,1)-1)*spc(1), Ktmp, 'bo-');

Ktmp =Knots2(2:end-1, 2, 2, 1);
plot(1 + (0:size(Ktmp,1)-1)*spc(1)/2, Ktmp, 'ks-');

repo: pTVreg


text:
import matlab.unittest.TestSuite
import matlab.unittest.TestRunner
import matlab.unittest.plugins.CodeCoveragePlugin
import matlab.unittest.plugins.TestReportPlugin

clc;
suite = TestSuite.fromFile('coverage/CMatrixTest.m');
runner = TestRunner.withTextOutput;
runner.addPlugin(CodeCoveragePlugin.forFile({'include/ddouble.m', 'include/AdaptiveChol.m'}))
result = runner.run(suite);

author: opencobra

repo: cobratoolbox


text:
function [weight] = ubah_kernel(data, weight_o)
% This file is the kernel core for the market strategy.
% If multiple experts are required, implement in this file.
% Here, only one expert is designed.
%
% function [weight] = ubah_kernel(data, weight_o)
% weight: generate portfolio for day t+1
% data: market sequence until day t
% weight_o: last (t^th) portfolio
%
% Example: [weight] ...
%          = ubah_kernel(data(1:t), weight_o);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

weight = ubah_expert(data, weight_o);

end
%%%%%%%%%%%%%%End%%%%%%%%%%%%%%%%%%%%%%

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function [f,str] = constructpolynomialmatrix2d(matrixsize,locs,degree)

% function [f,str] = constructpolynomialmatrix2d(matrixsize,locs,degree)
%
% <matrixsize> is a 2D matrix size like [100 50]
% <locs> is a row or column vector of indices into that matrix size
% <degree> is the maximum polynomial degree desired.
%   should be no greater than 10.
%
% return <f>, a matrix of dimensions length(<locs>) x N
% with polynomial basis functions evaluated at <locs> in
% the columns.  the polynomial basis functions are evaluated
% over the range [-1,1] which is presumed to correspond to
% the beginning and ending element along each of the two dimensions.
% (if a dimension has only one element, the values are all set to 1.)
% also, return <str>, the algebraic expression that corresponds to
% the columns of <f>.  'x' refers to the first matrix dimension; 'y'
% refers to the second matrix dimension.
%
% note that there may be various gain factors on the basis functions
% (e.g. don't assume that they are all unit-length).
%
% also, beware of numerical precision issues for high degrees...
%
% see also constructpolynomialmatrix3d.m.
%
% history:
% - 2013/05/23 - hard code the basis expansion to ensure consistent results across platforms.
%                this changes previous behavior.
%
% example:
% [f,str] = constructpolynomialmatrix2d([30 30],find(ones(30,30)),2);
% str
% f = reshape(f,30,30,[]);
% figure; imagesc(makeimagestack(f));

% prep
x = sym('x');
y = sym('y');

% do the algebra
%str = char(expand((x+y+1)^degree));
switch degree
case 0
  str = '1';
case 1
  str = 'x+y+1';
case 2
  str = 'x^2+2*x*y+2*x+y^2+2*y+1';
case 3
  str = 'x^3+3*x^2*y+3*x^2+3*x*y^2+6*x*y+3*x+y^3+3*y^2+3*y+1';
case 4
  str = '1+4*x+4*y+6*x^2+12*x*y+6*y^2+4*x^3*y+6*x^2*y^2+12*x^2*y+4*x*y^3+12*x*y^2+x^4+4*x^3+y^4+4*y^3';
case 5
  str = '1+5*x+5*y+10*x^2+20*x*y+10*y^2+20*x^3*y+30*x^2*y^2+30*x^2*y+20*x*y^3+30*x*y^2+5*x^4+10*x^3+5*y^4+10*y^3+y^5+5*x*y^4+10*x^2*y^3+10*x^3*y^2+5*x^4*y+x^5';
case 6
  str = '1+6*x+6*y+15*x^2+30*x*y+15*y^2+60*x^3*y+90*x^2*y^2+60*x^2*y+60*x*y^3+60*x*y^2+15*x^4+20*x^3+15*y^4+20*y^3+6*x*y^5+y^6+6*y^5+30*x*y^4+15*x^2*y^4+60*x^2*y^3+60*x^3*y^2+20*x^3*y^3+15*x^4*y^2+30*x^4*y+6*x^5+x^6+6*x^5*y';
case 7
  str = '1+7*x+7*y+21*x^2+42*x*y+21*y^2+140*x^3*y+210*x^2*y^2+105*x^2*y+140*x*y^3+105*x*y^2+35*x^4+35*x^3+35*y^4+35*y^3+7*x*y^6+42*x*y^5+7*y^6+21*y^5+y^7+105*x*y^4+105*x^2*y^4+210*x^2*y^3+21*x^2*y^5+210*x^3*y^2+35*x^3*y^4+140*x^3*y^3+105*x^4*y^2+105*x^4*y+21*x^5+7*x^6+x^7+35*x^4*y^3+42*x^5*y+21*x^5*y^2+7*x^6*y';
case 8
  str = '1+8*x+8*y+28*x^2+56*x*y+28*y^2+280*x^3*y+420*x^2*y^2+168*x^2*y+280*x*y^3+168*x*y^2+70*x^4+56*x^3+70*y^4+56*y^3+56*x*y^6+168*x*y^5+8*x*y^7+28*y^6+56*y^5+8*y^7+y^8+280*x*y^4+28*x^2*y^6+420*x^2*y^4+560*x^2*y^3+168*x^2*y^5+560*x^3*y^2+280*x^3*y^4+560*x^3*y^3+56*x^3*y^5+420*x^4*y^2+280*x^4*y+56*x^5+28*x^6+8*x^7+x^8+70*x^4*y^4+280*x^4*y^3+168*x^5*y+168*x^5*y^2+56*x^5*y^3+56*x^6*y+28*x^6*y^2+8*x^7*y';
case 9
  str = '1+9*x+9*y+36*x^2+72*x*y+36*y^2+504*x^3*y+756*x^2*y^2+252*x^2*y+504*x*y^3+252*x*y^2+126*x^4+84*x^3+126*y^4+84*y^3+252*x*y^6+504*x*y^5+72*x*y^7+84*y^6+126*y^5+36*y^7+9*y^8+630*x*y^4+252*x^2*y^6+1260*x^2*y^4+1260*x^2*y^3+756*x^2*y^5+1260*x^3*y^2+1260*x^3*y^4+1680*x^3*y^3+504*x^3*y^5+1260*x^4*y^2+630*x^4*y+126*x^5+84*x^6+36*x^7+9*x^8+630*x^4*y^4+1260*x^4*y^3+504*x^5*y+756*x^5*y^2+504*x^5*y^3+252*x^6*y+252*x^6*y^2+72*x^7*y+x^9+y^9+9*x^8*y+36*x^7*y^2+84*x^6*y^3+126*x^5*y^4+126*x^4*y^5+84*x^3*y^6+36*x^2*y^7+9*x*y^8';
case 10
  str = '1+10*x+10*y+45*x^2+90*x*y+45*y^2+840*x^3*y+1260*x^2*y^2+360*x^2*y+840*x*y^3+360*x*y^2+210*x^4+120*x^3+210*y^4+120*y^3+840*x*y^6+1260*x*y^5+360*x*y^7+210*y^6+252*y^5+120*y^7+45*y^8+1260*x*y^4+1260*x^2*y^6+3150*x^2*y^4+2520*x^2*y^3+2520*x^2*y^5+2520*x^3*y^2+4200*x^3*y^4+4200*x^3*y^3+2520*x^3*y^5+3150*x^4*y^2+1260*x^4*y+252*x^5+210*x^6+120*x^7+45*x^8+3150*x^4*y^4+4200*x^4*y^3+1260*x^5*y+2520*x^5*y^2+2520*x^5*y^3+840*x^6*y+1260*x^6*y^2+360*x^7*y+10*x^9+10*y^9+x^10+90*x^8*y+360*x^7*y^2+840*x^6*y^3+1260*x^5*y^4+1260*x^4*y^5+840*x^3*y^6+360*x^2*y^7+90*x*y^8+10*x^9*y+45*x^8*y^2+y^10+120*x^7*y^3+210*x^6*y^4+252*x^5*y^5+210*x^4*y^6+120*x^3*y^7+45*x^2*y^8+10*x*y^9';
otherwise
  die;
end

%REMOVE since hard-coded now
%str = sort(strsplit(str,'+'));%% sort the stuff in between + signs to try to ensure consistent ordering!!!
str = strsplit(str,'+');
str = cat(1,str,repmat({'+'},[1 length(str)]));
str = cat(2,str{:});
str = str(1:end-1);

% add a little padding so the 1 step below will work for degree 0
str = [' ' str ' '];

% change * to .*
  old = pwd;  % THIS IS A TOTAL HACK TO AVOID FUNCTION NAME CONFLICTS!
  cd(fullfile(matlabroot,'toolbox','matlab','funfun'));
str = vectorize(str);
  cd(old);

% remove +
str(str=='+') = ' ';

% change 1 to ones(size(x),1)
str0 = strrep(str,' 1 ',' ones(size(x)) '); assert(length(str0) ~= length(str));
str = str0;

% add brackets
str = [ '[' str ']' ];

% prep the linear coordinates
[x,y] = ind2sub(matrixsize,locs(:));
if matrixsize(1)~=1
  x = normalizerange(x,-1,1,1,matrixsize(1));
end
if matrixsize(2)~=1
  y = normalizerange(y,-1,1,1,matrixsize(2));
end

% do it
f = eval(str);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function g = findBestGammaAsCRF(img0, img1, exposure0, exposure1)
%
%       g = findBestGammaAsCRF(img0, img1, exposure0, exposure1)
%
%       This function computes the best gamma paramters for RGB that approximate
%       the camera response function.
%
%        Input:
%           -img0: an SDR image
%           -img1: an SDR image
%           -exposure0: the shutter speed of img0
%           -exposure1: the shutter speed of img1
%
%        Output:
%           -lin_fun: the inverse CRF
%           -max_lin_fun: maximum value of the inverse CRF
%
%     Copyright (C) 2014-15  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

delta_exposure = exposure1 / exposure0;
mask = ones(size(img0));
mask(img0>0.95) = 0.0;
mask(img1>0.95) = 0.0;
mask(img0<0.05) = 0.0;
mask(img1<0.05) = 0.0;

    function err = residualFunction(p)
        err = 0.0;
        for i=1:3
            gt = p(i);
            img0_lin(:,:,i) = (img0(:,:,i).^gt);
            img0_re(:,:,i) = img0_lin(:,:,i) * delta_exposure;
            img0_re(:,:,i) = img0_re(:,:,i).^(1.0/gt);

            delta = abs(img0_re(:,:,i) - img1(:,:,i)) .* mask(:,:,i);
            err = err + mean(delta(:));
        end
    end

    opts = optimset('Display', 'iter', 'TolFun', 1e-12, 'TolX', 1e-12, 'MaxIter', 1000);
    gi = 2.2 * ones(1,3);
    g = fminsearch(@residualFunction, gi, opts);
    for i=1:3
        img0_lin(:,:,i) = (img0(:,:,i).^g(i));
        img0_re(:,:,i) = img0_lin(:,:,i) * delta_exposure;
        img0_re(:,:,i) = img0_re(:,:,i).^(1.0/g(i));
    end

    imshow([img0, img0_re, img1]);

end

author: banterle


text:
function data_rx = extract_data(symbol_ofdm_rx,v_pilots)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                       %
%%       Name: extract_data.m                                            %
%%                                                                       %
%%       Description: In this function, the inverse process of the       %
%%        function "createsymbol.m" is realized. That is, we extract     %
%%        the known parts of the symbol that we have received:           %
%%        --> Pilot carriers                                             %
%%        --> Data carriers                                              %
%%        --> Guard carriers                                             %
%%                                                                       %
%%       Result: The data is given back only, since the rest does not    %
%%        interest us.                                                   %
%%                                                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% First, the position of the data is located
 v_data = setxor(1:length(symbol_ofdm_rx),v_pilots);

% Next, the values of the data and pilot carriers are extracted
 data_total = symbol_ofdm_rx (v_data);

% Now the dc carrier must be removed
 data_rx = [data_total(29:124) data_total(126:221)];


text:
function [label,depth] = RenderPCcameraMatlab(XYZ,RGB,XYZedges,RGBedges,K, Rt,width, height)

% XYZ should be 3*N double
% RGB should be 3*N uint8

%mex WarpMesh.cpp -lGLU -lOSMesa
%load debug.mat

if ~isempty(XYZ)
    XYZ = double(XYZ);
    XYZ = transformPointCloud(XYZ,Rt);
    XYZ(1,:) = -XYZ(1,:);
end

if ~isempty(XYZedges)
    XYZedges = double(XYZedges);
    XYZedges = transformPointCloud(XYZedges,Rt);
    XYZedges(1,:) = -XYZedges(1,:);
end

%{
figure(5); clf;
plot3(X(1,:),X(2,:),X(3,:),'.');
axis equal;
grid on;
xlabel('x');
ylabel('y');
zlabel('z');
set(gca,'ZDir','reverse');
set(gca,'YDir','reverse');
%}
P = [K [0;0;0]];
[label,depth]=RenderPCcamera(P,width,height,XYZ,RGB,XYZedges,RGBedges);

%{

% to deal with camera coodinate system
%XYZcamera(:,:,1) = -XYZcamera(:,:,1);
%P = K * Rt;
%P = P * [-1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];

%P = Rt;

XYZcamera(:,:,2) = -XYZcamera(:,:,2);
XYZcamera(:,:,3) = -XYZcamera(:,:,3);

F = [1 0 0 ; 0 -1 0; 0 0 -1];

P = [ F*Rt(:,1:3)*inv(F) F*Rt(:,4)];
%P = Rt; % * [1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 1];

%{
XYZcamera(:,:,2) = -XYZcamera(:,:,2);
XYZcamera(:,:,3) = -XYZcamera(:,:,3);
P = [1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 1] * [Rt; 0 0 0 1] * [1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 1];
P = P(1:3,:);
%}

good = XYZcamera(:,:,4)~=0;
good = good(:);

X = reshape(XYZcamera,[],4)';
X = X(1:3,good);

figure(5); clf;
plot3(X(1,:),X(2,:),X(3,:),'.');
axis equal;
grid on;
xlabel('x');
ylabel('y');
zlabel('z');
set(gca,'ZDir','reverse');
set(gca,'YDir','reverse');

X = transformPointCloud(X,P);

figure(3);clf;
plot3(X(1,:),X(2,:),X(3,:),'.');
axis equal;
grid on;
xlabel('x');
ylabel('y');
zlabel('z');
set(gca,'ZDir','reverse');
set(gca,'YDir','reverse');

figure(4);clf;
x= K*X;
x = [x(1,:)./x(3,:); x(2,:)./x(3,:)];
plot(x(1,:),x(2,:),'.')
axis equal
set(gca,'YDir','reverse');
hold on;
plot([1 640 640 1 1],[1 1 480 480 1],'-r');

P = K * P;

%{
XYZcamera(:,:,2) = -XYZcamera(:,:,2);
XYZcamera(:,:,3) = -XYZcamera(:,:,3);
%P = P * [1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 1];
Rt(2,4) = -Rt(2,4);
Rt(3,4) = -Rt(3,4);
P = K * Rt;
%}

[label,depth]=WarpMesh(P,640,480,XYZcamera);

%}

label = label';
label = label(:,end:-1:1);
%figure
%imagesc(label)

label = double(label);

B = mod(label,256);
G = mod(floor(label/256),256);
R = floor(label/(256*256));

clear label;
label(:,:,1) = R;
label(:,:,2) = G;
label(:,:,3) = B;
label = uint8(label);

depth = depth';
depth = depth(end:-1:1,end:-1:1);

if ~any(label(:))
    depth = single(zeros(size(depth)));
    return
end

z_near = 0.3;
z_far_ratio = 1.2;
depth = z_near./(1-single(depth)/2^32);
maxDepth = max(depth(abs(depth) < 100));

cropmask = (depth < z_near) | (depth > z_far_ratio * maxDepth);
depth(cropmask) = 0; %NaN;%z_far_ratio * maxDepth;

%figure
%imagesc(depth)

author: jianxiongxiao

repo: ProfXkit


text:
%if exist('images_read');
%   active_images = active_images & images_read;
%end;

var2fix = 'dX_default';

fixvariable;

var2fix = 'dY_default';

fixvariable;

var2fix = 'map';

fixvariable;

if ~exist('n_ima'),
    data_calib_no_read;
end;

check_active_images;

% Step used to clean the memory if a previous atttempt has been made to read the entire set of images into memory:
for kk = 1:n_ima,
    if (exist(['I_' num2str(kk)])==1),
        clear(['I_' num2str(kk)]);
    end;
end;

fprintf(1,'\nExtraction of the grid corners on the images\n');

if (exist('map')~=1), map = gray(256); end;

%disp('WARNING!!! Do not forget to change dX_default and dY_default in click_calib.m!!!')

if exist('dX'),
    dX_default = dX;
end;

if exist('dY'),
    dY_default = dY;
end;

if exist('n_sq_x'),
    n_sq_x_default = n_sq_x;
end;

if exist('n_sq_y'),
    n_sq_y_default = n_sq_y;
end;

if ~exist('dX_default')|~exist('dY_default');

    % Setup of JY - 3D calibration rig at Intel (new at Intel) - use units in mm to match Zhang
    dX_default = 30;
    dY_default = 30;

    % Setup of JY - 3D calibration of the Rosette (Google) - use units in mm to match Zhang
    dX_default = 0.1524; %100;
    dY_default = 0.1524; %100;

end;

if ~exist('n_sq_x_default')|~exist('n_sq_y_default'),
    n_sq_x_default = 10;
    n_sq_y_default = 10;
end;

if ~exist('wintx_default')|~exist('winty_default'),
    if ~exist('nx'),
        wintx_default = 20;
        winty_default = wintx_default;
        clear wintx winty
    else
        wintx_default = max(round(nx/128),round(ny/96));
        winty_default = wintx_default;
        clear wintx winty
    end;
end;

if ~exist('wintx') | ~exist('winty'),
    clear_windows; % Clear all the window sizes (to re-initiate)
end;

if ~exist('dont_ask'),
    dont_ask = 0;
end;

if ~dont_ask,
    ima_numbers = input('Number(s) of image(s) to process ([] = all images) = ');
else
    ima_numbers = [];
end;

if isempty(ima_numbers),
    ima_proc = 1:n_ima;
else
    ima_proc = ima_numbers;
end;

% Useful option to add images:
kk_first = ima_proc(1); %input('Start image number ([]=1=first): ');

%if isempty(kk_first), kk_first = 1; end;

if exist(['wintx_' num2str(kk_first)]),

    eval(['wintxkk = wintx_' num2str(kk_first) ';']);

    if isempty(wintxkk) | isnan(wintxkk),

        disp('Window size for corner finder (wintx and winty):');
        wintx = input(['wintx ([] = ' num2str(wintx_default) ') = ']);
        if isempty(wintx), wintx = wintx_default; end;
        wintx = round(wintx);
        winty = input(['winty ([] = ' num2str(winty_default) ') = ']);
        if isempty(winty), winty = winty_default; end;
        winty = round(winty);

        fprintf(1,'Window size = %dx%d\n',2*wintx+1,2*winty+1);

    end;

else

    disp('Window size for corner finder (wintx and winty):');
    wintx = input(['wintx ([] = ' num2str(wintx_default) ') = ']);
    if isempty(wintx), wintx = wintx_default; end;
    wintx = round(wintx);
    winty = input(['winty ([] = ' num2str(winty_default) ') = ']);
    if isempty(winty), winty = winty_default; end;
    winty = round(winty);

    fprintf(1,'Window size = %dx%d\n',2*wintx+1,2*winty+1);

end;

if ~dont_ask,
    fprintf(1,'Do you want to use the automatic square counting mechanism (0=[]=default)\n');
    manual_squares = input('  or do you always want to enter the number of squares manually (1,other)? ');
    if isempty(manual_squares),
        manual_squares = 0;
    else
        manual_squares = ~~manual_squares;
    end;
else
    manual_squares = 0;
end;

for kk = ima_proc,

    if ~type_numbering,
        number_ext =  num2str(image_numbers(kk));
    else
        number_ext = sprintf(['%.' num2str(N_slots) 'd'],image_numbers(kk));
    end;

    ima_name = [calib_name  number_ext '.' format_image];

    if exist(ima_name),

        fprintf(1,'\nProcessing image %d...\n',kk);

        fprintf(1,'Loading image %s...\n',ima_name);

        if format_image(1) == 'p',
            if format_image(2) == 'p',
                I = double(loadppm(ima_name));
            else
                I = double(loadpgm(ima_name));
            end;
        else
            if format_image(1) == 'r',
                I = readras(ima_name);
            else
                I = double(imread(ima_name));
            end;
        end;

        if size(I,3)>1,
            I = 0.299 * I(:,:,1) + 0.5870 * I(:,:,2) + 0.114 * I(:,:,3);
        end;

        [ny,nx,junk] = size(I);
        Wcal = nx; % to avoid errors later
        Hcal = ny; % to avoid errors later

        click_ima_calib_fisheye_no_read;

        active_images(kk) = 1;

    else
        eval(['dX_' num2str(kk) ' = NaN;']);
        eval(['dY_' num2str(kk) ' = NaN;']);

        eval(['wintx_' num2str(kk) ' = NaN;']);
        eval(['winty_' num2str(kk) ' = NaN;']);

        eval(['x_' num2str(kk) ' = NaN*ones(2,1);']);
        eval(['X_' num2str(kk) ' = NaN*ones(3,1);']);

        eval(['n_sq_x_' num2str(kk) ' = NaN;']);
        eval(['n_sq_y_' num2str(kk) ' = NaN;']);
    end;
end;

check_active_images;

% Fix potential non-existing variables:

for kk = 1:n_ima,
    if ~exist(['x_' num2str(kk)]),
        eval(['dX_' num2str(kk) ' = NaN;']);
        eval(['dY_' num2str(kk) ' = NaN;']);

        eval(['x_' num2str(kk) ' = NaN*ones(2,1);']);
        eval(['X_' num2str(kk) ' = NaN*ones(3,1);']);

        eval(['n_sq_x_' num2str(kk) ' = NaN;']);
        eval(['n_sq_y_' num2str(kk) ' = NaN;']);
    end;

    if ~exist(['wintx_' num2str(kk)]) | ~exist(['winty_' num2str(kk)]),

        eval(['wintx_' num2str(kk) ' = NaN;']);
        eval(['winty_' num2str(kk) ' = NaN;']);

    end;
end;

string_save = 'save calib_data active_images ind_active wintx winty n_ima type_numbering N_slots first_num image_numbers format_image calib_name Hcal Wcal nx ny map dX_default dY_default dX dY';

for kk = 1:n_ima,
    string_save = [string_save ' X_' num2str(kk) ' x_' num2str(kk) ' n_sq_x_' num2str(kk) ' n_sq_y_' num2str(kk) ' wintx_' num2str(kk) ' winty_' num2str(kk) ' dX_' num2str(kk) ' dY_' num2str(kk)];
end;

eval(string_save);

disp('done');

repo: ekfmonoslam


text:
function [xo,N]=gabglasso(ttype,xi,lambda,group);
%GABGLASSO   group lasso estimate (hard/soft) in time-frequency domain
%   Usage:  xo=gabglasso(ttype,x,lambda,group);
%           [xo,N]=gabglasso(ttype,x,lambda,group));
%
%   GABGLASSO('hard',x,lambda,'time') will perform
%   time hard group thresholding on x, i.e. all time-frequency
%   columns whose norm less than lambda will be set to zero.
%
%   GABGLASSO('soft',x,lambda,'time') will perform
%   time soft thresholding on x, i.e. all time-frequency
%   columns whose norm less than lambda will be set to zero,
%   and those whose norm exceeds lambda will be multiplied
%   by (1-lambda/norm).
%
%   GABGLASSO(ttype,x,lambda,'frequency') will perform
%   frequency thresholding on x, i.e. all time-frequency
%   rows whose norm less than lambda will be soft or hard thresholded
%   (see above).
%
%   [xo,N]=GABGLASSO(ttype,x,lambda,group) additionally returns
%   a number N specifying how many numbers where kept.
%
%   The function may meaningfully be applied to output from DGT, WMDCT or
%   from WIL2RECT(DWILT(...)).
%
%   See also:  gablasso, gabelasso
%
%   Demos: demo_audioshrink

%   AUTHOR : Bruno Torresani.
%   REFERENCE: OK

complainif_argnonotinrange(nargin,4,4,mfilename);

NbFreqBands = size(xi,1);
NbTimeSteps = size(xi,2);

xo = zeros(size(xi));

switch(lower(group))
 case {'time'}
  for t=1:NbTimeSteps,
    threshold = norm(xi(:,t));
    mask = (1-lambda/threshold);
    if(strcmp(ttype,'soft'))
      mask = mask * (mask>0);
    elseif(strcmp(ttype,'hard'))
      mask = (mask>0);
    end
    xo(:,t) = xi(:,t) * mask;
  end
 case {'frequency'}
  for f=1:NbFreqBands,
    threshold = norm(xi(f,:));
    mask = (1-lambda/threshold);
    mask = mask * (mask>0);
    if(strcmp(ttype,'soft'))
      mask = mask * (mask>0);
    elseif(strcmp(ttype,'hard'))
      mask = (mask>0);
    end
    xo(f,:) = xi(f,:) * mask;
  end
 otherwise
  error('"group" parameter must be either "time" or "frequency".');
end

if nargout==2
    signif_map = (abs(xo)>0);
    N = sum(signif_map(:));
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function q = vertcat(varargin)
% implements [q1;q2;q3..]

q = cat(1,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function vw = meshUpdateAll(vw)
% Update all mesh displays associated with a mrVista view.
%
% vw = meshUpdateAll([vw=cur gray view]);
%
% ras, 08/2008.
if notDefined('vw'), vw = getSelectedGray; end

if ~isfield(vw, 'mesh') || isempty(viewGet(vw, 'mesh'))
	return
end

% remember the currently selected mesh, to restore it later
curMesh = viewGet(vw, 'mesh num');

for ii = 1:length(viewGet(vw, 'allmeshes'))
	vw = viewSet(vw, 'mesh num', ii);
	vw = meshColorOverlay(vw);
end

% restore the originally-selected mesh
vw = viewSet(vw, 'mesh num', curMesh);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [samples, energies, diagn] = hmc(f, x, options, gradf, varargin)
%HMC	Hybrid Monte Carlo sampling.
%
%	Description
%	SAMPLES = HMC(F, X, OPTIONS, GRADF) uses a  hybrid Monte Carlo
%	algorithm to sample from the distribution P ~ EXP(-F), where F is the
%	first argument to HMC. The Markov chain starts at the point X, and
%	the function GRADF is the gradient of the `energy' function F.
%
%	HMC(F, X, OPTIONS, GRADF, P1, P2, ...) allows additional arguments to
%	be passed to F() and GRADF().
%
%	[SAMPLES, ENERGIES, DIAGN] = HMC(F, X, OPTIONS, GRADF) also returns a
%	log of the energy values (i.e. negative log probabilities) for the
%	samples in ENERGIES and DIAGN, a structure containing diagnostic
%	information (position, momentum and acceptance threshold) for each
%	step of the chain in DIAGN.POS, DIAGN.MOM and DIAGN.ACC respectively.
%	All candidate states (including rejected ones) are stored in
%	DIAGN.POS.
%
%	[SAMPLES, ENERGIES, DIAGN] = HMC(F, X, OPTIONS, GRADF) also returns
%	the ENERGIES (i.e. negative log probabilities) corresponding to the
%	samples.  The DIAGN structure contains three fields:
%
%	POS the position vectors of the dynamic process.
%
%	MOM the momentum vectors of the dynamic process.
%
%	ACC the acceptance thresholds.
%
%	S = HMC('STATE') returns a state structure that contains the state of
%	the two random number generators RAND and RANDN and the momentum of
%	the dynamic process.  These are contained in fields  randstate,
%	randnstate and mom respectively.  The momentum state is only used for
%	a persistent momentum update.
%
%	HMC('STATE', S) resets the state to S.  If S is an integer, then it
%	is passed to RAND and RANDN and the momentum variable is randomised.
%	If S is a structure returned by HMC('STATE') then it resets the
%	generator to exactly the same state.
%
%	The optional parameters in the OPTIONS vector have the following
%	interpretations.
%
%	OPTIONS(1) is set to 1 to display the energy values and rejection
%	threshold at each step of the Markov chain. If the value is 2, then
%	the position vectors at each step are also displayed.
%
%	OPTIONS(5) is set to 1 if momentum persistence is used; default 0,
%	for complete replacement of momentum variables.
%
%	OPTIONS(7) defines the trajectory length (i.e. the number of leap-
%	frog steps at each iteration).  Minimum value 1.
%
%	OPTIONS(9) is set to 1 to check the user defined gradient function.
%
%	OPTIONS(14) is the number of samples retained from the Markov chain;
%	default 100.
%
%	OPTIONS(15) is the number of samples omitted from the start of the
%	chain; default 0.
%
%	OPTIONS(17) defines the momentum used when a persistent update of
%	(leap-frog) momentum is used.  This is bounded to the interval [0,
%	1).
%
%	OPTIONS(18) is the step size used in leap-frogs; default 1/trajectory
%	length.
%
%	See also
%	METROP
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Global variable to store state of momentum variables: set by set_state
% Used to initialise variable if set
global HMC_MOM
if nargin <= 2
  if ~strcmp(f, 'state')
    error('Unknown argument to hmc');
  end
  switch nargin
    case 1
      samples = get_state(f);
      return;
    case 2
      set_state(f, x);
      return;
  end
end

display = options(1);
if (round(options(5) == 1))
  persistence = 1;
  % Set alpha to lie in [0, 1)
  alpha = max(0, options(17));
  alpha = min(1, alpha);
  salpha = sqrt(1-alpha*alpha);
else
  persistence = 0;
end
L = max(1, options(7)); % At least one step in leap-frogging
if options(14) > 0
  nsamples = options(14);
else
  nsamples = 100;	% Default
end
if options(15) >= 0
  nomit = options(15);
else
  nomit = 0;
end
if options(18) > 0
  step_size = options(18);	% Step size.
else
  step_size = 1/L;		% Default
end
x = x(:)';		% Force x to be a row vector
nparams = length(x);

% Set up strings for evaluating potential function and its gradient.
f = fcnchk(f, length(varargin));
gradf = fcnchk(gradf, length(varargin));

% Check the gradient evaluation.
if (options(9))
  % Check gradients
  feval('gradchek', x, f, gradf, varargin{:});
end

samples = zeros(nsamples, nparams);	% Matrix of returned samples.
if nargout >= 2
  en_save = 1;
  energies = zeros(nsamples, 1);
else
  en_save = 0;
end
if nargout >= 3
  diagnostics = 1;
  diagn_pos = zeros(nsamples, nparams);
  diagn_mom = zeros(nsamples, nparams);
  diagn_acc = zeros(nsamples, 1);
else
  diagnostics = 0;
end

n = - nomit + 1;
Eold = feval(f, x, varargin{:});	% Evaluate starting energy.
nreject = 0;
if (~persistence | isempty(HMC_MOM))
  p = randn(1, nparams);		% Initialise momenta at random
else
  p = HMC_MOM;				% Initialise momenta from stored state
end
lambda = 1;

% Main loop.
while n <= nsamples

  xold = x;		    % Store starting position.
  pold = p;		    % Store starting momenta
  Hold = Eold + 0.5*(p*p'); % Recalculate Hamiltonian as momenta have changed

  if ~persistence
    % Choose a direction at random
    if (rand < 0.5)
      lambda = -1;
    else
      lambda = 1;
    end
  end
  % Perturb step length.
  epsilon = lambda*step_size*(1.0 + 0.1*randn(1));

  % First half-step of leapfrog.
  p = p - 0.5*epsilon*feval(gradf, x, varargin{:});
  x = x + epsilon*p;

  % Full leapfrog steps.
  for m = 1 : L - 1
    p = p - epsilon*feval(gradf, x, varargin{:});
    x = x + epsilon*p;
  end

  % Final half-step of leapfrog.
  p = p - 0.5*epsilon*feval(gradf, x, varargin{:});

  % Now apply Metropolis algorithm.
  Enew = feval(f, x, varargin{:});	% Evaluate new energy.
  p = -p;				% Negate momentum
  Hnew = Enew + 0.5*p*p';		% Evaluate new Hamiltonian.
  a = exp(Hold - Hnew);			% Acceptance threshold.
  if (diagnostics & n > 0)
    diagn_pos(n,:) = x;
    diagn_mom(n,:) = p;
    diagn_acc(n,:) = a;
  end
  if (display > 1)
    fprintf(1, 'New position is\n');
    disp(x);
  end

  if a > rand(1)			% Accept the new state.
    Eold = Enew;			% Update energy
    if (display > 0)
      fprintf(1, 'Finished step %4d  Threshold: %g\n', n, a);
    end
  else					% Reject the new state.
    if n > 0
      nreject = nreject + 1;
    end
    x = xold;				% Reset position
    p = pold;   			% Reset momenta
    if (display > 0)
      fprintf(1, '  Sample rejected %4d.  Threshold: %g\n', n, a);
    end
  end
  if n > 0
    samples(n,:) = x;			% Store sample.
    if en_save
      energies(n) = Eold;		% Store energy.
    end
  end

  % Set momenta for next iteration
  if persistence
    p = -p;
    % Adjust momenta by a small random amount.
    p = alpha.*p + salpha.*randn(1, nparams);
  else
    p = randn(1, nparams);	% Replace all momenta.
  end

  n = n + 1;
end

if (display > 0)
  fprintf(1, '\nFraction of samples rejected:  %g\n', ...
    nreject/(nsamples));
end
if diagnostics
  diagn.pos = diagn_pos;
  diagn.mom = diagn_mom;
  diagn.acc = diagn_acc;
end
% Store final momentum value in global so that it can be retrieved later
HMC_MOM = p;
return

% Return complete state of sampler (including momentum)
function state = get_state(f)

global HMC_MOM
state.randstate = rand('state');
state.randnstate = randn('state');
state.mom = HMC_MOM;
return

% Set complete state of sampler (including momentum) or just set randn
% and rand with integer argument.
function set_state(f, x)

global HMC_MOM
if isnumeric(x)
  rand('state', x);
  randn('state', x);
  HMC_MOM = [];
else
  if ~isstruct(x)
    error('Second argument to hmc must be number or state structure');
  end
  if (~isfield(x, 'randstate') | ~isfield(x, 'randnstate') ...
      | ~isfield(x, 'mom'))
    error('Second argument to hmc must contain correct fields')
  end
  rand('state', x.randstate);
  randn('state', x.randnstate);
  HMC_MOM = x.mom;
end
return

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function triangle_svg_test01 ( )

%*****************************************************************************80
%
%% TRIANGLE_SVG_TEST01 calls TRIANGLE_SVG to plot a triangle and some points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2014
%
%  Author:
%
%    John Burkardt
%
  p_num = 13;
  p = [ ...
       0.333333333333333, 0.333333333333333; ...
       0.479308067841923, 0.260345966079038; ...
       0.260345966079038, 0.479308067841923; ...
       0.260345966079038, 0.260345966079038; ...
       0.869739794195568, 0.065130102902216; ...
       0.065130102902216, 0.869739794195568; ...
       0.065130102902216, 0.065130102902216; ...
       0.638444188569809, 0.312865496004875; ...
       0.638444188569809, 0.048690315425316; ...
       0.312865496004875, 0.638444188569809; ...
       0.312865496004875, 0.048690315425316; ...
       0.048690315425316, 0.638444188569809; ...
       0.048690315425316, 0.312865496004875 ]';
  t = [ ...
       0.0, 0.0; ...
       1.0, 0.0; ...
       0.0, 1.0 ]';

  plot_filename = 'test01.svg';

  triangle_svg ( plot_filename, t, p_num, p );

  return
end

author: johannesgerer

repo: jburkardt-m