text:
classdef PerimeterInteriorCircleExperiment < handle

    methods (Access = public)

        function obj = PerimeterInteriorCircleExperiment()
            obj.computeGradientVariationWithRadius();
            obj.computeGradientVariationWithBoundary();
            obj.computeGradientSurfPerimeterComputer();
        end

    end

    methods (Access = private)

        function computeGradientVariationWithBoundary(obj)
            s = obj.createParams();
            g = GradientVariationWithBoundaryExperiment(s);
            g.compute();
        end

        function computeGradientVariationWithRadius(obj)
            s = obj.createParams();
            g = GradientVariationWithRadiusExperiment(s);
            g.compute();
        end

        function computeGradientSurfPerimeterComputer(obj)
            s = obj.createParams();
            g = GradientSurfPerimeterComputer(s);
            g.compute();
        end

    end

    methods (Access = private, Static)

        function s = createParams()
            radius = 0.25;
            s.circleCase = 'interior';
            s.levelSetParams.fracRadius = radius;
            s.levelSetParams.type = 'circle';
            L = 1;
            halfSide = L/2;
            s.levelSetParams.fracRadius = radius/halfSide;
            s.curvature = 1/radius;
            s.circleCase = 'interior';
            s.nameCase = 'GradientCirclePerimeterExperiment';
            nameRoot   = 'SquareMacroTriangle';
            s.inputFiles = {nameRoot; [nameRoot,'Fine'];[nameRoot,'FineFine']};
            s.outputFolder = '/home/alex/git-repos/Perimeter/AllImages/';
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%Help file for INTLAB Version 5.1
%
% - corrections in long toolbox (thanks to Nobito Yamamoto and Nozomu Matsuda)
% - new function 'isspd': verify positive definiteness
% - revision of verifylss for sparse matrices (omit eigs etc.)
% - warning when defining complex by infsup(zinf,zsup) [overestimation, use midrad]
% - abss replaced by mag (thanks to Arnold for proposing better naming)
% - generation of extremely ill-conditioned matrices (randmat) improved
% - norm delivers interval result
% - rounding unchanged after call of any INTLAB routine
%

repo: OpenMAS


text:
function value = p15_i4 ( action, name, value )

%*****************************************************************************80
%
%% P15_I4 sets or gets I4 parameters for problem 15.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'D' to set a parameter to its default value.
%    'S' to set a parameter,
%    'G' to get a parameter,
%    'I' to increment a parameter.
%
%    Input, string NAME, the name of the variable.
%    '#' is the number of calls to the integrand routine.
%    'N' is the number of factors to multiply in the integrand.
%
%    Input/output, integer VALUE.
%    * If ACTION = 'I', then VALUE is an input quantity, and is the
%      new value to be added to NAME.
%    * If ACTION = 'S', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'G', then VALUE is an output quantity, and is the
%      current value of NAME.
%
  persistent calls;
  persistent n;

  if ( action == 'D' | action == 'd' )

    if ( name == '#' | name == '*' )
      calls = 0;
    end

    if ( name == 'n' | name == 'N' | name == '*' )
      n = 1;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == '#' )
      value = calls;
    elseif ( name == 'n' | name == 'N' )
      value = n;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P15_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P15_I4 - Fatal error!' );
    end

  elseif ( action == 'I' | action == 'i' )

    if ( name == '#' )
      calls = calls + value;
    elseif ( name == 'n' | name == 'N' )
      n = n + value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P15_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P15_I4 - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == '#' )
      calls = value;
    elseif ( name == 'n' | name == 'N' )
      n = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P15_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P15_I4 - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P15_I4 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P15_I4 - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [] = ResumeExperiment( jobID, taskID, Niter )
%INPUT
%  jobID : integer/name of job
%  taskID : integer/name of task
%  Niter  : # iterations of MCMC

resumeBPHMM( {jobID, taskID}, {'Niter', Niter}, {'jobID', jobID, 'taskID', taskID} );

author: michaelchughes

repo: NPBayesHMM


text:
function term_mat = zpk2term(z,p,k,T)
%
% Utility Function: ZPK2TERM
%
% The purpose of this function is to convert zero/pole/gain format into
% the term matrix format used by the environments

% Author: Craig Borghesani
% Date: 9/5/94
% Copyright (c) 1999, Prentice-Hall

narg_vals = nargin;
if narg_vals==4,
 if ~length(T),
  narg_vals=3;
 end
end

ctz=length(z);
ctp=length(p);
ct=1;
term_mat(1,:) = [k,0,NaN,1];
term_mat(2,:) = [0,NaN,NaN,2];

if narg_vals==3,
 while ct<=ctz,
  if z(ct)==0,
   term_mat(2,1)=term_mat(2,1)-1;
   ct=ct+1;
  elseif imag(z(ct))==0,
   term_mat=[term_mat;-z(ct),NaN,NaN,5];
   ct=ct+1;
  else
   re=real(z(ct)); im=imag(z(ct));
   wn=sqrt(re^2+im^2); zta=-re/wn;
   term_mat=[term_mat;zta wn NaN 7];
   ct=ct+2;
  end
 end
 ct=1;
 while ct<=ctp,
  if p(ct)==0,
   term_mat(2,1)=term_mat(2,1)+1;
   ct=ct+1;
  elseif imag(p(ct))==0,
   term_mat=[term_mat;-p(ct),NaN,NaN,4];
   ct=ct+1;
  else
   re=real(p(ct)); im=imag(p(ct));
   wn=sqrt(re^2+im^2); zta=-re/wn;
   term_mat=[term_mat;zta wn NaN 6];
   ct=ct+2;
  end
 end
else
 while ct<=ctz,
  if z(ct)==0,
   term_mat(3,1)=term_mat(3,1)-1;
   ct=ct+1;
  elseif imag(z(ct))==0,
   if z(ct)==1,
    term_mat(2,2)=term_mat(2,2)+1;
   else
    term_mat=[term_mat;exp(-z(ct)*T),NaN,NaN,2];
   end
   ct=ct+1;
  else
   z(ct)=log(z(ct))/T;
   re=real(z(ct)); im=imag(z(ct));
   wn=sqrt(re^2+im^2); zta=-re/wn;
   term_mat=[term_mat;zta wn NaN 4];
   ct=ct+2;
  end
 end
 ct=1;
 while ct<=ctp,
  if p(ct)==0,
   term_mat(3,1)=term_mat(3,1)+1;
   ct=ct+1;
  elseif imag(p(ct))==0,
   if p(ct)==1,
    term_mat(2,1)=term_mat(2,1)+1;
   else
    term_mat=[term_mat;exp(-p(ct)*T),NaN,NaN,1];
   end
   ct=ct+1;
  else
   p(ct)=log(p(ct))/T;
   re=real(p(ct)); im=imag(p(ct));
   wn=sqrt(re^2+im^2); zta=-re/wn;
   term_mat=[term_mat;zta wn NaN 3];
   ct=ct+2;
  end
 end
end

author: Sable

repo: mcbench-benchmarks


text:
function [xUpdate,SUpdate,innov,Szz,W]=sqrtCubKalUpdateWithPred(z,SR,zPred,otherInfo)
%%SQRTCUBKALUPDATEWITHPRED Given the output of the measurement prediction
%           step from sqrtCubKalMeasPred and a measurement, complete the
%           measurement update step of the square root cubature Kalman
%           filter. Separating the measurement prediction step from the
%           rest of the update step can make the creation of multiple
%           measurement association hypotheses from a single target
%           prediction more efficient. The full measurement update function
%           is sqrtCubKalUpdate.
%
%INPUTS: z The zDimX1 measurement vector.
%       SR The zDimXzDim lower-triangular square root of the measurement
%          covariance matrix in the native coordinate system of the
%          measurement.
%    zPred The zDimXnumComp measurement predictions from the filter.
% otherInfo The intermediate results returned in the otherInfo output of
%          the sqrtCubKalMeasPred function.
%
%OUTPUTS: xUpdate The xDimXnumComp updated state vectors.
%         SUpdate The updated xDimXxDimXnumComp lower-triangular square-
%                 root state covariance matrices.
%      innov, Szz The zDimXnumComp innovations and the zDimXzDimXnumComp
%                 square-root innovation covariance matrices are returned
%                 in case one wishes to analyze the consistency of the
%                 estimator or use those values in gating or likelihood
%                 evaluation.
%               W The xDimXzDimXnumComp gains used in the update.
%
%See the comments to the function sqrtCubKalMeasPred for an example of
%usage of this function. See the comments to sqrtCubKalUpdate for more
%information on the algorithm.
%
%June 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    xPred=otherInfo.xPred;
    innovTrans=otherInfo.innovTrans;
    stateDiffTrans=otherInfo.stateDiffTrans;
    xPredCenPoints=otherInfo.xPredCenPoints;
    stateTrans=otherInfo.stateTrans;
    zPredCenPoints=otherInfo.zPredCenPoints;
    Pxz=otherInfo.Pxz;

    xDim=size(xPred,1);
    numComp=size(xPred,2);
    zDim=size(z,1);

    xUpdate=zeros(xDim,numComp);
    SUpdate=zeros(xDim,xDim,numComp);
    innov=zeros(zDim,numComp);
    Szz=zeros(zDim,zDim,numComp);
    W=zeros(xDim,zDim,numComp);

    for k=1:numComp
        Szz(:,:,k)=tria([zPredCenPoints(:,:,k),SR]);

        %The filter gain
        W(:,:,k)=(Pxz(:,:,k)/Szz(:,:,k)')/Szz(:,:,k);

        %The innovation, transformed as necessary to keep values in a
        %desired range.
        innov(:,k)=innovTrans(z,zPred(:,:,k));

        %Updated state estimate
        xUpdate(:,k)=stateTrans(xPred(:,k)+W(:,:,k)*innov(:,k));

        %Updated state root covariance
        SUpdate(:,:,k)=tria([stateDiffTrans(xPredCenPoints(:,:,k)-W(:,:,k)*zPredCenPoints(:,:,k)),W(:,:,k)*SR]);
    end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%% Copyright (C) 2016-2022 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 csch (@var{x})
%% Symbolic csch function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = csch (x)
%%   @result{} y = (sym) csch(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = csch(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('csch', x);
end

%!error csch (sym(1), 2)
%!assert (isequaln (csch (sym(nan)), sym(nan)))

%!shared x, d
%! d = 1;
%! x = sym('1');

%!test
%! f1 = csch(x);
%! f2 = csch(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = csch(A);
%! f2 = csch(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = csch (d);
%! f = csch (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)

repo: octsympy


text:
%Arrival the blueprint for Arrival objects in GISMO
% An Arrival object is a container for phase arrival metadata
% See also Catalog
classdef Arrival
    properties(Dependent)
    %properties
        channelinfo
        time
        %arid
        %jdate
        iphase
        %deltim
        %azimuth
        %delaz
        %slow
        %delslo
        %ema
        %rect
        amp
        per
        %clip
        %fm
        signal2noise
        %qual
        %auth
    end
    properties
        waveforms
    end
    properties(Hidden)
        table
    end
    methods
        function obj = Arrival(sta, chan, time, iphase, varargin)
            % Parse required, optional and param-value pair arguments,
            % set default values, and add validation conditions
            p = inputParser;
            p.addRequired('sta', @iscell);
            p.addRequired('chan', @iscell);
            %p.addRequired('time', @(t) t>0 & t<now+1);
            p.addRequired('time', @isnumeric);
            p.addRequired('iphase', @iscell);
            p.addParameter('amp', NaN, @isnumeric);
            p.addParameter('per', NaN, @isnumeric);
            p.addParameter('signal2noise', NaN, @isnumeric);

            % Missed several properties out here just because of laziness.
            % Add them as needed.
            p.parse(sta, chan, time, iphase, varargin{:});
            fields = fieldnames(p.Results);
            for i=1:length(fields)
                field=fields{i};
                val = p.Results.(field);
                eval(sprintf('%s = val;',field));
            end
            ctag = ChannelTag.array('',sta,'',chan)';
            obj.table = table(time, datestr(time,26), datestr(time,'HH:MM'), datestr(time,'SS.FFF'), ctag.string(), iphase, amp, per, signal2noise, ...
                'VariableNames', {'time' 'date' 'hour_minute' 'second' 'channelinfo' 'iphase' 'amp' 'per' 'signal2noise'})
            obj.table = sortrows(obj.table, 'time', 'ascend');
            fprintf('\nGot %d arrivals\n',height(obj.table));edit
            misc_fields = struct;

        end

        function val = get.time(obj)
            val = obj.table.time;
        end

        function val = get.channelinfo(obj)
            val = obj.table.channelinfo;
        end

        function val = get.iphase(obj)
            val = obj.table.iphase;
        end

        function val = get.amp(obj)
            val = obj.table.amp;
        end

        function val = get.signal2noise(obj)
            val = obj.table.signal2noise;
        end

        function obj = set.amp(obj, amp)
            obj.table.amp = amp;
        end

        function summary(obj, showall)
        % ARRIVAL.SUMMARY Summarise Arrival object
            for c=1:numel(obj)
                obj(c)
                numarrs = height(obj(c).table)
                fprintf('Number of arrivals: %d\n',numarrs);
                if numarrs > 0
                    if ~exist('showall','var')
                            showall = false;
                    end
                    if numel(obj) == 1
                        if height(obj.table) <= 50 || showall
                            disp(obj.table)
                        else
                            disp(obj.table([1:50],:))

                            disp('* Only showing first 50 rows/arrivals - to see all rows/arrivals use:')
                            disp('*      arrivalObject.disp(true)')
                        end
                    end
                end
            end
        end

        function self2 = subset(self, columnname, findval)
            self2 = self;
            N = numel(self.time);
            indexes = [];
            if ~exist('findval','var')
                % assume columnname is actually row numbers
                indexes = columnname;
            else

                for c=1:N
                    gotval = eval(sprintf('self.%s(c);',columnname));
                    if isa(gotval,'cell')
                        gotval = cell2mat(gotval);
                    end
                    if isnumeric(gotval)
                        if gotval==findval
                            indexes = [indexes c];
                        end
                    else
                        if strcmp(gotval,findval)
                            indexes = [indexes c];
                        end
                    end
                end
            end
            self2.table = self.table(indexes,:);

%             % now go into misc_fields and apply same index subset to
%             % anything with N elements
%             fields = fieldnames(self.misc_fields);
%             for fieldnum=1:numel(fields)
%                 fieldval = getfield(self.misc_fields, fields{fieldnum});
%                 if numel(fieldval)==N
%                     self2.misc_fields = setfield(self2.misc_fields, fields{fieldnum}, fieldval(indexes));
%                 end
%             end

            self2.waveforms = self.waveforms(indexes);

        end

        % prototypes
        catalogobj = associate(self, maxTimeDiff)
        %arrivalobj = setminman(self, w, pretrig, posttrig, maxtimediff)
        arrivalobj = addmetrics(self, maxtimediff)
        arrivalobj = addwaveforms(self, datasourceobj, pretrigsecs, posttrigsecs);
    end
    methods(Static)
        function self = retrieve(dataformat, varargin)
        %ARRIVAL.RETRIEVE Read arrivals from common file formats & data sources.
        % retrieve can read phase arrivals from different earthquake catalog file
        % formats (e.g. Seisan, Antelope) and data sources (e.g. IRIS DMC) into a
        % GISMO Catalog object.
        %
        % Usage:
        %       arrivalObject = ARRIVAL.RETRIEVE(dataformat, 'param1', _value1_, ...
        %                                                   'paramN', _valueN_)
        %
        % dataformat may be:
        %
        %   * 'iris' (for IRIS DMC, using irisFetch.m),
        %   * 'antelope' (for a CSS3.0 Antelope/Datascope database)
        %   * 'seisan' (for a Seisan database with a REA/YYYY/MM/ directory structure)
        %   * 'zmap' (converts a Zmap data strcture to a Catalog object)
        %
        % See also CATALOG, IRISFETCH, CATALOG_COOKBOOK

        % Author: Glenn Thompson (glennthompson1971@gmail.com)

        %% To do:
        % Implement name-value parameter pairs for all methods
        % Test the Antelope method still works after factoring out db_load_origins
        % Test the Seisan method more
        % Add in support for 'get_arrivals'

            debug.printfunctionstack('>')
            self = [];
            switch lower(dataformat)
                case {'css3.0','antelope', 'datascope'}
                    if admin.antelope_exists()
                        try
                            self = Arrival.read_arrivals.antelope(varargin{:});
                        catch
                            % no arrivals
                        end
                    else
                        warning('Antelope toolbox for MATLAB not found')
                    end
                case 'hypoellipse'
                    self = read_hypoellipse(varargin{:});
                otherwise
                    self = NaN;
                    fprintf('format %s unknown\n\n',dataformat);
            end

            debug.printfunctionstack('<')
        end

        %cookbook()
    end
end

author: geoscience-community-codes

repo: GISMO


text:
function [g]= spm_gx_adem_cue(x,v,a,P)
% returns the prediction for cued responses (proprioception and vision)
% FORMAT [g]= spm_gx_adem_cue(x,v,a,P)
%
% x    - hidden states:
%   x.o  - intrinsic motor state (proprioceptive)
%
% v    - hidden causes
%
% P    - target locations (visual) - extrinsic coordinates (Cartesian)
%
% g    - sensations:
%   g.o  - motor angle (proprioception)
%   g.p  - finger location (visual)
%   g.c  - target contrast (visual)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_gx_adem_cue.m 4230 2011-03-07 20:58:38Z karl $

% evaluate positions in intrinsic (polar) coordinates
%--------------------------------------------------------------------------
g.o = x.o;
g.p = tan(x.o);
g.c = x.a*4;

author: spm


text:
function projections = findProjections(filePath,vecs,meanValues,pixels,parameters)
%findPosturalEigenmodes finds the projection of a set of images onto
%postural eigenmodes.
%
%   Input variables:
%
%       filePath -> cell array of VideoReader objects or a directory
%                       containing aligned .avi files
%       vecs -> postural eignmodes (L x (M<L) array)
%       meanValues -> mean value for each of the pixels
%       pixels -> radon-transform space pixels to use (Lx1 or 1xL array)
%       parameters -> struct containing non-default choices for parameters
%
%
%   Output variables:
%
%       projections -> N x d array of projection values
%
%
% (C) Gordon J. Berman, 2014
%     Princeton University

    if nargin < 5
        parameters = [];
    end
    parameters = setRunParameters(parameters);

    setup_parpool(parameters.numProcessors)

    %files = findAllImagesInFolders(filePath,'tiff');
    %N = length(files);

    if iscell(filePath)

        vidObjs = filePath;

    else

        files = findAllImagesInFolders(filePath,'avi');
        N = length(files);
        vidObjs = cell(N,1);
        parfor i=1:N
           vidObjs{i} = VideoReader(files{i});
        end

    end

    numThetas = parameters.num_Radon_Thetas;
    spacing = 180/numThetas;
    thetas = linspace(0,180-spacing,numThetas);
    scale = parameters.rescaleSize;
    numProjections = parameters.numProjections;
    batchSize = parameters.pca_batchSize;

    projections = find_PCA_projections(vidObjs,vecs(:,1:numProjections),...
        meanValues,pixels,thetas,numProjections,scale,batchSize);

    if parameters.numProcessors > 1  && parameters.closeMatPool
        close_parpool
    end

author: gordonberman

repo: MotionMapper


text:
function filename = fetch_url(filename)

% FETCH_URL checks the filename and downloads the file to a local copy in
% case it is specified as an Universal Resource Locator. It returns the
% name of the temporary file on the local filesystem.
%
% Use as
%   filename = fetch_url(filename)
%
% In case the filename does not specify an URL, it just returns the original
% filename.

% Copyright (C) 2012 Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if filetype_check_uri(filename, 'sftp')
  [user, host, filename] = filetype_check_uri(filename);
  [p, f, x] = fileparts(filename);
  p = tempdir;
  try
    mkdir(p);
    cmd = sprintf('sftp %s@%s:%s %s', user, host, filename, fullfile(p, [f x]));
    system(cmd);
    filename = fullfile(p, [f x]);
  end
  % elseif filetype_check_uri(filename, 'http')
  % FIXME the http scheme should be supported using default MATLAB
  % elseif filetype_check_uri(filename, 'ftp')
  % FIXME the http scheme should be supported using default MATLAB
  % elseif filetype_check_uri(filename, 'smb')
  % FIXME the smb scheme can be supported using smbclient
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


function  matcaffe_init(use_gpu, model_def_file, model_file)
% matcaffe_init(model_def_file, model_file, use_gpu)
% Initilize matcaffe wrapper

if nargin < 1
  % By default use CPU
  use_gpu = 0;
end
if nargin < 2 || isempty(model_def_file)
  % By default use imagenet_deploy
  model_def_file = '../../examples/imagenet/imagenet_deploy.prototxt';
end
if nargin < 3 || isempty(model_file)
  % By default use caffe reference model
  model_file = '../../examples/imagenet/caffe_reference_imagenet_model';
end

if caffe('is_initialized') == 0
  if exist(model_file, 'file') == 0
    % NOTE: you'll have to get the pre-trained ILSVRC network
    error('You need a network model file');
  end
  if ~exist(model_def_file,'file')
    % NOTE: you'll have to get network definition
    error('You need the network prototxt definition');
  end
  caffe('init', model_def_file, model_file)
end
fprintf('Done with init\n');

% set to use GPU or CPU
if use_gpu
  fprintf('Using GPU Mode\n');
  caffe('set_mode_gpu');
else
  fprintf('Using CPU Mode\n');
  caffe('set_mode_cpu');
end
fprintf('Done with set_mode\n');

% put into test mode
caffe('set_phase_test');
fprintf('Done with set_phase_test\n');


text:
function [c,G] = PAES(Problem,c,G,H,N,l_fails,l_opt,div)
% PAES local search

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    fails = 0;
    moves = 0;
    while fails < l_fails && moves < l_opt
        m = OperatorFEP(Problem,c);
        if all(c.obj<=m.obj)
            fails = fails + 1;
        else
            [H,dominated,~,mCrowd,cCrowd] = UpdateArchive(H,m,c,N,div);
            if all(c.obj>=m.obj)
                c = m;
                fails = 0;
            elseif ~dominated && mCrowd < cCrowd
                c = m;
            end
        end
        G     = UpdateArchive(G,m,SOLUTION(),N,div);
        moves = moves + 1;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef unbndfun < classicfun
%UNBNDFUN    Represent global functions on an unbounded interval [-inf inf] or
% a semi-infinite interval [-inf b] or [a inf].
%
% Constructor inputs:
%   UNBNDFUN(OP) constructs a UNBNDFUN object from the function handle OP on the
%   domain determined by the default in CHEBFUNPREF by mapping the domain to
%   [-1, 1] and constructing a ONEFUN object to represent the function
%   prescribed by OP.  OP should be vectorised (i.e., accept a vector input)
%   and output a vector of the same length as its input.
%
%   UNBNDFUN(OP, DATA) constructs an UNBNDFUN object using the data supplied in
%   the DATA structure.  DATA fields used by UNBNDFUN are:
%     DATA.DOMAIN       (Default:  Determined by CHEBFUNPREF)
%         A row vector with two elements in increasing order defining the
%         construction domain.  At least one element must be infinite.
%   In addition, UNBNDFUN may modify the following DATA fields before passing
%   them on to the ONEFUN constructor:
%     DATA.HSCALE       (Default:  Determined by DATA.DOMAIN)
%         Horizontal construction scale.
%     DATA.EXPONENTS    (Default:  Empty)
%         Exponents describing the behavior of the function at the endpoints of
%         the construction domain.  (See SINGFUN.)
%     DATA.SINGTYPE     (Default:  Empty)
%         Types of endpoint singularities to search for.  (See SINGFUN.)
%   If any fields in DATA are empty or not supplied, or if DATA itself is empty
%   or not supplied, appropriate default values are set.  Any fields in DATA
%   which are not recognized will be passed as-is to the ONEFUN constructor.
%
%   UNBNDFUN(OP, DATA, PREF) overrides the default behavior with that given by
%   the preferences in the structure or CHEBFUNPREF object PREF. See
%   CHEBFUNPREF for details.
%
% See also CLASSICFUN, CHEBFUNPREF, ONEFUN.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% UNBNDFUN Class Description:
%
% The UNBNDFUN class represents global functions on the infinite or
% semi-infinite interval [-inf, b], [a, inf], or [-inf, inf]. It achieves this
% by taking a onefun on [-1, 1] and applying a nonlinear mapping.
%
% Note that all binary UNBNDFUN operators (methods which can take two UNBNDFUN
% arguments) assume that the domains of the UNBNDFUN objects agree. The methods
% will not throw warnings if assumption is violated, but the results will not be
% meaningful under that circumstance.
%
% Class diagram: [<<CLASSICFUN>>] <>-- [<<onefun>>]
%                       ^
%                       |
%                   [unbndfun]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function obj = unbndfun(op, data, pref)

            % Parse inputs.
            if ( (nargin < 1) || isempty(op) )
                obj.domain = [];
                obj.mapping = [];
                obj.onefun = [];
                return
            end

            if ( (nargin < 2) || isempty(data) )
                data = struct();
            end

            if ( (nargin < 3) || isempty(pref) )
                pref = chebfunpref();
            else
                pref = chebfunpref(pref);
            end

            data = parseDataInputs(data, pref);

            % Check domain
            if ( ~all(size(data.domain) == [1, 2]) || (diff(data.domain) <= 0) )
                error('CHEBFUN:UNBNDFUN:unbndfun:domain',...
                    ['Domain argument should be a row vector with two ', ...
                    'entries in increasing order.']);
            elseif ( ~any(isinf(data.domain)) )
                error('CHEBFUN:UNBNDFUN:unbndfun:boundedDomain',...
                    'Domain argument should be unbounded.');
            end

            % Remap the OP to be a function on [-1, 1].
            unbndmap = unbndfun.createMap(data.domain);
            if ( isa(op, 'function_handle') )
                op = @(x) op(unbndmap.For(x));
            elseif ( isnumeric(op) )
                if ( ~any(op(:)) )
                    op = @(x) zeros(length(x), size(op, 2));
                else
                    %[TODO]: Implement this.
                    error('CHEBFUN:UNBNDFUN:unbndfun:inputValues',...
                        ['UNBNDFUN does not support non-zero construction ' ...
                         'from values.']);
                end
            end

            % Deal with exponents for singular functions.
            if ( isempty(data.exponents) )
                % The user hasn't supplied exponents, but functions on
                % unbounded domains are often singular once remapped to [-1,
                % 1].  Try to detect this by seeing if the function is infinite
                % at either endpoint and, if so, enable singularity detection.
                lVal = feval(op, -1);
                rVal = feval(op, 1);
                if ( any(isinf([lVal rVal])) )
                    pref.blowup = true;
                    singType = pref.blowupPrefs.defaultSingType;
                    data.singType = {singType, singType};
                end
            else
                % Remapping to [-1, 1] negates exponents, which are given with
                % respect to the function on the infinite domain.
                ind = isinf(data.domain);
                pref.blowup = true;
                data.exponents(ind) = -data.exponents(ind);
            end

            % Call the ONEFUN constructor:
            obj.onefun = onefun.constructor(op, data, pref);

            % Add the domain and mapping:
            obj.domain = data.domain;
            obj.mapping = unbndmap;
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )

        % Retrieve and modify preferences for this class.
        prefs = pref(varargin);

        % Noninear map from [-1, 1] to the domain of the UNBNDFUN.
        m = createMap(domain);

        % Make a UNBNDFUN (constructor shortcut):
        f = make(varargin);

    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PRIVATE METHODS IMPLEMENTED IN THIE FILE:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function data = parseDataInputs(data, pref)
%PARSEDATAINPUTS   Parse inputs from the DATA structure and assign defaults.

if ( ~isfield(data, 'domain') || isempty(data.domain) )
    data.domain = pref.domain;
end

if ( ~isfield(data, 'hscale') || isempty(data.hscale) )
    % TODO:  Why is the hscale of an unbounded domain always 1?
    data.hscale = 1;
end

if ( ~isfield(data, 'exponents') || isempty(data.exponents) )
    data.exponents = [];
end

if ( ~isfield(data, 'singType') || isempty(data.singType) )
    data.singType = [];
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ net,opts ] = selective_sgd( net,opts )
%NET_APPLY_SELECTIVE_SGD Summary of this function goes here
%   Detailed explanation goes here

    if (mod(opts.parameters.current_ep,opts.parameters.ssgd_search_freq)==1||opts.parameters.ssgd_search_freq==1)

        if ~isfield(opts.parameters,'selection_count')
            opts.parameters.selection_count=0;
        end

        [lr_best]=select_learning_rate(net,opts);
        opts.parameters.lr=lr_best;

        if ~isfield(opts.parameters,'selected_lr')
            opts.parameters.selected_lr(1)=lr_best;
        else
            opts.parameters.selected_lr(end+1)=lr_best;
        end

        opts.parameters.selection_count=opts.parameters.selection_count+1;

        disp(['Selected learning rate: ',num2str(opts.parameters.lr)]);

    end

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
%  Save NIFTI header extension.
%
%  Usage: save_nii_ext(ext, fid)
%
%  ext - struct with NIFTI header extension fields.
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function save_nii_ext(ext, fid)

   if ~exist('ext','var') || ~exist('fid','var')
      error('Usage: save_nii_ext(ext, fid)');
   end

   if ~isfield(ext,'extension') || ~isfield(ext,'section') || ~isfield(ext,'num_ext')
      error('Wrong header extension');
   end

   write_ext(ext, fid);

   return;                                      % save_nii_ext

%---------------------------------------------------------------------
function write_ext(ext, fid)

   fwrite(fid, ext.extension, 'uchar');

   for i=1:ext.num_ext
      fwrite(fid, ext.section(i).esize, 'int32');
      fwrite(fid, ext.section(i).ecode, 'int32');
      fwrite(fid, ext.section(i).edata, 'uchar');
   end

   return;                                      % write_ext

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%%*****************************************************************************
%% HSDsqlp: solve an semidefinite-quadratic-linear program
%%   by infeasible path-following method on the homogeneous self-dual model.
%%
%%  [obj,X,y,Z,info,runhist] =
%%      HSDsqlp(blk,At,C,b,OPTIONS,X0,y0,Z0);
%%
%%  Input: blk: a cell array describing the block diagonal structure of SQL data.
%%          At: a cell array with At{p} = [svec(Ap1) ... svec(Apm)]
%%         b,C: data for the SQL instance.
%%     OPTIONS: a structure that specifies parameters required in HSDsqlp.m,
%%              (if it is not given, the default in sqlparameters.m is used).
%%
%%  (X0,y0,Z0): an initial iterate (if it is not given, the default is used).
%%
%%  Output: obj  = [<C,X> <b,y>].
%%          (X,y,Z): an approximately optimal solution or a primal or dual
%%                   infeasibility certificate.
%%          info.termcode = termination-code
%%          info.iter     = number of iterations
%%          info.obj      = [primal-obj, dual-obj]
%%          info.cputime  = total-time
%%          info.gap      = gap
%%          info.pinfeas  = primal_infeas
%%          info.dinfeas  = dual_infeas
%%          runhist.pobj    = history of primal objective value.
%%          runhist.dobj    = history of dual   objective value.
%%          runhist.gap     = history of <X,Z>.
%%          runhist.pinfeas = history of primal infeasibility.
%%          runhist.dinfeas = history of dual   infeasibility.
%%          runhist.cputime = history of cputime spent.
%%-----------------------------------------------------------------
%%  The OPTIONS structure specifies the required parameters:
%%      vers  gam  predcorr  expon  gaptol  inftol  steptol
%%      maxit  printlevel  ...
%%      (all have default values set in sqlparameters.m).
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function [obj,X,y,Z,info,runhist] = ...
            HSDsqlp(blk,At,C,b,OPTIONS,X0,y0,Z0,kap0,tau0,theta0);

    if (nargin < 5); OPTIONS = []; end

    isemptyAtb = 0;
    if isempty(At) & isempty(b);
       %% Add redundant constraint: <-I,X> <= 0
       b = 0;
       At = ops(ops(blk,'identity'),'*',-1);
       numblk = size(blk,1);
       blk{numblk+1,1} = 'l'; blk{numblk+1,2} = 1;
       At{numblk+1,1} = 1; C{numblk+1,1} = 0;
       isemptyAtb = 1;
    end
%%
%%-----------------------------------------
%% get parameters from the OPTIONS structure.
%%-----------------------------------------
%%
   warning off;

   matlabversion = sscanf(version,'%f');
   if strcmp(computer,'PCWIN64') | strcmp(computer,'GLNXA64')
      par.computer = 64;
   else
      par.computer = 32;
   end
   par.matlabversion = matlabversion(1);
   par.vers        = 0;
   par.predcorr    = 1;
   par.gam         = 0;
   par.expon       = 1;
   par.gaptol      = 1e-8;
   par.inftol      = 1e-8;
   par.steptol     = 1e-6;
   par.maxit       = 100;
   par.printlevel  = 3;
   par.stoplevel   = 1;
   par.scale_data  = 0;
   par.spdensity   = 0.4;
   par.rmdepconstr = 0;
   par.smallblkdim = 50;
   par.schurfun     = cell(size(blk,1),1);
   par.schurfun_par = cell(size(blk,1),1);
%%
   parbarrier  = cell(size(blk,1),1);
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      if strcmp(pblk{1},'s') | strcmp(pblk{1},'q')
         parbarrier{p} = zeros(1,length(pblk{2}));
      elseif strcmp(pblk{1},'l') | strcmp(pblk{1},'u' )
         parbarrier{p} = zeros(1,sum(pblk{2}));
      end
   end
   parbarrier_0 = parbarrier;
%%
   if exist('OPTIONS')
      if isfield(OPTIONS,'vers');        par.vers     = OPTIONS.vers; end
      if isfield(OPTIONS,'predcorr');    par.predcorr = OPTIONS.predcorr; end
      if isfield(OPTIONS,'gam');         par.gam      = OPTIONS.gam; end
      if isfield(OPTIONS,'expon');       par.expon    = OPTIONS.expon; end
      if isfield(OPTIONS,'gaptol');      par.gaptol   = OPTIONS.gaptol; end
      if isfield(OPTIONS,'inftol');      par.inftol   = OPTIONS.inftol; end
      if isfield(OPTIONS,'steptol');     par.steptol  = OPTIONS.steptol; end
      if isfield(OPTIONS,'maxit');       par.maxit    = OPTIONS.maxit; end
      if isfield(OPTIONS,'printlevel');  par.printlevel  = OPTIONS.printlevel; end
      if isfield(OPTIONS,'stoplevel');   par.stoplevel   = OPTIONS.stoplevel; end
      if isfield(OPTIONS,'scale_data');  par.scale_data  = OPTIONS.scale_data; end
      if isfield(OPTIONS,'spdensity');   par.spdensity   = OPTIONS.spdensity; end
      if isfield(OPTIONS,'rmdepconstr'); par.rmdepconstr = OPTIONS.rmdepconstr; end
      if isfield(OPTIONS,'smallblkdim'); par.smallblkdim = OPTIONS.smallblkdim; end
      if isfield(OPTIONS,'parbarrier');
         parbarrier = OPTIONS.parbarrier;
         if isempty(parbarrier); parbarrier = parbarrier_0; end
         if ~iscell(parbarrier);
            tmp = parbarrier; clear parbarrier; parbarrier{1} = tmp;
         end
         if (length(parbarrier) < size(blk,1))
            len = length(parbarrier);
            parbarrier(len+1:size(blk,1)) = parbarrier_0(len+1:size(blk,1));
         end
      end
      if isfield(OPTIONS,'schurfun');
         par.schurfun = OPTIONS.schurfun;
         if ~isempty(par.schurfun); par.scale_data = 0; end
      end
      if isfield(OPTIONS,'schurfun_par'); par.schurfun_par = OPTIONS.schurfun_par; end
      if isempty(par.schurfun);     par.schurfun = cell(size(blk,1),1); end
      if isempty(par.schurfun_par); par.schurfun_par = cell(size(blk,1),1); end
   end
   if (size(blk,2) > 2); par.smallblkdim = 0; end
%%
%%-----------------------------------------
%% convert matrices to cell arrays.
%%-----------------------------------------
%%
   if ~iscell(At); At = {At}; end;
   if ~iscell(C);  C = {C}; end;
   if all(size(At) == [size(blk,1), length(b)]);
      convertyes = zeros(size(blk,1),1);
      for p = 1:size(blk,1)
         if strcmp(blk{p,1},'s') & all(size(At{p,1}) == sum(blk{p,2}))
            convertyes(p) = 1;
         end
      end
      if any(convertyes)
         if (par.printlevel);
            fprintf('\n sqlp: converting At into required format');
         end
         At = svec(blk,At,ones(size(blk,1),1));
      end
   end
%%
%%-----------------------------------------
%% validate SQLP data.
%%-----------------------------------------
%%
   tstart = cputime;
   [blk,At,C,b,blkdim,numblk] = validate(blk,At,C,b,par);
   [blk,At,C,b,iscmp] = convertcmpsdp(blk,At,C,b);
   if (iscmp) & (par.printlevel>=2)
      fprintf('\n SQLP has complex data');
   end
   if (nargin <= 5) | (isempty(X0) | isempty(y0) | isempty(Z0));
      if (max([ops(At,'norm'),ops(C,'norm'),norm(b)]) > 1e2)
         [X0,y0,Z0] = infeaspt(blk,At,C,b,1);
      else
         [X0,y0,Z0] = infeaspt(blk,At,C,b,2,1);
      end
   end
   if ~iscell(X0);  X0 = {X0}; end;
   if ~iscell(Z0);  Z0 = {Z0}; end;
   if (length(y0) ~= length(b))
      error('HSDsqlp: length of b and y0 not compatible');
   end
   [X0,Z0] = validate_startpoint(blk,X0,Z0,par.spdensity);
%%
   if (nargin <= 8) | (isempty(kap0) | isempty(tau0) | isempty(theta0))
      if (max([ops(At,'norm'),ops(C,'norm'),norm(b)]) > 1e6)
         kap0 = 10*blktrace(blk,X0,Z0);
      else
         kap0 = blktrace(blk,X0,Z0);
      end
      tau0 = 1; theta0 = 1;
   end
%%
   if (par.printlevel>=2)
      fprintf('\n num. of constraints = %2.0d',length(b));
      if blkdim(1);
         fprintf('\n dim. of sdp    var  = %2.0d,',blkdim(1));
         fprintf('   num. of sdp  blk  = %2.0d',numblk(1));
      end
      if blkdim(2);
         fprintf('\n dim. of socp   var  = %2.0d,',blkdim(2));
         fprintf('   num. of socp blk  = %2.0d',numblk(2));
      end
      if blkdim(3); fprintf('\n dim. of linear var  = %2.0d',blkdim(3)); end
      if blkdim(4); fprintf('\n dim. of free   var  = %2.0d',blkdim(4)); end
   end
%%
%%-----------------------------------------
%% detect unrestricted blocks in linear blocks
%%-----------------------------------------
%%
   user_supplied_schurfun = 0;
   for p = 1:size(blk,1)
      if ~isempty(par.schurfun{p}); user_supplied_schurfun = 1; end
   end
   if (user_supplied_schurfun == 0)
      [blk2,At2,C2,ublkinfo,parbarrier2,X02,Z02] = ...
      detect_ublk(blk,At,C,parbarrier,X0,Z0,par.printlevel);
   else
      blk2 = blk; At2 = At; C2 = C;
      parbarrier2 = parbarrier; X02 = X0; Z02 = Z0;
      ublkinfo = cell(size(blk2,1),1);
   end
   ublksize = blkdim(4);
   for p = 1:size(ublkinfo,1)
      ublksize = ublksize + length(ublkinfo{p});
   end
%%
%%-----------------------------------------
%% detect diagonal blocks in semidefinite blocks
%%-----------------------------------------
%%
   if (user_supplied_schurfun==0)
      [blk3,At3,C3,diagblkinfo,diagblkchange,parbarrier3,X03,Z03] = ...
      detect_lblk(blk2,At2,C2,b,parbarrier2,X02,Z02,par.printlevel);
   else
      blk3 = blk2; At3 = At2; C3 = C2;
      parbarrier3 = parbarrier2; X03 = X02; Z03 = Z02;
      diagblkchange = 0;
      diagblkinfo = cell(size(blk3,1),1);
   end
%%
%%-----------------------------------------
%% main solver
%%-----------------------------------------
%%
   exist_analytic_term = 0;
   for p = 1:size(blk3,1);
      idx = find(parbarrier3{p} > 0);
      if ~isempty(idx); exist_analytic_term = 1; end
   end
%%
   if (par.vers == 0);
      if blkdim(1); par.vers = 1; else; par.vers = 2; end
   end
   par.blkdim = blkdim;
   par.ublksize = ublksize;
   [obj,X3,y,Z3,info,runhist] = ...
        HSDsqlpmain(blk3,At3,C3,b,par,X03,y0,Z03,kap0,tau0,theta0);
%%
%%-----------------------------------------
%% recover semidefinite blocks from linear blocks
%%-----------------------------------------
%%
   if any(diagblkchange)
      X2 = cell(size(blk2,1),1); Z2 = cell(size(blk2,1),1);
      count = 0;
      for p = 1:size(blk2,1)
         pblk = blk2(p,:);
         n = sum(pblk{2});
         blkno      = diagblkinfo{p,1};
         idxdiag    = diagblkinfo{p,2};
         idxnondiag = diagblkinfo{p,3};
         if ~isempty(idxdiag)
            len = length(idxdiag);
            Xtmp = [idxdiag,idxdiag,X3{end}(count+[1:len]); n, n, 0];
            Ztmp = [idxdiag,idxdiag,Z3{end}(count+[1:len]); n, n, 0];
            if ~isempty(idxnondiag)
               [ii,jj,vv] = find(X3{blkno});
               Xtmp = [Xtmp; idxnondiag(ii),idxnondiag(jj),vv];
               [ii,jj,vv] = find(Z3{blkno});
               Ztmp = [Ztmp; idxnondiag(ii),idxnondiag(jj),vv];
            end
            X2{p} = spconvert(Xtmp);
            Z2{p} = spconvert(Ztmp);
	    count = count + len;
         else
	    X2(p) = X3(blkno); Z2(p) = Z3(blkno);
         end
      end
   else
      X2 = X3; Z2 = Z3;
   end
%%
%%-----------------------------------------
%% recover linear block from unrestricted block
%%-----------------------------------------
%%
   numblk = size(blk,1);
   numblknew = numblk;
   X = cell(numblk,1); Z = cell(numblk,1);
   for p = 1:numblk
      n = blk{p,2};
      if isempty(ublkinfo{p,1})
         X{p} = X2{p}; Z{p} = Z2{p};
      else
	 Xtmp = zeros(n,1); Ztmp = zeros(n,1);
         Xtmp(ublkinfo{p,1}) = max(0,X2{p});
         Xtmp(ublkinfo{p,2}) = max(0,-X2{p});
         Ztmp(ublkinfo{p,1}) = max(0,Z2{p});
         Ztmp(ublkinfo{p,2}) = max(0,-Z2{p});
         if ~isempty(ublkinfo{p,3})
            numblknew = numblknew + 1;
            Xtmp(ublkinfo{p,3}) = X2{numblknew};
            Ztmp(ublkinfo{p,3}) = Z2{numblknew};
         end
	 X{p} = Xtmp; Z{p} = Ztmp;
      end
   end
%%
%%-----------------------------------------
%% recover complex solution
%%-----------------------------------------
%%
   if (iscmp)
      for p = 1:numblk
         pblk = blk(p,:);
         n = sum(pblk{2})/2;
         if strcmp(pblk{1},'s');
            X{p} = X{p}(1:n,1:n) + sqrt(-1)*X{p}(n+[1:n],1:n);
            Z{p} = Z{p}(1:n,1:n) + sqrt(-1)*Z{p}(n+[1:n],1:n);
            X{p} = 0.5*(X{p}+X{p}');
            Z{p} = 0.5*(Z{p}+Z{p}');
         end
      end
   end
   if (isemptyAtb)
      X = X(1:end-1); Z = Z(1:end-1);
   end
%%*****************************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
% DEMO2 - Subsetting ncdataset data
echo('on')
% STARTING DEMO2 ----------------------------------------------------------
% An example of subsetting data using ncdataset

%% ---- Open the dataset
ds = ncdataset('http://geoport.whoi.edu/thredds/dodsC/examples/OS_M1_20081008_TS.nc');

% You can view the variables available to you
ds.variables

%% Plot all the data
figure;
plot(ds.time('TIME'), ds.data('TEMP', [1 1 1 1], [max(ds.size('TEMP')) 1 1 1], [1 1 1 1]))
hold('on')

%% ---- Lets fetch a subset of time in Matlab's native format
startIdx = 100;
endIdx = max(ds.size('TIME'));
stride = 10;
t = ds.data('TIME', startIdx, endIdx, stride);
t = ds.time('TIME', t); % Convert time data to matlab format. See help ncdataset.time

%% ---- Now lets get a subset of the temperature data.
% NOTE: The shape of the variables size is important for subsetting
ds.size('TEMP')
% Use variable, start, end, stride to subset
temp = ds.data('TEMP', [startIdx 1 1 1], [endIdx 1 1 1], [stride 1 1 1]);

%% ---- Add Subsetted Data to Plot
plot(t, temp, 'r.');...
datetick('x', 2);grid;...
legend('All Data', 'Decimated Data');...
title({'Surface Temperature at M1 Mooring in Monterey Bay',ds.attribute('references'),ds.location},'interpreter','none');...
ylabel('Temperature [^oC]');
hold('off');
shg

echo('off') % ENDING DEMO2 ------------------------------------------------

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function [theta_gibbs,sigma_gibbs,sigmatilde_gibbs,sig_gibbs]=panel5gibbs(y,Y,Xtilde,Xdot,N,n,T,d,theta0,Theta0,alpha0,delta0,It,Bu,pick,pickf)

% compute first  preliminary elements
% compute alphabar
alphabar=N*n*T+alpha0;
% compute the inverse Theta0
invTheta0=sparse(diag(1./diag(Theta0)));
% initiate the Gibbs sampler
% initiate the counting of iterations
count=1;
pickcount=1;
% initiate the record matrices
theta_gibbs=zeros(d,It-Bu);
sigmatilde_gibbs=zeros((N*n)^2,It-Bu);
sig_gibbs=zeros(1,It-Bu);
sigma_gibbs=zeros((N*n)^2,It-Bu);

% step 1: compute initial values
% initial value for theta (use OLS values)
theta=(Xtilde*Xtilde')\(Xtilde*y);
% initial value for sigmatilde (use residuals form OLS values)
eps=y-Xtilde'*theta;
eps=reshape(eps,T,N*n);
sigmatilde=eps'*eps;
% initiate value for the sigma, the scaling term for the errors
sig=1;
% initiate value for the matrix sigma, the residual variance-covariance matrix
sigma=sig*sigmatilde;
% initiate value for eyesigma
% compute the inverse of sigma
C=bear.trns(chol(bear.nspd(sigma),'Lower'));
invC=C\speye(N*n);
invsigma=invC*invC';
% then compute eyesigma
eyesigma=kron(speye(T),invsigma);
% finally, initiate eyetheta
eyetheta=kron(speye(T),theta);

hbar = bear.parfor_progressbar(It-Bu,'Progress of Panel BVAR Gibbs Sampler.');  %create the progress bar

% run the Gibbs sampler
while count<=It

% step 2: obtain sigmatilde
% compute Sbar
Sbar=(1/sig)*(Y-Xdot*eyetheta)*(Y-Xdot*eyetheta)';
sigmatilde=bear.iwdraw(Sbar,T);

% step 3: obtain sig
% compute the inverse of sigmatilde
C=bear.trns(chol(bear.nspd(sigmatilde),'Lower'));
invC=C\speye(N*n);
invsigmatilde=invC*invC';
% compute deltabar
deltabar=trace((Y-Xdot*eyetheta)*(Y-Xdot*eyetheta)'*invsigmatilde)+delta0;
% draw sig
sig=bear.igrandn(alphabar/2,deltabar/2);

% step 4: compute sigma and eyesigma
sigma=sig*sigmatilde;
C=bear.trns(chol(bear.nspd(sigma),'Lower'));
invC=C\speye(N*n);
invsigma=invC*invC';
eyesigma=kron(speye(T),invsigma);

% step 5: obtain theta
% compute Thetabar
invThetabar=full((Xtilde*eyesigma*Xtilde'+invTheta0));
C=bear.trns(chol(bear.nspd(invThetabar),'Lower'));
invC=C\speye(d);
Thetabar=invC*invC';
% compute thetabar
thetabar=Thetabar*(Xtilde*eyesigma*y+invTheta0*theta0);
% draw theta
theta=thetabar+chol(bear.nspd(Thetabar),'lower')*mvnrnd(zeros(d,1),eye(d))';

% step 6: obtain eyetheta
eyetheta=kron(speye(T),theta);

   % record phase
   % if the burn-in sample phase is not yet over
   if count<=Bu
   % simply add 1 to the iteration count
   count=count+1;
   % on the other hand, if the burn-in sample phase is over
   elseif count>Bu
   % adding one iteration to the count will depend on wether post-burn selection applies
      % if there is no post burn selection
      if pick==0
      % record the draw
      theta_gibbs(:,count-Bu)=theta;
      sigmatilde_gibbs(:,count-Bu)=bear.vec(sigmatilde);
      sig_gibbs(1,count-Bu)=sig;
      sigma_gibbs(:,count-Bu)=bear.vec(sigma);
      % and add one to the count
      count=count+1;
      % if there is post burn selection, only one draw over 'fpick' draws will be retained
      elseif pick==1
         % if the iteration does not correspond to fpick, don't record the results, don't increase the regular count, but do increase pickcount by 1
         if pickcount~=pickf
         pickcount=pickcount+1;
         % on the other hand, if the iteration does correspond to fpick
         elseif pickcount==pickf
         % do record the results
         theta_gibbs(:,count-Bu)=theta;
         sigmatilde_gibbs(:,count-Bu)=bear.vec(sigmatilde);
         sig_gibbs(1,count-Bu)=sig;
         sigma_gibbs(:,count-Bu)=bear.vec(sigma);
         % then increase the regular count by 1 and re-initialise pickcount
         count=count+1;
         pickcount=1;
         end
      end
   end

   hbar.iterate(1);   % update progress by one iteration

end

close(hbar);   %close progress bar

author: european-central-bank

repo: BEAR-toolbox


text:
% eeglab2eegtoolbox() - convert EEGLAB structure to EEG toolbox
%                       structure (for ERP peak picking mainly).
%
% Usage:
%   >> p = eeglab2eegtoolbox( EEG );
%
% Inputs:
%   EEG - EEG structure
%
% Outputs:
%   p   - EEG toolbox structure
%
% Author: Arnaud Delorme, SCCN, INC, UCSD
%
% See also: eeglab()

% Copyright (C) 2006 Arnaud Delorme, Salk Institute, arno@salk.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 p = eeglab2eegtoolbox( EEG );

    if nargin < 1
        help eeglab2eegtoolbox;
        return;
    end;

    % default structure
    % -----------------
    p = eeg_toolbox_defaults;

    % exporting importing electrodes
    % ------------------------------
    fid = fopen('tmp.txt', 'w');
    % fiducials should also be added here
    for index = 1:length(EEG.chanlocs)
        tmp = EEG.chanlocs(index);
        fprintf(fid, '%s\t%d\t%f\t%f\t%f\n', tmp.labels, 69, tmp.X, tmp.Y, tmp.Z);
    end;
    fprintf(fid, 'Centroid\t%d\t%f\t%f\t%f\n\n', 99, 0, 0, 0);
    fclose(fid);

    p.elec.path = pwd;
    p.elec.file = 'tmp.txt';
    p.elec.n    = length(EEG.chanlocs)+1;
    p = elec_open(p);
    delete(fullfile(pwd, 'tmp.txt'));

    % importing data
    % --------------
    p.volt.data = double(mean(EEG.data,3)');
    p.volt.timeArray = EEG.times';
    p.volt.points    = EEG.pnts;
    p.volt.sampleHz  = EEG.srate;
    p.volt.sampleTime = EEG.srate;
    p.volt.sampleMsec = 1000/EEG.srate;
    p.volt.channels   = length(EEG.chanlocs);
    p.volt.epochStart = EEG.xmin*1000;
    p.volt.epochEnd   = EEG.xmax*1000;

author: PatternRecognition

repo: OpenBMI


text:
function n = cfn_e_size_total ( m, ell_max )

%*****************************************************************************80
%
%% CFN_E_SIZE_TOTAL: Closed Fully Nested, Exponential Growth.
%
%  Discussion:
%
%    This calculation assumes that an exponential growth rule is being used,
%    that is, that the 1D rules have orders 1, 3, 7, 15, 31, and so on.
%
%    It counts the number of points in all the product rules that compose
%    the sparse grid, and it does not reduce the count in any way to account
%    for repeated points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Fabio Nobile, Raul Tempone, Clayton Webster,
%    A Sparse Grid Stochastic Collocation Method for Partial Differential
%    Equations with Random Input Data,
%    SIAM Journal on Numerical Analysis,
%    Volume 46, Number 5, 2008, pages 2309-2345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer ELL_MAX, the sparse grid level.
%
%    Output, integer N, the total number of points in the grids.
%
  n = 0;

  lvec = zeros ( m, 1 );

  ell_min = max ( ell_max + 1 - m, 0 );

  for ell = ell_min : ell_max

    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ lvec, more, h, t ] = comp_next ( ell, m, lvec, more, h, t );

      nvec = cfn_e_nvec_from_lvec ( lvec );
      ell_num = prod ( nvec(1:m) );
      n = n + ell_num;

      if ( ~ more )
        break
      end

    end
  end

  return
end


text:
function variance = normal_01_variance ( )

%*****************************************************************************80
%
%% NORMAL_01_VARIANCE returns the variance of the Normal 01 PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = 1.0;

  return
end


text:
%function feat = rgb_fun(filename, fun)
% apply function handler fun separately to all rgb channel of image located
% in filename
function feat = rgb_fun(filename, fun)
	rgbx = imread(filename);
	feat = [];
	if (numel(size(rgbx)) == 3)
		for c = 1:3
			cx = single(squeeze(rgbx(:,:,c)));
			feat = [feat, fun(cx)];
		end
	else
		cx = single(squeeze(rgbx(:,:)));
		feat = repmat(fun(cx),[1,3]);
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: SSD versus rotations, linearInter, HNSP level=8
%
%==============================================================================

clear, close all, help(mfilename);

setup2DHNSPData;
imgModel('set','imgModel','linearInter');
level = 8; omega = ML{level}.omega; m = ML{level}.m;
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega);
xc    = getCellCenteredGrid(omega,m);
Rc    = imgModel(R,omega,xc);

center = (omega(2:2:end)-omega(1:2:end))'/2;
trafo('set','trafo','rotation2D','c',center);

wc = pi/2*linspace(-1,1,101);  dc = zeros(size(wc));
figure(1); clf;
for j=1:length(wc),
  yc = trafo(wc(j),xc);
  Tc = imgModel(T,omega,yc);
  dc(j) = SSD(Tc,Rc,omega,m);
  viewImage(255-abs(Tc-Rc),omega,m); drawnow; FAIRpause(1/60)
end;
figure(2); clf; p1 = plot(wc,dc);
%==============================================================================


text:
%% Import data from text file.
% Script for importing data from the following text file:
%
%    C:\ArizCU12.5m.dat
%
% To extend the code to different selected data or a different text file,
% generate a function instead of a script.

% Auto-generated by MATLAB on 2013/01/03 09:36:00
% David Willingham February 2013
% Copyright 2013 The MathWorks, Inc
%% Initialize variables.
filename = 'ArizCU12.5m.dat';
delimiter = '\t';
startRow = 7;

%% Format string for each line of text:
%   column1: double (%f)
%	column2: double (%f)
%   column3: double (%f)
%	column4: double (%f)
% For more information, see the TEXTSCAN documentation.
formatSpec = '%f%f%f%f%[^\n\r]';

%% Open the text file.
fileID = fopen(filename,'r');

%% Read columns of data according to format string.
% This call is based on the structure of the file used to generate this
% code. If an error occurs for a different file, try regenerating the code
% from the Import Tool.
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'HeaderLines' ,startRow-1, 'ReturnOnError', false);

%% Close the text file.
fclose(fileID);

%% Post processing for unimportable data.
% No unimportable data rules were applied during the import, so no post
% processing code is included. To generate code which works for
% unimportable data, select unimportable cells in a file and regenerate the
% script.

%% Allocate imported array to column variable names
x = dataArray{:, 1};
y = dataArray{:, 2};
z = dataArray{:, 3};
c = dataArray{:, 4};

%% Clear temporary variables
clearvars filename delimiter startRow formatSpec fileID dataArray ans;


text:
function [obj, rawDb, filteredDb] =  get_antelope_traces(obj)
   % GET_ANTELOPE_TRACE gets interesting data from database, and returns the tracebuf object
   % [tr(1:end)] =  get_antelope_trace(startdates, enddates, criteriaList, database)
   %    STARTDATE is a matlab datenum
   %    ENDDATE is also in matlab datenum
   %    CRITERIALIST is a cell of statements used to filter the antelope
   %        database.  example: {'sta=="OKCF"','chan=="BHZ"'}
   %    DATABASE can either be a database name, or an open antelope database
   %    pointer.
   %
   %    TR is a tracebuf object containing the data requested.
   %
   %  In this usage, with a single output argument, TR, the database is opened
   %  and then closed before the function returns.
   %
   % You can opt to keep the database open, which can make subsequent calls to
   % get_antelope_trace much faster.  To do this, ask for the database pointer
   % as one of the output arguments.  The database will be closed by any
   % function call that uses the database pointer as the 'database' parameter,
   % and does not ask for a database pointer back.
   %
   % [tr, mydb] = get_antelope_trace(...)
   %    if mydb is asked for, then the database will not be closed, allowing for
   %    more efficient reuse.  mydb is the unfiltered (and open) wfdisc table
   %
   % [tr, mydb, filteredDB] = get_antelope_trace(...)
   %   returns wfdisc table in mydb, and a filtered wfdisc table in filteredDB
   %
   % if no records are found, then tr will be set to []
   %
   % will generate error message if no records found, so consider using in a
   % TRY-CATCH loop.
   %
   %
   % Example:
   %  % set up the data
   %  sDate = datenum('7/12/2008 05:00:00.00');
   %  eDate = datenum('7/12/2008 05:10:00.00');
   %  stationsToGet = {'OKFG','OKSO','OMG'};
   %  channelsToGet = {'BHZ', 'BHE'};
   %  mydb = '/mydir/my_antelope_database';
   %
   %  % get our criteria together
   %  for st = 1:numel(stationsToGet)
   %    staCriteria(st) = {sprintf('sta=="%s"',stationsToGet(st)};
   %  end
   %  for ch = 1:numel(channelsToGet)
   %    chanCriteria(ch) = {sprintf('chan=="%s"',channelsToGet(ch)};
   %  end
   %
   %  % open the database once, and read in each tracebuf. Then translate into
   %  % waveforms.
   %  n = 1;
   %  for sta = 1:numel(staCriteria)
   %    for ch = 1:numel(chanCriteria)
   %      critera = [chanCriteria(ch),staCriteria(sta)];
   %      [tr, mydb] = get_antelope_trace(sDate,eDate, criteria, mydb)
   %      w(n) = traceToWaveform(tr);
   %      n = n+1;
   %    end
   %  end
   %  % finally, close up shop.
   %  dbclose(mydb)

   % Modifications
   %%% Glenn Thompson 2012/02/06: Occasionally the C program trload_css cannot even load the trace data. Added try...catch..end to handle this.

   useExistingDatabasePtr =  obj.pointsToAntelopeDatabase;

   if useExistingDatabasePtr
      try
         dbnrecs(obj.dbpointer);
      catch
         warning('Waveform:load_antelope:databaseNotOpen', ...
            'a Database Pointer was passed to trace, but the database was not open');
         %tr = trnew; %return a new object, forcing the ability to destroy it later
         obj.trpointer = { -1 };
         return;
      end
   end

   % do not close the database if a pointer is asked for in the return arguments
   onFinishCloseDB = nargout < 2;

   antelope_starts = mep2dep(obj.startdate);
   antelope_ends = mep2dep(obj.enddates);

   if ~obj.isopen()
      obj = obj.open();
   end

   %if the database isn't already open, then open it for reading
   if ~obj.pointsToAntelopeDatabase
      obj = obj.open();
   else
      mydb = obj.dbpointer;
   end
   %check to ensure wfdisk table exists and is populated

   if obj.nrecs == 0,
      cleanUpFail('Waveform:load_antelope:databaseNotFound', ...
         'Database not found: %s', dbquery(mydb,'dbTABLE_FILENAME'));
      return;
   end;

   rawDb = mydb;   % keep a copy of the pre-subset (raw) database

   % subset the data based upon the desired criteria
   %criteriaList = keepCriteriaThatMatchDatabaseFields(criteriaList, mydb);
   %allExp = getAsExpressions(criteriaList);

   %subset the database based on this particular criterion
   mydb = dbsubset(mydb,allExp);
   if safe_dbnrecs(mydb) == 0
      cleanUpFail('Waveform:load_antelope:dataNotFound', 'No records found for criteria [%s].', allExp);
      return;
   end;

   filteredDb = mydb;

   [st, ed] = dbgetv(mydb,'time','endtime');
   %% Get the tracebuf object for this starttime, endtime
   % Loop through all times.  Result is tr(1:numel(starttimes) of all tracebuffers.
   for mytimeIDX = 1:numel(antelope_starts)
      someDataExists = any(antelope_starts(mytimeIDX)<= (ed) & antelope_ends(mytimeIDX) >= (st));
      if someDataExists
         %%% Glenn Thompson 2012/02/06: Occasionally the C program trload_css cannot even load the trace data.
         % This error needs to be handled. So adding a try..catch..end around the original instruction.
         trStructs = safe_trload(db, st, et);
         T = trstruct2SeismicTrace(trStructs);
         try
            obj.trpointer{mytimeIDX} = trload_css(mydb, antelope_starts(mytimeIDX), antelope_ends(mytimeIDX)); %#ok<AGROW>
         catch
            cleanUpFail('Waveform:load_antelope:trload_css failed', ...
               'Database not found: %s', dbquery(mydb,'dbTABLE_FILENAME'));
            disp(obj.dbname)
            disp(allExp)
            starttimes = antelope_starts(mytimeIDX);
            endtimes = antelope_ends(mytimeIDX);
            fprintf('%.0f %.0f\n ',starttimes,endtimes);
            disp(['starttimes: ', datestr(epoch2datenum(starttimes))]);
            disp(['endtimes: ', datestr(epoch2datenum(endtimes))]);
            fprintf('trload_css(mydb, starttimes, endtimes))\n');
            return
         end
         trsplice(obj.trpointer{mytimeIDX},20);
      else
         obj.trpointer{mytimeIDX} = -1; %#ok<AGROW>
      end
   end %mytimeIDX
   closeIfAppropriate(mydb, onFinishCloseDB);

   function cleanUpFail(varargin)
      %cleanUpFail   tidy up database and records and display warning
      closeIfAppropriate(mydb, onFinishCloseDB);
      obj.trpointer = { -1 };
      filteredDb = dbinvalid;
      warning(varargin{:});
   end
end


text:
function bnet = mk_asia_bnet(CPD_type, p, arity)
% MK_ASIA_BNET Make the 'Asia' bayes net.
%
% BNET = MK_ASIA_BNET uses the parameters specified on p21 of Cowell et al,
% "Probabilistic networks and expert systems", Springer Verlag 1999.
%
% BNET = MK_ASIA_BNET('cpt', p) uses random parameters drawn from a Dirichlet(p,p,...)
% distribution. If p << 1, this is nearly deterministic; if p >> 1, this is nearly uniform.
%
% BNET = MK_ASIA_BNET('bool') makes each CPT a random boolean function.
%
% BNET = MK_ASIA_BNET('gauss') makes each CPT a random linear Gaussian distribution.
%
% BNET = MK_ASIA_BNET('orig') is the same as MK_ASIA_BNET.
%
% BNET = MK_ASIA_BNET('cpt', p, arity) can specify non-binary nodes.

if nargin == 0, CPD_type = 'orig'; end
if nargin < 3, arity = 2; end

Smoking = 1;
Bronchitis = 2;
LungCancer = 3;
VisitToAsia = 4;
TB = 5;
TBorCancer = 6;
Dys = 7;
Xray = 8;

n = 8;
dag = zeros(n);
dag(Smoking, [Bronchitis LungCancer]) = 1;
dag(Bronchitis, Dys) = 1;
dag(LungCancer, TBorCancer) = 1;
dag(VisitToAsia, TB) = 1;
dag(TB, TBorCancer) = 1;
dag(TBorCancer, [Dys Xray]) = 1;

ns = arity*ones(1,n);
if strcmp(CPD_type, 'gauss')
    dnodes = [];
else
    dnodes = 1:n;
end
bnet = mk_bnet(dag, ns, 'discrete', dnodes);

switch CPD_type
    case 'orig'
        % true is 2, false is 1
        bnet.CPD{VisitToAsia} = tabular_CPD(bnet, VisitToAsia, [0.99   0.01]);
        bnet.CPD{Bronchitis} = tabular_CPD(bnet, Bronchitis, [0.7 0.4   0.3 0.6]);
        % minka: bug fix
        bnet.CPD{Dys} = tabular_CPD(bnet, Dys, [0.9 0.2 0.3 0.1   0.1 0.8 0.7 0.9]);
        bnet.CPD{TBorCancer} = tabular_CPD(bnet, TBorCancer, [1 0 0 0   0 1 1 1]);
        % minka: bug fix
        bnet.CPD{LungCancer} = tabular_CPD(bnet, LungCancer, [0.99 0.9  0.01 0.1]);
        bnet.CPD{Smoking} = tabular_CPD(bnet, Smoking, [0.5 0.5]);
        bnet.CPD{TB} = tabular_CPD(bnet, TB, [0.99 0.95  0.01 0.05]);
        bnet.CPD{Xray} = tabular_CPD(bnet, Xray, [0.95 0.02  0.05 0.98]);
    case 'bool'
        for i=1:n
            bnet.CPD{i} = boolean_CPD(bnet, i, 'rnd');
        end
    case 'gauss'
        for i=1:n
            bnet.CPD{i} = gaussian_CPD(bnet, i, 'cov', 1*eye(ns(i)));
        end
    case 'cpt'
        for i=1:n
            bnet.CPD{i} = tabular_CPD(bnet, i, p);
        end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [x_g,w_g,phi,p_x,p_y] = twod_shapeiso(x_local,r,s,w)
%-----------------------------------------------------------------------
%  twod_shapeiso.m - computes test functions and derivatives on an
%                    element given element coordinates and Gauss points.
%
%                    ! Note: specifically set up for isoparametric
%                    ! elements.  Use `twod_shape' for linear elements.
%
%  Copyright (c) 2001, Jeff Borggaard, Virginia Tech
%  Version: 1.0
%
%  Usage:    [x_g,w_g,phi,p_x,p_y] = twod_shapeiso(x_local,r,w)
%
%  Variables:     x_local
%                        Coordinates of the element nodes
%                 (r,s)
%                        Coordinates of Gauss points in unit triangle
%                 w
%                        Gauss weights associated with (r,s)
%
%                 x_g
%                        Coordinates of Gauss points in the element
%                 w_g
%                        Gauss weights scaled by the element Jacobian
%                 phi
%                        Value of element shape functions at x_g
%                 p_x
%                 p_y
%                        First spatial derivatives of phi
%                 p_xx
%                        Second spatial derivatives of phi
%                        (only applicable with 2nd order or higher)
%                        (currently not implemented)
%-----------------------------------------------------------------------
  x      = x_local;
  [n,t1] = size(x);  % t1 had better be 2
  rule   = length(r);

  if (n == 6)
    % Transform coordinates for quadratic element
    c0 =    x(1,:);
    c1 =-3.*x(1,:) -    x(2,:) + 4.*x(4,:);
    c2 =-3.*x(1,:) -    x(3,:) + 4.*x(6,:);
    c3 = 2.*x(1,:) + 2.*x(2,:) - 4.*x(4,:);
    c4 = 4.*x(1,:) - 4.*x(4,:) + 4.*x(5,:) - 4.*x(6,:);
    c5 = 2.*x(1,:) + 2.*x(3,:) - 4.*x(6,:);

    x_g(:,1) = c0(1) + c1(1)*r + c2(1)*s + c3(1)*r.*r +...
               c4(1)*r.*s + c5(1)*s.*s;
    xr  = c1(1) + 2.*c3(1)*r + c4(1)*s;
    xs  = c2(1) + c4(1)*r + 2.*c5(1)*s;

    x_g(:,2) = c0(2) + c1(2)*r + c2(2)*s + c3(2)*r.*r +...
               c4(2)*r.*s + c5(2)*s.*s;
    yr  = c1(2) + 2.*c3(2)*r + c4(2)*s;
    ys  = c2(2) + c4(2)*r + 2.*c5(2)*s;

    w_g = (xr.*ys-yr.*xs).*w;

    % Compute the Jacobian of the (r,s) -> (x,y) transformation
    jac = 1./( xr.*ys - yr.*xs );
    rx  = ys.*jac;
    sx  =-yr.*jac;
    ry  =-xs.*jac;
    sy  = xr.*jac;

    % Compute shape function and derivatives at Gauss points
    phi = zeros(rule,n);
    phi(:,1) = 1. - 3.*r - 3.*s + 2.*r.*r + 4.*r.*s + 2.*s.*s;
    phi(:,2) =    - 1.*r        + 2.*r.*r                    ;
    phi(:,3) =           - 1.*s                     + 2.*s.*s;
    phi(:,4) =      4.*r        - 4.*r.*r - 4.*r.*s          ;
    phi(:,5) =                              4.*r.*s          ;
    phi(:,6) =             4.*s           - 4.*r.*s - 4.*s.*s;

    p_x = zeros(rule,n);
    p_x(:,1) = ( -3. + 4.*r + 4.*s ).*rx + ( -3. + 4.*r + 4.*s ).*sx;
    p_x(:,2) = ( -1. + 4.*r        ).*rx                            ;
    p_x(:,3) =                             ( -1.        + 4.*s ).*sx;
    p_x(:,4) = (  4. - 8.*r - 4.*s ).*rx + (     - 4.*r        ).*sx;
    p_x(:,5) = (              4.*s ).*rx + (       4.*r        ).*sx;
    p_x(:,6) = (            - 4.*s ).*rx + (  4. - 4.*r - 8.*s ).*sx;

    p_y = zeros(rule,n);
    p_y(:,1) = ( -3. + 4.*r + 4.*s ).*ry + ( -3. + 4.*r + 4.*s ).*sy;
    p_y(:,2) = ( -1. + 4.*r        ).*ry                            ;
    p_y(:,3) =                             ( -1.        + 4.*s ).*sy;
    p_y(:,4) = (  4. - 8.*r - 4.*s ).*ry + (     - 4.*r        ).*sy;
    p_y(:,5) = (              4.*s ).*ry + (       4.*r        ).*sy;
    p_y(:,6) = (            - 4.*s ).*ry + (  4. - 4.*r - 8.*s ).*sy;
  else
    fprintf('element not supported');
  end

author: johannesgerer

repo: jburkardt-m


text:
function B = bernoulli(N)
%BERNOULLI   Bernoulli polynomials as chebfuns.
%   B = BERNOULLI(N) returns a quasimatrix of the first N+1 Bernoulli
%   polynomials on [0,1].
%
%   Example (Bernolli numbers):
%      B = cheb.bernoulli(4);
%      format rat, B(0,:)

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

x = chebfun('x', [0,1]);

% Initalize the constant 1 as the first Bernoulli polynomial.
B = 0*x + 1;
% Compute the requested polynomials.
for j = 1:N
    B(:,j+1) = j*cumsum(B(:,j));
    B(:,j+1) = B(:,j+1) - sum(B(:,j+1));
end

end


text:
function value = cube_monomial ( a, b, expon )

%*****************************************************************************80
%
%% CUBE_MONOMIAL integrates a monomial over a cube in 3D.
%
%  Discussion:
%
%    This routine integrates a monomial of the form
%
%      product ( 1 <= dim <= 3 ) x(dim)^expon(dim)
%
%    The combination 0^0 should be treated as 1.
%
%    The integration region is:
%      A(1) <= X <= B(1)
%      A(2) <= Y <= B(2)
%      A(3) <= Z <= B(3)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(3), B(3), the lower and upper limits.
%
%    Input, integer EXPON(3), the exponents.
%
%    Output, real VALUE, the integral of the monomial.
%
  for i = 1 : 3

    if ( expon(i) == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CUBE_MONOMIAL - Fatal error!\n' );
      fprintf ( 1, '  Exponent of -1 encountered.\n' );
      error ( 'CUBE_MONOMIAL - Fatal error!' );
    end

  end

  value = 1.0;

  for i = 1 : 3

    value = value * ( b(i) ^ ( expon(i) + 1 ) - a(i) ^ ( expon(i) + 1 ) ) ...
      / ( expon(i) + 1 );

  end

  return
end


text:
function varargout = drawDigraph(varargin)
%DRAWDIGRAPH Draw a directed graph, given as a set of vertices and edges
%
%   drawDigraph(NODES1, NODES2, EDGES)
%   NODES1 are originating vertices
%   NODES2 are destination vertices
%   EDGES is an array, with first column containing index of origin vertex
%   (index in NODES1), and second column containing index of destination
%   vertex (index in NODES2).
%   Edges are drawn with arrows.
%
%   H = drawDigraph(...)
%   return handle to the set of edges.
%
%
%   -----
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 17/08/2004.
%

%% Initialisations

% check number of arguments
if nargin < 3
    help drawDigraph;
    return;
end

% initialisations
sn1 = 'bo';     % nodes are red circles
sn2 = 'ro';     % nodes are red circles

%% process input arguments

% First extract the graph structure
n1 = varargin{1};
n2 = varargin{2};
e  = varargin{3};
varargin = varargin(4:length(varargin));

% extract drawing style
if ~isempty(varargin)
    sn1 = varargin{1};
end

if length(varargin)>1
    sn2 = varargin{2};
end

%% main drawing processing

hold on;

nodes = [n1 ; n2];
e(:,2) = e(:,2)+length(n1);

% Draw a 2 dimensional directed graph ----------------------

% Draw 2D Edges ----------------------
%if ~strcmp(se, 'none') & size(e, 1)>0
%    he = plot([n1(e(:,1),1) n2(e(:,2),1)]', [n1(e(:,1),2) n2(e(:,2),2)]', se);
%end
he = drawDirectedEdges(nodes, e);

hn = [];
% Draw 2D nodes ----------------------
if ~strcmp(sn1, 'none')
    hn = plot(n1(:,1), n1(:,2), sn1);
end

 % Draw 2D nodes ----------------------
if ~strcmp(sn2, 'none')
    hn = plot(n2(:,1), n2(:,2), sn2);
end

%% format output arguments

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

if nargout==2
    varargout{1} = hn;
    varargout{2} = he;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [mindcf,pmiss,pfa] = plot_mindcf_point(plot_obj,target_prior,plot_args,legend_string)
% Places the mindcf point for the current system.
% Inputs:
%   target_prior: The effective target prior.
%   plot_args: A cell array of arguments to be passed to 'plot' that control
%     the appearance of the curve. See Matlab's help on 'plot' for information.
%   legend_string: Optional. A string to describe this curve in the legend.

if ischar(plot_args)
    plot_args = {plot_args};
end

[mindcf,pmiss,pfa] = fast_minDCF(plot_obj.tar,plot_obj.non,logit(target_prior),true);
if (pfa < plot_obj.pfa_limits(1)) || (pfa > plot_obj.pfa_limits(2))
    log_warning('pfa of %f is not between %f and %f. The mindcf point will not be plotted.\n',pfa,plot_obj.pfa_limits(1),plot_obj.pfa_limits(2))
elseif (pmiss < plot_obj.pmiss_limits(1)) || (pmiss > plot_obj.pmiss_limits(2))
    log_warning('pmiss of %f is not between %f and %f. The mindcf point will not be plotted.\n',pmiss,plot_obj.pmiss_limits(1),plot_obj.pmiss_limits(2))
else
    figure(plot_obj.fh);
    assert(iscell(plot_args))
    lh = plot(probit(pfa),probit(pmiss),plot_args{:});
    if exist('legend_string','var') && ~isempty(legend_string)
	assert(ischar(legend_string))
	plot_obj.add_legend_entry(lh,legend_string,true);
    end
end

end

author: nesl


text:
function D = driving_function_mono_sdm_kx_ls(kx,xs,f,conf)
%DRIVING_FUNCTION_MONO_SDM_KX_LS driving signal for a line source in SDM in
%the kx-domain
%
%   Usage: D = driving_function_mono_sdm_kx_ps(kx,xs,f,conf)
%
%   Input parameters:
%       kx          - kx dimension [nx1]
%       xs          - position of line source / m [1x3]
%       f           - frequency of the monochromatic source / Hz
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function signal [nx1]
%
%   See also: driving_function_mono_sdm_kx

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 4;
nargmax = 4;
narginchk(nargmin,nargmax);
isargmatrix(kx);
isargxs(xs);
isargpositivescalar(f);
isargstruct(conf);

%% ===== Configuration ==================================================
xref = conf.xref;
c = conf.c;
dimension = conf.dimension;
driving_functions = conf.driving_functions;
x0 = conf.secondary_sources.center;
withev = conf.sdm.withev;

%% ===== Computation ====================================================
% Calculate the driving function in time-frequency domain

% Frequency
omega = 2*pi*f;
% Indexes for evanescent contributions and propagating part of the wave field
idxpr = (( abs(kx) <= (omega/c) ));
idxev = (( abs(kx) > (omega/c) ));
D = zeros(1,length(kx));

if strcmp('2D',dimension)

    % === 2-Dimensional ==================================================

    % Ensure 2D
    xs = xs(1:2);

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 2D line source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('2.5D',dimension)

    % === 2.5-Dimensional ================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 2.5D line source.'],upper(mfilename),driving_functions);
    end

elseif strcmp('3D',dimension)

    % === 3-Dimensional ==================================================

    switch driving_functions
    case 'default'
        % --- SFS Toolbox ------------------------------------------------
        to_be_implemented;
    otherwise
        error(['%s: %s, this type of driving function is not implemented ', ...
            'for a 3D line source.'],upper(mfilename),driving_functions);
    end

else
    error('%s: the dimension %s is unknown.',upper(mfilename),dimension);
end

author: sfstoolbox

repo: sfs-matlab


text:
function varargout = ghostscript(cmd)
%GHOSTSCRIPT  Calls a local GhostScript executable with the input command
%
% Example:
%   [status result] = ghostscript(cmd)
%
% Attempts to locate a ghostscript executable, finally asking the user to
% specify the directory ghostcript was installed into. The resulting path
% is stored for future reference.
%
% Once found, the executable is called with the input command string.
%
% This function requires that you have Ghostscript installed on your
% system. You can download this from: http://www.ghostscript.com
%
% IN:
%   cmd - Command string to be passed into ghostscript.
%
% OUT:
%   status - 0 iff command ran without problem.
%   result - Output from ghostscript.

% Copyright: Oliver Woodford, 2009-2010

% Thanks to Jonas Dorn for the fix for the title of the uigetdir window on
% Mac OS.

% Thanks to Nathan Childress for the fix to the default location on 64-bit
% Windows systems.

% 27/4/11 - Find 64-bit Ghostscript on Windows. Thanks to Paul Durack and
% Shaun Kline for pointing out the issue

% 4/5/11 - Thanks to David Chorlian for pointing out an alternative
% location for gs on linux.

% Call ghostscript
[varargout{1:nargout}] = system(sprintf('"%s" %s', gs_path, cmd));
return

function path = gs_path
% Return a valid path
% Start with the currently set path
path = user_string('ghostscript');
% Check the path works
if check_gs_path(path)
    return
end
% Check whether the binary is on the path
if ispc
    bin = {'gswin32c.exe', 'gswin64c.exe'};
else
    bin = {'gs'};
end
for a = 1:numel(bin)
    path = bin{a};
    if check_store_gs_path(path)
        return
    end
end
% Search the obvious places
if ispc
    default_location = 'C:\Program Files\gs\';
    dir_list = dir(default_location);
    if isempty(dir_list)
        default_location = 'C:\Program Files (x86)\gs\'; % Possible location on 64-bit systems
        dir_list = dir(default_location);
    end
    executable = {'\bin\gswin32c.exe', '\bin\gswin64c.exe'};
    ver_num = 0;
    % If there are multiple versions, use the newest
    for a = 1:numel(dir_list)
        ver_num2 = sscanf(dir_list(a).name, 'gs%g');
        if ~isempty(ver_num2) && ver_num2 > ver_num
            for b = 1:numel(executable)
                path2 = [default_location dir_list(a).name executable{b}];
                if exist(path2, 'file') == 2
                    path = path2;
                    ver_num = ver_num2;
                end
            end
        end
    end
    if check_store_gs_path(path)
        return
    end
else
    bin = {'/usr/bin/gs', '/usr/local/bin/gs'};
    for a = 1:numel(bin)
        path = bin{a};
        if check_store_gs_path(path)
            return
        end
    end
end
% Ask the user to enter the path
while 1
    if strncmp(computer, 'MAC', 3) % Is a Mac
        % Give separate warning as the uigetdir dialogue box doesn't have a
        % title
        uiwait(warndlg('Ghostscript not found. Please locate the program.'))
    end
    base = uigetdir('/', 'Ghostcript not found. Please locate the program.');
    if isequal(base, 0)
        % User hit cancel or closed window
        break;
    end
    base = [base filesep];
    bin_dir = {'', ['bin' filesep], ['lib' filesep]};
    for a = 1:numel(bin_dir)
        for b = 1:numel(bin)
            path = [base bin_dir{a} bin{b}];
            if exist(path, 'file') == 2
                if check_store_gs_path(path)
                    return
                end
            end
        end
    end
end
error('Ghostscript not found. Have you installed it from www.ghostscript.com?');

function good = check_store_gs_path(path)
% Check the path is valid
good = check_gs_path(path);
if ~good
    return
end
% Update the current default path to the path found
if ~user_string('ghostscript', path)
    warning('Path to ghostscript installation could not be saved. Enter it manually in ghostscript.txt.');
    return
end
return

function good = check_gs_path(path)
% Check the path is valid
[good message] = system(sprintf('"%s" -h', path));
good = good == 0;
return

author: opencobra

repo: cobratoolbox


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtRayCube.m                                  |
%|    #    |   VERSION    : 0.41                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 01.04.2018                                    |
%| ( === ) |   SYNOPSIS   : Ray tracing with absorbing cube               |
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Parameters
L     = [5 4 3]
% Xsrc  = [3.9 2.8 1.7]
% Xmes  = [2.1 2.2 1.3]
Xsrc  = [4 2 1.7];
Xmes  = [2 2 1.7];
Nray  = 1e5
rad   = 0.2

% Read mesh
mesh = msh('cube1e1.mesh');

% Mesh reshape to feat dimensions in L
mesh.vtx = 0.5 .* (1 + mesh.vtx);
mesh.vtx = (ones(size(mesh.vtx,1),1)*L) .* mesh.vtx;

% Initialize ray
ray = ray(mesh,Xsrc,Nray);

% Graphical sphere
[X,Y,Z] = sphere(50);
X = Xmes(1) + rad*X; Y = Xmes(2) + rad*Y; Z = Xmes(3) + rad*Z;

% Graphical representation
figure
plot(mesh)
hold on
% plot(ray)
surf(X,Y,Z,ones(size(X)))
axis equal;
xlabel('X');   ylabel('Y');   zlabel('Z');
alpha(0.3)
colorbar

% Material properties (http://www.odeon.dk/material-manufactures)
ray.msh.col(mesh.col==1) = 10;    % X = 1
ray.msh.col(mesh.col==2) = 20;    % Y = 1
ray.msh.col(mesh.col==3) = 30;    % X = 0
ray.msh.col(mesh.col==4) = 40;    % Y = 0
ray.msh.col(mesh.col==5) = 50;    % Z = 0
ray.msh.col(mesh.col==6) = 60;    % Z = 1
% ray.msh.col(:) = 60;

% Maximum distances
r1000 = rad/2 * sqrt(Nray/1000);
r100  = rad/2 * sqrt(Nray/100);
r10   = rad/2 * sqrt(Nray/10);
rMax  = rad/2 * sqrt(Nray/2);

% Ray-tracing
tic
ray = ray.tracer(100,rMax);
toc
% plot(ray)

% Images sources
tic
[img,nrg] = ray.image(Xmes,rad,rMax);
toc

% Graphical representation for images
n   = size(img,1);
tmp = ones(n,1)*Xmes + img;
tmp = msh(tmp,(1:n)');
plot(tmp,10*log10(nrg(:,1)));
axis equal
colorbar

% Analytical solution
load('odeon.mat')
num = [30 ; % X = 0
    10 ;    % X = 1
    40 ;    % Y = 0
    20 ;    % Y = 1
    50 ;    % Z = 0
    60];    % Z = 1
% num = 60 * ones(6,1);
mat = odeon.mat(num,:);
rhm = 30;
air = 5.5 * (50/rhm) .* (odeon.frq/1000).^1.7 * 1e-4;
tic
[imgRef,nrgRef] = rayCubeAnalytic(L,mat,air,Xsrc,Xmes,rMax);
toc

% Graphical representation
plot3(imgRef(:,1)+Xmes(1),imgRef(:,2)+Xmes(2),imgRef(:,3)+Xmes(3),'ok')

% Data
r    = sqrt(sum(img.^2,2));
rRef = sqrt(sum(imgRef.^2,2));
sol  = mean(nrg,2);
ref  = mean(nrgRef,2);

% Energy error in dB
figure
plot(rRef,10*log10(ref),'or',r,10*log10(sol),'+b')
hold on
plot([r1000 r1000],[-60,0],'k--')
text(r1000,-55,' n = 1000')
plot([r100 r100],[-60,0],'k--')
text(r100,-57,' n = 100')
plot([r10 r10],[-60,0],'k--')
text(r10,-55,' n = 10')
plot([rMax rMax],[-60,0],'k--')
text(rMax,-55,' n = 1')
xlabel('Distance source mesure')
ylabel('Energie mesuree (dB)')
title('Energie mesuree selon la distance de mesure')
legend({'Analytique','Statistique'})
grid on

% Audio file
[audio,fs] = audioread('anechoicVoice.wav');

% Fir from images
T = floor(r/340*fs) + 1;
for i = 1:size(nrg,2)
    fir8(:,i) = accumarray(T,nrg(:,i),[max(T) 1]);
end

% Bank filtering
dir8 = ray.bank(256,fs);
rir  = 0;
abs(sum(sum(dir8,2)) - 1)
for i = 1:size(dir8,2)
    rir = rir + fftfilt(dir8(:,i),fir8(:,i));
%     figure
%     freqz(dir8(:,i),1,2*256)
end

% Graphical representation
figure
plot(rir)

% Audio rendering (5 seconds)
out = fftfilt(rir,audio(1:5*fs));
% sound(out,fs)

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan


text:
clear
% define the root name of database
root = '../data_preparation/prepared_data/';

% which scales we're doing
sigma = 1;
num_samples = 5e5;

scales = [0.25,0.35,0.5];
frontalView = 1;

profileViewInds = [2,3,4];

version = 'multi_pie';
ratio_neg = 5;
norm = 1;

data_loc = 'mpie_';
rng(0);

% where to save generated patches for external training
% (e.g. for training CEN patches)
patches_loc = './patches/';
patch_folder = [patches_loc version '/'];

for s=scales
    Save_all_patches(root, frontalView, profileViewInds,...
        s, sigma, version, patch_folder, 'ratio_neg', ratio_neg,...
        'num_samples', num_samples, 'data_loc', data_loc,...
        'normalisation_size', 19);
end

author: TadasBaltrusaitis

repo: OpenFace