text:
classdef Settings %< handle%& matlab.mixin.Copyable

    properties %optmizer access
        optimizerSettings
        plotting = true
        showBC = true
        BCscale_factor = 0.10
        rotation_per_it = 0
        printing = true
        printing_physics = false
        monitoring = true
        monitoring_interval = 10
        maxiter = 2000
        constraint_case = {'EQUALITY'}
        HJiter0
        e2
        ub = 1;
        lb = 0;
        volumeMicro
        superEllipseRatio
    end

    properties %target parameters
        Vfrac_initial
        optimality_initial
        constr_initial
        Vfrac_final
        optimality_final
        constr_final
        epsilon_initial
        epsilon_final
        Perimeter_target
        epsilon_isotropy_initial
        epsilon_isotropy_final
        stressNormExponent_initial
        stressNormExponent_final
    end

    properties
        perimeter = struct;
    end

    properties %LevelSetCreator
        levelSetDataBase
    end

    properties    %topopt access
        ptype
        pdim
        case_file
        filename
        material
        initial_case
        cost
        weights
        constraint
        optimizer
        optimizerUnconstrained
        line_search_initiator
        incrementFactor
        rate
        filter
        unfitted_mesh_algorithm='DELAUNAY'
        TOL = struct;
        target_parameters = struct;
        nsteps
        micro = struct;
        selectiveC_Cstar
        nconstr
        warningHoleBC
        printIncrementalIter
        printChangingFilter
        printMode = 'DesignAndShapes';
        homegenizedVariablesComputer
        materialInterpolation
        designVariable
        vademecumFileName
        nelem
        m1
        m2
        alpha0
        rho0
        isDesignVariableFixed
        costDomainNotOptimizable
        constraintDomainNotOptimizable
    end

    properties %exploring tests
        shFuncParamsName
    end

    methods
        function obj = Settings(case_file)
            run(case_file)
            obj.case_file=case_file;
            obj.filename = filename;
            obj.ptype = ptype;

            if exist('method','var')
                obj.materialInterpolation = method;
            end

            if exist('materialType','var')
                obj.material = materialType;
            end

            if exist('initial_case','var')
                 obj.initial_case = initial_case;
                 if isequal(initial_case,'full')
                    obj.levelSetDataBase.type = 'initial_case';
                 end
            end

            if exist('isDesignVariableFixed','var')
                femReader = FemInputReader_GiD();
                s = femReader.read(obj.filename);
                coord = s.mesh.coord;
                obj.isDesignVariableFixed.nodes  = isDesignVariableFixed.nodes(coord);
                obj.isDesignVariableFixed.values = isDesignVariableFixed.values(coord);
            end

            if exist('costDomainNotOptimizable','var')
                femReader = FemInputReader_GiD();
                s = femReader.read(obj.filename);
                coord = transpose(s.mesh.computeBaricenter);
                for i = 1:numel(costDomainNotOptimizable)
                    cD = costDomainNotOptimizable{i};
                    if ~isempty(cD)
                        obj.costDomainNotOptimizable{i} = cD(coord);
                    else
                        obj.costDomainNotOptimizable{i} = false(size(coord,1));
                    end
                end
            end
            if exist('constraintDomainNotOptimizable','var')
                femReader = FemInputReader_GiD();
                s = femReader.read(obj.filename);
                coord = transpose(s.mesh.computeBaricenter);
                for i = 1:numel(constraintDomainNotOptimizable)
                    cD = constraintDomainNotOptimizable{i};
                    if ~isempty(cD)
                        obj.constraintDomainNotOptimizable{i} = cD(coord);
                    else
                        obj.constraintDomainNotOptimizable{i} = false(size(coord,1),1);
                    end
                end
            end

            obj.cost = cost;
            obj.weights = weights;
            obj.constraint = constraint;
            obj.nconstr = length(constraint);
            obj.optimizer = optimizer;
            if exist('incrementFactor','var')
               obj.incrementFactor = incrementFactor;
            end
            obj.filter = filterType;
            obj.nsteps = nsteps;
            if exist('TOL','var')
                obj.TOL.rho_plus = TOL.rho_plus;
                obj.TOL.rho_minus = TOL.rho_minus;
                obj.TOL.E_plus = TOL.E_plus;
                obj.TOL.E_minus = TOL.E_minus;
                obj.TOL.nu_plus = TOL.nu_plus;
                obj.TOL.nu_minus = TOL.nu_minus;
            end

            obj.Vfrac_initial = Vfrac_initial;
            obj.optimality_initial  = optimality_initial;
            obj.constr_initial = constr_initial;
            obj.optimality_final = optimality_final;
            obj.constr_final = constr_final;

            if exist('line_search_initiator','var')
                obj.line_search_initiator = line_search_initiator;
            else
                if strcmp(obj.optimizer,'PROJECTED GRADIENT')
                    obj.line_search_initiator = 'INCREASING LAST STEP';
                else
                    obj.line_search_initiator = 'STANDARD';
                end
            end

            if exist('constraint_case','var')
                obj.constraint_case = constraint_case;
            end

            if exist('plotting','var')
                obj.plotting = plotting;
            end
            if exist('showBC','var')
                obj.showBC = showBC;
            end
            if exist('BCscale_factor','var')
                obj.BCscale_factor = BCscale_factor;
            end
            if exist('rotation_per_it','var')
                obj.rotation_per_it = rotation_per_it;
            end
            if exist('printing','var')
                obj.printing = printing;
            end
            if exist('printing_physics','var')
                obj.printing_physics = printing_physics;
            end
            if ~obj.printing && obj.printing_physics
                warning('Physical variables will not be printed.')
            end
            if exist('monitoring','var')
                obj.monitoring = monitoring;
            end
            if exist('monitoring_interval','var')
                obj.monitoring_interval = monitoring_interval;
            end

%             if ~contains(case_file,'test','IgnoreCase',true)
%                 fprintf('Loaded %s: \n -Optimizer: %s \n -Cost: ',case_file,obj.optimizer)
%                 fprintf('%s, ',obj.cost{:})
%                 fprintf('\n -Constraint: ')
%                 fprintf('%s, ', obj.constraint{:})
%                 fprintf('\n -Incremental Steps: %f \n ',obj.nsteps)
%             end

            if exist('maxiter','var')
                obj.maxiter = maxiter;
                if ~contains(case_file,'test','IgnoreCase',true)
%                     fprintf('-Max iters: %f \n ',obj.maxiter)
                end
            end

            if exist('Vfrac_final','var')
                obj.Vfrac_final = Vfrac_final;
                if ~contains(case_file,'test','IgnoreCase',true)
%                     fprintf('-Volume target: %f \n ',obj.Vfrac_final)
                end
            end
            if exist('Perimeter_target','var')
                obj.Perimeter_target = Perimeter_target;
                if ~contains(case_file,'test','IgnoreCase',true)
%                     fprintf('-Perimeter target: %f \n',obj.Perimeter_target)
                end
            end
            if exist('epsilon_initial','var')
                obj.epsilon_initial = epsilon_initial;
                obj.epsilon_final = epsilon_final;
            end
            if exist('HJiter0','var')
                obj.HJiter0 = HJiter0;
            end
            if exist('e2','var')
                obj.e2 = e2;
            else
                obj.e2 = 1;
            end

            if exist('micro','var')
                obj.micro.alpha = micro.alpha;
                obj.micro.beta = micro.beta;
            end
            if exist('epsilon_isotropy_initial','var')
                obj.epsilon_isotropy_initial = epsilon_isotropy_initial;
                obj.epsilon_isotropy_final = epsilon_isotropy_final;
            end
            if exist('selectiveC_Cstar','var')
                obj.selectiveC_Cstar = selectiveC_Cstar;
            end

            if exist('widthH','var')
                obj.levelSetDataBase.widthH = widthH;
            end

            if exist('widthV','var')
                obj.levelSetDataBase.widthV = widthV;
            end

            if exist('widthSquare','var')
                obj.levelSetDataBase.widthSquare = widthSquare;
            end

            if exist('N_holes','var')
                obj.levelSetDataBase.nHoles = N_holes;
            end
            if exist('R_holes','var')
                obj.levelSetDataBase.rHoles = R_holes;
            end
            if exist('phase_holes','var')
                obj.levelSetDataBase.phaseHoles = phase_holes;
            end

            if exist('warningHoleBC','var')
                obj.levelSetDataBase.warningHoleBC = warningHoleBC;
            end

            if exist('fracRadius','var')
                obj.levelSetDataBase.fracRadius = fracRadius;
            end

            if exist('levFib','var')
                obj.levelSetDataBase.levFib = levFib;
            end

            if exist('yn','var')
                obj.levelSetDataBase.yn = yn;
            end

            if exist('levelFibers','var')
               obj.levelSetDataBase.levelFibers = levelFibers;
            end

            if exist('volumeFibers','var')
               obj.levelSetDataBase.volumeFibers = volumeFibers;
            end

            if exist('shFuncParamsName','var')
               obj.shFuncParamsName = shFuncParamsName;
            end

            if exist('designVariable','var')
               obj.designVariable = designVariable;
            end

            if exist('homegenizedVariablesComputer','var')
               obj.homegenizedVariablesComputer = homegenizedVariablesComputer;
            else
               obj.homegenizedVariablesComputer = 'ByInterpolation';

            end

            if exist('vademecumFileName','var')
                obj.vademecumFileName = vademecumFileName;
            end

            if  ~(contains(filename,'test','IgnoreCase',true) || contains(filename,'RVE') || obj.hasToAddSpaceBecauseOfIncremental())
                fprintf('\n')
            end

            if ~exist('optimizerUnconstrained','var')
                switch obj.optimizer
                    case {'SLERP','HAMILTON-JACOBI','PROJECTED GRADIENT'}
                        obj.optimizerUnconstrained = obj.optimizer;
                        obj.optimizer = 'AlternatingPrimalDual';
                end
            else
                obj.optimizerUnconstrained = optimizerUnconstrained;
            end

            if exist('rate','var')
                obj.rate = rate;
            else
                obj.rate = 0.5;
            end

            if exist('ub','var')
                obj.ub = ub;
            else
                obj.ub = 1;
            end

            if exist('lb','var')
                obj.lb = lb;
            else
                obj.lb = 0;
            end

            if exist('superEllipseRatio','var')
               obj.levelSetDataBase.vigdergauzDataBase.superEllipseRatio = superEllipseRatio;
            end

            if exist('volumeMicro','var')
               obj.levelSetDataBase.vigdergauzDataBase.volumeMicro = volumeMicro;
            end

            if exist('vigdergauzType','var')
               obj.levelSetDataBase.vigdergauzDataBase.type =  vigdergauzType;
            end

            if exist('vigdergauzStrainMacro','var')
               obj.levelSetDataBase.vigdergauzDataBase.strain =  vigdergauzStrainMacro;
               if exist('TOL','var')
                   obj.levelSetDataBase.vigdergauzDataBase.E1 = TOL.E_plus;
                   obj.levelSetDataBase.vigdergauzDataBase.E0 = TOL.E_minus;
                   obj.levelSetDataBase.vigdergauzDataBase.nu1 = TOL.nu_plus;
                   obj.levelSetDataBase.vigdergauzDataBase.nu0 = TOL.nu_minus;
               end
            end

            if exist('rho0','var')
                obj.rho0 = rho0;
            end

            if exist('m1','var')
               obj.m1 = m1;
            end

            if exist('m2','var')
               obj.m2 = m2;
            end

            if exist('alpha0','var')
               obj.alpha0 = alpha0;
            end

            if exist('stressNormExponent_initial','var')
                obj.stressNormExponent_initial = stressNormExponent_initial;
            end

            if exist('stressNormExponent_final','var')
                obj.stressNormExponent_final = stressNormExponent_final;
            end

        end

    end

    methods (Access = private)

        function itHas = hasToAddSpaceBecauseOfIncremental(obj)
            itHas = true;
            if ~isempty(obj.printIncrementalIter)
                itHas = obj.printIncrementalIter;
            end
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function H = sims_bfgsi(H0,dg,dx)
% H = bfgsi(H0,dg,dx)
% dg is previous change in gradient; dx is previous change in x;
% 6/8/93 version that updates inverse hessian instead of hessian
% itself.
% Copyright by Christopher Sims 1996.  This material may be freely
% reproduced and modified.
if size(dg,2)>1
   dg=dg';
end
if size(dx,2)>1
   dx=dx';
end
Hdg = H0*dg;
dgdx = dg'*dx;
if (abs(dgdx) >1e-12)
   H = H0 + (1+(dg'*Hdg)/dgdx)*(dx*dx')/dgdx - (dx*Hdg'+Hdg*dx')/dgdx;
else

   H=H0;
end
save H.dat H

author: european-central-bank

repo: BEAR-toolbox


text:
function pass = test_lu( pref )
% Test for LU decomposition of a chebfun2.

if ( nargin == 0 )
    pref = chebfunpref;
end

tol = 100*pref.cheb2Prefs.chebfun2eps;

% Decomposition on [-1,1,-1,1]:
x = chebpts(100);
f = chebfun2( @(x,y) cos( x.*y ) );
[L, U] = lu( f );
g = L * U;
pivPos = f.pivotLocations;

% Accurate decomposition:
pass(1) = norm( fevalm(f, x, x) - fevalm(g, x, x) ) < tol;
pass(2) = norm( diag( L( pivPos( :, 2 ) , :) ) - ones( length(f), 1) ) < tol;
% pass(j) = norm( triu( L( pivPos( :, 2 ) , :) , 1) ) < sqrt(tol); j = j + 1;
pass(3) = norm( tril( U( :, pivPos( :, 1 ) ) , -1) ) < sqrt(tol);

% Try the same thing on different domain:
x = chebpts(100, [-2.1, 4.3]);
y = chebpts(100, [-1, 2.7]);

f = chebfun2( @(x,y) cos( x.*y ), [-2.1, 4.3, -1, 2.7]);
[L, U] = lu( f );
g = L * U;
pivPos = f.pivotLocations;

% Accurate decomposition:
pass(4) = norm( fevalm(f, x, y) - fevalm(g, x, y) ) < 2*tol;
pass(5) = norm( diag( L( pivPos( :, 2 ) , :) ) - ones( length(f), 1) ) < tol;
% pass(j) = norm( triu( L( pivPos( :, 2 ) , :) , 1) ) < sqrt(tol); j = j + 1;
pass(6) = norm( tril( U( :, pivPos( :, 1 ) ) , -1) ) < sqrt(tol);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% This is a stub with an example of how to attach voxel data, which will
% trigger the calculation of derivatives (region averages)
% This function will likely not be necessary, as the functions can be accomplished
% with listeners - this is an example really.

% img = fmri_data('/Users/torwager/Downloads/smoothed_residual.nii.gz')
% b = brainpathway();

% Resample image space of data to the atlas space
% (We don't want to do it the other way around, because it will make it
% difficult to have a consistent atlas definition and aggregate across
% different datasets)

img_resampled = resample_space(img, b.region_atlas);

% Assign the resampled data to the brainpathway object
% This triggers several steps:
% Updating node and region averages, and node and region connectivity
b.voxel_dat = img_resampled.dat;

% Plot the inter-region connectivity
% plot_connectivity(b);

% Load a different atlas instead, with a different parcellation:
% This triggers updates of the regions and region connectivity
b.region_atlas = load_atlas('yeo17networks');

b = brainpathway(load_atlas('yeo17networks'));
b.voxel_dat = img_resampled.dat;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [wavespec] = bz_WaveSpec(lfp,varargin)
%[wavespec] = bz_WaveSpec(lfp) calculates the
%wavelet transform of a signal with nfreqs frequencies in the range frange
%[fmin fmax]. Spacing between frequencies can be 'lin' or 'log'.
%Time-frequency resolution is defined by ncyc, the number of cycles in each
%wavelet. Uses Morlet (Gabor) wavelet.
%
%
%INPUT
%    lfp            a buzcode structure with fields lfp.data,
%                                                   lfp.timestamps
%                                                   lfp.samplingRate
%                   -lfp can also be a [t x 1] timeseries signal. in which
%                   case you need to input 'samplingRate'
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%       'frange'    [low frequency, high frequency]     (default: [1 128])
%       'nfreqs'    number of frequencies               (default: 100
%       'roundfreqs' round freqs to unique integer vals (default: false
%                       *Note this may decrease number of freqs
%       'nfreqs'    number of frequencies               (default: 100
%       'ncyc'      number of cycles in your wavelet    (default: 5)
%       'fvector'   predefined vector of frequencies
%       'space'     'log' or 'lin'  spacing of f's      (default: 'log')
%       'samplingRate' (only if input is not a buzcode structure)
%       'intervals'  restrict your spectrogram to timestamps in specific
%                   intervals
%       'chanID'    if lfp structure has multiple channels, which one would
%                   you like to calcaulte the wavelet transform of?
%                   (note: requires field lfp.channels)
%       'showprogress' true/false (default:false)
%       'saveMat '   put the basePath to save an LFP file
%       'MatNameExtraText'   text(X) to add to name as in: 'basename.wavespec(text).lfp.mat'
%       'downsampleout' factor by which to downsample output (default: 1)
%    =========================================================================
%
%OUTPUT
%   wavespec            buzcode-style structure
%       .data           [t x nfreqs] your spectrogram
%       .timestamps     [t x 1] timestamps
%       .freqs          frequencies of each column
%       .samplingRate   (Hz)
%       .channels       Channel indices of channels filtered... taken from lfp input
%       .filterparms    a structure that holds the parameters used for
%                       filtering, for future reference
%

%TO DO:
%   -Matricise For loop?
%   -Don't need to FFT every freq...
%   -update t output for if LFP is a cell array
%
%
%Dependencies
%   WaveFilt
%   MorletWavelet
%   FConv
%
%
%Last Updated: 10/9/15
%DLevenstein
%Modified by Antonio FR, 7/18/18

%% Parse the inputs

%Parameters
parms = inputParser;
addParameter(parms,'frange',[1 128],@isnumeric);
addParameter(parms,'nfreqs',100,@isnumeric);
addParameter(parms,'ncyc',5,@isnumeric);
addParameter(parms,'space','log');
addParameter(parms,'samplingRate',[]);
addParameter(parms,'showprogress',false,@islogical);
addParameter(parms,'roundfreqs',false,@islogical);
addParameter(parms,'saveMatPath',[]);
addParameter(parms,'MatNameExtraText',[]);
addParameter(parms,'fvector',[]);
addParameter(parms,'intervals',[-Inf Inf])
addParameter(parms,'chanID',[])
addParameter(parms,'downsampleout',1)

parse(parms,varargin{:})
frange = parms.Results.frange;
nfreqs = parms.Results.nfreqs;
ncyc = parms.Results.ncyc;
space = parms.Results.space;
samplingRate = parms.Results.samplingRate;
showprogress = parms.Results.showprogress;
roundfreqs = parms.Results.roundfreqs;
saveMatPath = parms.Results.saveMatPath;
MatNameExtraText = parms.Results.MatNameExtraText;
fvector = parms.Results.fvector;
intervals = parms.Results.intervals;
chanID = parms.Results.chanID;
downsampleout = parms.Results.downsampleout;

%Channel restrict
if ~isempty(chanID)
    usechannel = ismember(lfp.channels,chanID);
    lfp.data = lfp.data(:,usechannel);
    lfp.channels = lfp.channels(usechannel);
end

%lfp input
if isstruct(lfp)
    samplingRate = lfp.samplingRate;
elseif isempty(lfp)
    wavespec = lfp;
    return
elseif isnumeric(lfp)
    data_temp = lfp;
    clear lfp
    lfp.data = data_temp;
    lfp.timestamps = [1:length(lfp.data)]'./samplingRate;
end

si = 1./samplingRate;

%Restrict to intervals, with overhang to remove edge effects at transitions
%(then remove later)
overhang = (ncyc)./frange(1);
overint = bsxfun(@(X,Y) X+Y,intervals,overhang.*[-1 1]);
keepIDX = InIntervals(lfp.timestamps,overint);
lfp.data = lfp.data(keepIDX,:);
lfp.timestamps = lfp.timestamps(keepIDX);

%%
if ~isa(lfp.data,'single') || ~isa(lfp.data,'double')
    lfp.data = single(lfp.data);
end

%Frequencies
if ~isempty(fvector)
    freqs = fvector;
else
    fmin = frange(1);
    fmax = frange(2);
    if strcmp(space,'log')
        assert(fmin~=0,'Log-spaced frequencies cannot have min of 0')
        freqs = logspace(log10(fmin),log10(fmax),nfreqs);
    elseif strcmp(space,'lin')
        freqs = linspace(fmin,fmax,nfreqs);
    else
        display('Frequency spacing must be "lin" or "log".')
    end
end

if roundfreqs
    freqs = unique(round(freqs));
end

%Filter with wavelets
nfreqs = size(freqs,2);
nchan = size(lfp.data,2);
ntime = ceil(size(lfp.data,1)./downsampleout);
wavespec.data = nan(ntime,nfreqs,nchan);
wavespec.timestamps = downsample(lfp.timestamps,downsampleout);
for cidx = 1:nchan
    for f_i = 1:nfreqs
        if showprogress
            bz_Counter(f_i,nfreqs,'Wavelet Frequency')
        end
        wavelet = MorletWavelet(freqs(f_i),ncyc,si);
         wavespec.data(:,f_i,cidx) = ...
             downsample(FConv(wavelet',lfp.data(:,cidx)),downsampleout);
    end
end

%% Output in buzcode format
%Remove the overhang from intervals
keepIDX = InIntervals(wavespec.timestamps,intervals);
wavespec.data = wavespec.data(keepIDX,:);
wavespec.timestamps = wavespec.timestamps(keepIDX);

wavespec.freqs = freqs;
wavespec.nfreqs = nfreqs;
wavespec.samplingRate = samplingRate./downsampleout;
if isstruct(lfp) && isfield(lfp,'channels')
    wavespec.channels = lfp.channels;
end
wavespec.filterparms.ncyc = ncyc;
wavespec.filterparms.nfreqs = nfreqs;
wavespec.filterparms.frange = frange;
wavespec.filterparms.space = space;

clear lfp

if saveMatPath
    baseName = bz_BasenameFromBasepath(saveMatPath);
    if ~isempty(MatNameExtraText)
        lfpfilename = fullfile(saveMatPath,[baseName,'.wavespec' MatNameExtraText '.lfp.mat']);
    else
        lfpfilename = fullfile(saveMatPath,[baseName,'.wavespec.lfp.mat']);
    end

    s = whos('wavespec');
    if s.bytes>=1073741824%if greater than 2GB
        disp('wavespec variable greater than 2GB, saving as v7.3 .mat file')
        save(lfpfilename,'wavespec','-v7.3')
    else
        save(lfpfilename,'wavespec')
    end
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
a = Twist(SE3([1 2 3]))
%a.SE
% a.Ad
x = inv(a.SE);
x.Ad

Xtrans([1 2 3])

a = Twist(SE3.Rx(pi/2))
%a.SE
x = inv(a.SE);
x.Ad

Xrotx(pi/2)

author: petercorke

repo: spatialmath-matlab


text:
function [Rsp,U,V]=standardUVBeamPattern(T,xyPoints,normVals,Sigma,lineParams,bounds,numPoints)
%%STANDARDUVBEAMPATTERN Compute the beam pattern of a tapered array as a
%                       function of direction given in terms of direction
%                       cosines for a narrowband linear or planar array.
%                       The beam pattern is just the sum of all of the
%                       tapered elements taken for signals in different
%                       directions. This function can output the response
%                       over all (-1,+1) u-v values or in a rectangular
%                       subset. It can also output a linear cut across the
%                       values.
%
%INPUTS: T The numSubarraysXnumElements tapering matrix of the array. If
%          there are no subarrays, then this can be a numElementsX1
%          vector of weights for every element. The weights can be complex,
%          which means that difference beams can be formed and steering can
%          be taken into account. If an empty matrix is passed, then it is
%          assumed that no tapering is used so T will be a 1XnumElements
%          vector of all ones.
% xyPoints A 1XNumDim or 2XnumDim matrix of the [x;y] locations of the
%          points in the linear or planar array. The units of the distances
%          are in terms of wavelengths for the narrowband model.
% normVals This indicates how the sum beam value should be normalized.
%          Possible values are:
%          'ArrayGain' Return the ratio of the output power to the noise
%                      power. This is the default if this parameter is
%                      omitted or an empty matrix is passed.
%          'NormPowGain' Return the power of the output normalized such
%                        that the highest value is 1.
%          'AbsPowGain' Return the power of the output. This assumes that
%                       the tapering matrix contains the true gain or loss
%                       values for the tapering and is thus not scaled by
%                       any constant value.
%          'NormRealVal' Display the normalized real component of the
%                       output value. It is normalized such that the
%                       largest absolute value is one. Note that this is
%                       not squared, so it corresponds to an output
%                       voltage, not a power.
%           'RawOutput' Provide the raw sum beam output.
%   Sigma If the array gain is desired (normVals='ArrayGain'), then this
%         parameter is used. This is the numElsXnumEls covariance matrix of
%         the noise at the individual elements in the array. If this
%         parameter is omitted or an empty matrix is passed, then the
%         identity matrix will be used.
% lineParams If xyPoints is a 2XnumDims set of points (for a planar array),
%         then if this parameter is provided and is not an empty matrix,
%         Rsp will be a 1D cut across u and v values rather than all u and
%         v values. The equation for the line along with the beam pattern
%         will be evaluated is v=lineParams.intercept+lineParams.slope*u if
%         the value lineParams.vIndep=false or the vIndep component of
%         lineParams is omitted. Otherwise, the line is
%         u=lineParams.intercept+lineParams.slope*v .
%  bounds A 2X1 (or 1X2) or a 4X1 (or 1X4) vector with the bounds in u and
%         v of the plot. For 1-dimensional plots, which are the case if
%         xyPoints is 1XnumDim or lineParams is provided and xyPoints is
%         2XnumDim, then bounds=[minVal;maxVal] for the independent
%         variable. For two-dimensional plots, then
%         bounds=[minU;maxU;minV;maxV]. If this parameter is omitted or an
%         empty matrix is passed, then [-1;1;-1;1] is used to go over all
%         possible values in u and v (or just u for a 1D plot).
% numPoints A parameter determining the size of the output matrix. For 1D
%         plots, Rsp is numPointsX1. For 2D plots, Rsp is
%         numPointsXnumPoints. If this parameter is omitted or an empty
%         matrix is passed, then a default of numPoints=125 points is used.
%
%OUTPUTS: Rsp The array beam pattern over the selected region, normalized
%             as specified. The value in entry Rsp(i,j) corresponds to the
%             U and V values U(i,j) and V(i,j). For a 1D response for a
%             linear array, Rsp(i) corresponds to U(i) and if xyPoints is
%             2XnumDims, V(i) is the corresponding v value. The value of
%             Rsp for points outside of the visible region (u^2+v^2>1)
%             is set to zero.
%           U A matrix of points corresponding to the u values of the
%             responses in Rsp.
%           V A matrix of points corresponding to the v values of the
%             responses in Rsp. For 1D responses,
%
%The idea of a beam pattern for a narrowband array is discussed in Chapter
%2.2 of [1]. The array gain is discussed in Chapter 2.6.2.
%
%EXAMPLE 1:
%Here, we find the array gain beam pattern of a 20 element 1D linear array
%with Taylor tapering and half-wavelength element spacing. When consdiering
%the array gain, it does not matter that the elements are not provided
%centered about the origin. However, the points must be centered to get the
%proper tapering values from the TaylorLinearTapering function.
% xPoints=0:0.5:9.5;
% xPoints=xPoints-mean(xPoints);
% nBar=3;
% sidelobedB=-25;
% T=TaylorLinearTapering(nBar,sidelobedB,xPoints).';
% [Rsp,U]=standardUVBeamPattern(T,xPoints);
% figure(1)
% clf
% plot(U,10*log10(Rsp),'linewidth',2)
% axis([-1 1 -30 15])
% axis square
% h1=xlabel('u');
% h2=ylabel('Array Gain');
% title('Array Power Gain in Decibels')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%EXAMPLE 2:
%In this example, we plot the normalized power gain of a 2D hexagonal array
%without any tapering.
% xyPoints=getShaped2DLattice([7;14],'hexagonal');
% [Rsp,U,V]=standardUVBeamPattern([],xyPoints,'NormPowGain');
% figure(2)
% clf
% surface(U,V,10*log10(Rsp),'EdgeColor','None')
% axis([-1 1 -1 1 -40 0])
% caxis([-40 0])
% colormap(jet(256))
% colorbar()
% view(45,45)
% light()
% h1=xlabel('u');
% h2=ylabel('v');
% h3=zlabel('Array Gain');
% title('Array Power Gain in Decibels')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h3,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%EXAMPLES 3:
%In this example, we plot the difference pattern of a circular array.
%First, we plot the normalized real value of the pattern (the pattern is
%real anyway so this is not an issue) in 2D and then we take 1D cuts of it
%in two different directions. Bayliss tapering are used to obtain the
%difference pattern.
% xyPoints=getShaped2DLattice([30;30],'circular');
% sidelobedB=-30;
% N=17;
% T=BaylissTapering(sidelobedB,N,xyPoints).';
% [Rsp,U,V]=standardUVBeamPattern(T,xyPoints,'NormRealVal');
% figure(3)
% clf
% surface(U,V,Rsp,'EdgeColor','None')
% axis([-1 1 -1 1 -1 1])
% caxis([-1 1])
% colormap(jet(256))
% colorbar()
% view(45,10)
% light()
% h1=xlabel('u');
% h2=ylabel('v');
% h3=zlabel('Real Response');
% title('Array Difference  Beam Pattern')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h3,'FontSize',14,'FontWeight','bold','FontName','Times')
%
% %Now, we take a cut along to v=0 line.
% lineParams=[];
% lineParams.intercept=0;
% lineParams.slope=0;
% lineParams.vIndep=false;
%
% [Rsp,U]=standardUVBeamPattern(T,xyPoints,'NormRealVal',[],lineParams);
% figure(4)
% clf
% plot(U,Rsp,'linewidth',2)
% h1=xlabel('u');
% h2=ylabel('Real Response');
% title('Array Difference Beam Pattern')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
% %And here we take a cut along the u=0 line.
% lineParams.vIndep=true;
%
% [Rsp,~,V]=standardUVBeamPattern(T,xyPoints,'NormRealVal',[],lineParams);
% figure(5)
% clf
% plot(V,Rsp,'linewidth',2)
% h1=xlabel('v');
% h2=ylabel('Real Response');
% title('Array Difference Beam Pattern')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%EXAMPLE 4:
%This is an example of using elements of a circular array with a fixed
%Taylor tapering that have been broken into subarrays. We then try to form
%the best approximation to a Bayliss difference pattern modifying only the
%subarray outputs. We also demonstrate how steering can be used to move the
%difference beam away from the center of the array.
% xyPoints=getShaped2DLattice([30;30],'circular');
% numEls=size(xyPoints,2);
% %Fixed element-level tapering.
% sidelobedB=-30;
% nBar=4;
% g=TaylorTapering(nBar,sidelobedB,xyPoints);
% %It is assumed the disjoint subarrays can be formed (no adjacency matrix
% %used).
% numLevels=3;
% N=17;
% T=double(findBaylissSubarrays(xyPoints,numLevels,sidelobedB,N));
% %We now apply the tapering to the matrix
% T=bsxfun(@times,T,g.');
%
% %The desired element-level tapering
% g=BaylissTapering(sidelobedB,N,xyPoints);
%
% %Next, we try to find the best subarray weights to approximate the desired
% %tapering. We specifically add a null at the boresight.
% g=findSubarrayWeights(T,g,ones(numEls,1));
%
% %Apply the tapering to the elements
% T=bsxfun(@times,g,T);
%
% %We also steer the array off boresight to u0. This is element-level
% %steering.
% u0=[-0.5;0.5];
% D=diag(exp(1j*2*pi*sum(bsxfun(@times,xyPoints,u0),1)));
% T=T*D;
%
% [Rsp,U,V]=standardUVBeamPattern(T,xyPoints,'NormRealVal');
% figure(3)
% clf
% surface(U,V,Rsp,'EdgeColor','None')
% axis([-1 1 -1 1 -1 1])
% caxis([-1 1])
% colormap(jet(256))
% colorbar()
% view(45,10)
% light()
% h1=xlabel('u');
% h2=ylabel('v');
% h3=zlabel('Real Response');
% title('Array Difference  Beam Pattern')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h3,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] H. L. Van Trees, Optimum Array Processing. New York: Wiley-
%    Interscience, 2002.
%
%August 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(xyPoints,1);
numEls=size(xyPoints,2);

if(isempty(T))
   T=ones(1,numEls);
end

if(nargin<7||isempty(numPoints))
   numPoints=125;
end

if(nargin<6||isempty(bounds))
    if(numDim==1)
        bounds=[-1;1];
    else
        bounds=[-1;1;-1;1];
    end
end

if(nargin<5||isempty(lineParams))
    if(numDim>1)
        %If a 2D plot is desired.
        uVals=linspace(bounds(1),bounds(2),numPoints);
        vVals=linspace(bounds(3),bounds(4),numPoints);
        [U,V]=meshgrid(uVals,vVals);
        Rsp=zeros(numPoints,numPoints);

        if(isvector(T))
            for curVal=1:(numPoints*numPoints)
                u=[U(curVal);V(curVal)];
                if(u'*u>1)
                    continue;
                end

                Rsp(curVal)=T*exp(-1j*2*pi*sum(bsxfun(@times,xyPoints,u),1)).';
            end
        else
            for curVal=1:(numPoints*numPoints)
                u=[U(curVal);V(curVal)];
                if(u'*u>1)
                    continue;
                end

                Rsp(curVal)=sum(sum(T*exp(-1j*2*pi*sum(bsxfun(@times,xyPoints,u),1)).'));
            end
        end
    else
        %If a 1D plot is required.
        U=linspace(bounds(1),bounds(2),numPoints)';
        V=[];

        if(isvector(T))
            Rsp=zeros(numPoints,1);
            for curVal=1:numPoints
                if(abs(U(curVal))>1)
                    continue;
                end

                Rsp(curVal)=T*exp(-1j*2*pi*xyPoints*U(curVal)).';
            end
        else
            Rsp=zeros(numPoints,1);
            for curVal=1:numPoints
                if(abs(U(curVal))>1)
                    continue;
                end

                Rsp(curVal)=sum(T*exp(-1j*2*pi*xyPoints*U(curVal)).');
            end
        end
    end
else
    %If a 1D cut across the 2D surface is desired. The
    slope=lineParams.slope;
    intercept=lineParams.intercept;
    if(isfield(lineParams,'vIndep'))
        vIndep=lineParams.vIndep;
    else
        vIndep=false;
    end

    if(vIndep==false)
        %The independent variable is u.
        U=linspace(bounds(1),bounds(2),numPoints)';
        V=zeros(numPoints,1);
        Rsp=zeros(numPoints,1);

        if(isvector(T))
            for curVal=1:numPoints
                V(curVal)=intercept+slope*U(curVal);
                uVec=[U(curVal);V(curVal)];

                if(uVec'*uVec>1)
                    continue;
                end

                Rsp(curVal)=T*exp(-1j*2*pi*sum(bsxfun(@times,xyPoints,uVec),1)).';
            end
        else
            for curVal=1:numPoints
                V(curVal)=intercept+slope*U(curVal);
                uVec=[U(curVal);V(curVal)];

                if(uVec'*uVec>1)
                    continue;
                end

                Rsp(curVal)=sum(sum(T*exp(-1j*2*pi*sum(bsxfun(@times,xyPoints,uVec),1)).'));
            end
        end
    else
        %The independent variable is v.
        V=linspace(bounds(1),bounds(2),numPoints)';
        U=zeros(numPoints,1);
        Rsp=zeros(numPoints,1);

        if(isvector(T))
             for curVal=1:numPoints
                U(curVal)=intercept+slope*V(curVal);
                uVec=[U(curVal);V(curVal)];

                if(uVec'*uVec>1)
                    continue;
                end

                Rsp(curVal)=T*exp(-1j*2*pi*sum(bsxfun(@times,xyPoints,uVec),1)).';
            end
        else
            for curVal=1:numPoints
                U(curVal)=intercept+slope*V(curVal);
                uVec=[U(curVal);V(curVal)];

                if(uVec'*uVec>1)
                    continue;
                end

                Rsp(curVal)=sum(sum(T*exp(-1j*2*pi*sum(bsxfun(@times,xyPoints,uVec),1)).'));
            end
        end
    end
end

if(nargin<3||isempty(normVals))
    normVals='ArrayGain';
end

switch(normVals)
    case 'ArrayGain'%Array power gain versus spatially white noise.
        if(nargin<4||isempty(Sigma))
           Sigma=eye(numEls,numEls);
        end

        if(isvector(T))
            Rsp=abs(Rsp).^2/sum(T*Sigma*T');
        else

            numSubarrays=size(T,1);
            e=ones(numSubarrays,1);
            %The sum sums up all of the subarray outputs.
            Rsp=abs(Rsp).^2/sum(e'*T*Sigma*T'*e);
        end
    case 'NormPowGain'%Normalized absolute array power gain
        Rsp=abs(Rsp);
        Rsp=(Rsp/max(Rsp(:))).^2;%Make it relative to the peak.
    case 'AbsPowGain'
        Rsp=abs(Rsp)^2;
    case 'NormRealVal'%Display the real component of the response.
        Rsp=real(Rsp);
        Rsp=Rsp/max(Rsp(:));%Make it relative to the peak.
    case 'RawOutput'
    otherwise
        error('Unknown plot type requested')
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:
function [abpl, ubpl] = tapas_ti_bpl(llh, t)
%% Computes the Bayesiand predictive log likelihood.
%
%   Input
%       llh     Multidimensional array of preditive likelihood. Dimensions are
%               Number of subjects x number of chains x number of samples
%       t       Temperature schedule of size 1 x number of chains
%   Output
%       abpl    Adjusted Bayesian preditive likelihood
%       ubpl    Unadjusted Bayesian predictive likelihood
%

% aponteeduardo@gmail.com
% copyright (C) 2017
%

[ns, nc, np] = size(llh);
%if all([1, nc] == size(t))
%    error('tapas:ti:bpl', 'Dimensions of t and llh do not match');
%end

% Compute the free energy

lnz = trapz(t, sum(mean(llh, 3), 1));

% Make sure that the first dimension is the temperature
ubpl = trapz(t', mean(llh, 3)')';

% Do post hoc importance sampling
abpl = zeros(ns, 1);
for s = 1:ns
    sllh = bsxfun(@times, -t', squeeze(llh(s, :, :)));
    sllh = bsxfun(@minus, sllh, mean(sllh, 2));
    sllh = bsxfun(@minus, sllh, log(sum(exp(sllh), 2)));
    % Change the dynamic range to max being 100
    adja = -max(sllh')' + log(1000);
    sllh = bsxfun(@plus, sllh, adja);
    sllh = exp(sllh) ;
    nllh = squeeze(sum(skip_column(llh, s), 1));
    nllh = sum(bsxfun(@times, nllh, sllh), 2) .* exp(-adja);
    abpl(s) = lnz - trapz(t, nllh);
end

end

function [nllh] = skip_column(llh, sbj)
% Removes an slice and return array minux slice

[ns, nc, np] = size(llh);
i = 1:ns;
i(sbj) = [];

nllh = llh(i, :, :);

end

author: translationalneuromodeling

repo: tapas


text:
function [prob]=mhprob(val,cand,term,n,N)

% obtain the first exponential term
exp1=exp(-0.5*(exp(-cand)-exp(-val))*term);
% obtain the second exponential term
exp2=exp(-0.5*N*n*(cand-val));

% obtain the acceptance probability
prob=min(1,exp1*exp2);


text:
function [compression_ratio] = flacwrite(varargin);
%flacWRITE Write FLAC (".FLAC") sound file. (Free Lossless Audio Codec)
%    FLACWRITE(Y,FS,NBITS,FLACFILE,ENCODING) writes data Y to a FLAC Vorbis
%    file specified by the file name FLACFILE, with a sample rate
%    of FS Hz and with NBITS number of bits. Stereo data should
%    be specified as a matrix with two columns.
%    ENCODING must be specified as an integer number from 1 to 9
%
%   1 = Fast Encoding.
%   .
%   .
%   .
%   9 = Best Encoding.
%
%   Y,FS and NBITS are mandatory fields. If FLACFILE is not defined the file
%   name will be 'Default_name.flac'. If ENCODING is not defined encoding
%   type '2' will be used by deault.
%
%   [compression_ratio] = flacwrite(Y,FS,NBITS,FLACFILE,ENCODING); returns
%   the compression ratio obtained.
%
%    See also FLACREAD, WAVREAD, WAVWRITE.
if length(varargin) < 3 | length(varargin) > 5
    error('Unsopported number of argument inputs')
end
Y = varargin{1};
FS = varargin{2};
NBITS = varargin{3};
if NBITS~=8 & NBITS~=16 & NBITS~=24 & NBITS~=32
    error('Unsopported bit depth')
end
if length(varargin) >= 4
    FLACFILE = varargin{4};
    if ischar(FLACFILE) ~= 1
        error('File name is not a string')
    end
else
    FLACFILE = 'Default_name.flac';
    disp('File name = Default_name.flac')
end

if isempty(findstr(FLACFILE,'.flac'))
    FLACFILE = strcat(FLACFILE,'.flac');
end

if length(varargin) == 5
    ENCODING = varargin{5};
else
    ENCODING = '2';
    disp('Variable bit rate, joint-stereo, 128 kb/s encoding')
end
app_dir = which('flacwrite.m');
app_dir_base = findstr('flacwrite.m',app_dir);
encoder_dir = app_dir(1:app_dir_base-2);
%%%%% Temporary File %%%%%
[data,channels] = size(Y);
wavwrite(Y,FS,NBITS,strcat(encoder_dir,'\temp.wav'));
tmpfile = strcat(encoder_dir,'\temp.wav');
FLACFILE = strcat(pwd,'\',FLACFILE);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% Data Encoding  using "flac.exe"  %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ENCODING =  num2str(ENCODING);
switch ENCODING
    case {'1'}
        [stat,info_out] = dos([encoder_dir,'\flac.exe', ' ', tmpfile, ' ', ' -0 ', ' -f ',' -o ',FLACFILE]);
    case {'2'}
        [stat,info_out] = dos([encoder_dir,'\flac.exe', ' ', tmpfile, ' ', ' -1 ', ' -f ',' -o ',FLACFILE]);
    case {'3'}
        [stat,info_out] = dos([encoder_dir,'\flac.exe', ' ', tmpfile, ' ', ' -2 ', ' -f ',' -o ',FLACFILE]);
    case {'4'}
        [stat,info_out] = dos([encoder_dir,'\flac.exe', ' ', tmpfile, ' ', ' -3 ', ' -f ',' -o ',FLACFILE]);
    case {'5'}
        [stat,info_out] = dos([encoder_dir,'\flac.exe', ' ', tmpfile, ' ', ' -4 ', ' -f ',' -o ',FLACFILE]);
    case {'6'}
        [stat,info_out] = dos([encoder_dir,'\flac.exe', ' ', tmpfile, ' ', ' -5 ', ' -f ',' -o ',FLACFILE]);
    case {'7'}
        [stat,info_out] = dos([encoder_dir,'\flac.exe', ' ', tmpfile, ' ', ' -6 ', ' -f ',' -o ',FLACFILE]);
    case {'8'}
        [stat,info_out] = dos([encoder_dir,'\flac.exe', ' ', tmpfile, ' ', ' -7 ', ' -f ',' -o ',FLACFILE]);
    case {'9'}
        [stat,info_out] = dos([encoder_dir,'\flac.exe', ' ', tmpfile, ' ', ' -8 ', ' -f ',' -o ',FLACFILE]);
    otherwise
        error('Encoding parameters not suported')
        stat = 2;
end
compression_ratio_ini = findstr(info_out,'bytes, ratio=');
compression_ratio_end = compression_ratio_ini+18;
compression_ratio = info_out(compression_ratio_ini+6:compression_ratio_end);
if stat == 1
    delete(FLACFILE);
    delete(tmpfile);
    error('Error while encoding')
end
%Delete temporary file
delete(tmpfile);

author: Sable

repo: mcbench-benchmarks


text:
%% design path from points by spline

point = [0, 0;
    1, 0;
    2,0
    3,0.5
    4,1.5
    4.8, 1.5
    5,0.8
    6, 0.5
    6.5, 0
    7.5, 0.5
    7,2
    6, 3
    5, 4
    4., 2.5
    3, 3
    2., 3.5;
    1.3, 2.2
    0.5, 2.
    0,3];

s = 1:1:length(point);

px_spline = spline(s, point(:,1), 1:0.01:length(point));
py_spline = spline(s, point(:,2), 1:0.01:length(point));

p_spline = [px_spline', py_spline'];

%% insert yaw

yaw = zeros(length(px_spline), 1);
for i = 2:length(px_spline)-1
    x_forward = px_spline(i+1);
    x_backward = px_spline(i-1);
    y_forward = py_spline(i+1);
    y_backward = py_spline(i-1);
    yaw(i) = atan2(y_forward-y_backward, x_forward-x_backward);
end
yaw(1) = yaw(2);
yaw(end) = yaw(end-1);

%% plot with attitude

arrow_scale = 0.01;

figure(101);
plot(point(:,1), point(:,2), 'bo-'); hold on;
quiver(px_spline', py_spline', cos(yaw)*arrow_scale, sin(yaw)*arrow_scale);
plot(px_spline, py_spline,'r-'); grid on; hold off;

%% save path
path = [px_spline', py_spline', yaw];

save('path', 'path')

author: TakaHoribe


text:
function [ actions ] = BuildActionList
%BuildActionList

%actions for the mountain car problem
actions = [-1.0; 1.0];

author: epfl-lasa


text:
%{
    Demonstration of computing the kernel function matrix.
%}

clc
clear all
close all
addpath(genpath(pwd))

% data
X = rand(33, 2);
Y = rand(23, 2);

% kernel setting
kernel1 = Kernel('type', 'gaussian', 'gamma', 0.01);
kernel2 = Kernel('type', 'polynomial', 'degree', 3);
kernel3 = Kernel('type', 'linear');
kernel4 = Kernel('type', 'sigmoid', 'gamma', 0.01);
kernel5 = Kernel('type', 'laplacian', 'gamma', 0.1);

% compute the kernel function matrix
matrix{1} = kernel1.computeMatrix(X, Y);
matrix{2} = kernel2.computeMatrix(X, Y);
matrix{3} = kernel3.computeMatrix(X, Y);
matrix{4} = kernel4.computeMatrix(X, Y);
matrix{5} = kernel5.computeMatrix(X, Y);

author: iqiukp

repo: KPCA-MATLAB

save_path: github-repos/MATLAB/iqiukp-KPCA-MATLAB


text:
function [ y, m, d, ierror ] = ymd_check_eg_civil ( y, m, d )

%*****************************************************************************80
%
%% YMD_CHECK_EG_CIVIL checks an Egyptian Civil YMD date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, the YMD date, which may
%    be corrected if necessary and possible.
%
%    Output, integer  IERROR, is 0 if the date is legal.
%

%
%  Check the year.
%
  if ( y <= 0 )
    ierror = 1;
    return
  end
%
%  Check the month.
%
  [ y, m, ierror ] = ym_check_eg_civil ( y, m );

  if ( ierror ~= 0 )
    return
  end
%
%  Check the day.
%
  [ y, m, d ] = day_borrow_eg_civil ( y, m, d );

  [ y, m, d ] = day_carry_eg_civil ( y, m, d );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [edge, isInside] = clipRay(ray, bb)
%CLIPRAY Clip a ray with a box.
%
%   EDGE = clipRay(RAY, BOX);
%   RAY is a straight ray given as a 4 element row vector: [x0 y0 dx dy],
%   with (x0 y0) being the origin of the ray and (dx dy) its direction
%   vector, BOX is the clipping box, given by its extreme coordinates:
%   [xmin xmax ymin ymax].
%   The result is given as an edge, defined by the coordinates of its 2
%   extreme points: [x1 y1 x2 y2].
%   If the ray does not intersect the box, [NaN NaN NaN NaN] is returned.
%
%   Function works also if RAY is a N-by-4 array, if BOX is a Nx4 array, or
%   if both RAY and BOX are N-by-4 arrays. In these cases, EDGE is a N-by-4
%   array.
%
%   See also
%     rays2d, boxes2d, edges2d, clipLine, drawRay
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2010-05-13, using Matlab 7.4.0.287 (R2007a)
% Copyright 2010-2022 INRA - Cepia Software Platform

% adjust size of two input arguments
if size(ray, 1) == 1
    ray = repmat(ray, size(bb, 1), 1);
elseif size(bb, 1) == 1
    bb = repmat(bb, size(ray, 1), 1);
elseif size(ray, 1) ~= size(bb, 1)
    error('bad sizes for input');
end

% first compute clipping of supporting line
edge = clipLine(ray, bb);

% detects valid edges (edges outside box are all NaN)
inds = find(isfinite(edge(:, 1)));

% compute position of edge extremities relative to the ray
pos1 = linePosition(edge(inds,1:2), ray(inds,:), 'diag');
pos2 = linePosition(edge(inds,3:4), ray(inds,:), 'diag');

% if first point is before ray origin, replace by origin
edge(inds(pos1 < 0), 1:2) = ray(inds(pos1 < 0), 1:2);

% if last point of edge is before origin, set all edge to NaN
edge(inds(pos2 < 0), :) = NaN;

% eventually returns result about inside or outside
if nargout > 1
    isInside = isfinite(edge(:,1));
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function FBL = ffblasius(RE)
% FFBLASIUS Blasius skin friction factor
%  FFCW(RE) returns the skin friction loss coefficient
%  using the Blasius correlation.
%  CALLED FUNCTION: none
%   Required Inputs are:
%    RE  - Reynolds number (-)
% ---------------------------------------------------------------
% The MATLAB function was created by Tibor Balint, December 1998
% TBoreal Research Corporation, Toronto, Ont. Canada
% (tibor@netcom.ca) and also, University of Warwick, UK
% ---------------------------------------------------------------

format long g;                  % set the format of the calculations

if or((RE<4000),(RE>100000))
   error('The valid Reynolds number range is (4000 < RE < 100000)')
end

%set the correlation constants
a=0.3164;
b=-0.25;

%calculate the Blasius skin friction factor
FBL=a*RE^b;

return                          %end of the function
% -------------- end of the function ----------------------


text:
classdef TestBitwiseXor
    %TestBitwiseXor

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','img001.jpg');
    end

    methods (Static)
        function test_rgb_images
            img1 = cv.imread(TestBitwiseXor.im, 'ReduceScale',2);
            img2 = cv.cvtColor(img1, 'RGB2HSV');

            % rectangular mask
            [h,w,~] = size(img1);
            mask = false([h w]);
            mask(100:h-100,100:w-100) = true;

            out = cv.bitwise_xor(img1, img2);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_xor(img1, img2);
            assert(isequal(out, expected));

            out = cv.bitwise_xor(img1, img2, 'Mask',mask);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_xor(img1, img2, mask);
            assert(isequal(out, expected));

            out = cv.bitwise_xor(img1, img2, 'Mask',mask, 'Dest',img1);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_xor(img1, img2, mask, img1);
            assert(isequal(out, expected));
        end

        function test_float_images
            img1 = cv.imread(TestBitwiseXor.im, 'ReduceScale',2);
            img1 = single(img1) ./ 255;
            img2 = cv.cvtColor(img1, 'RGB2HSV');
            img2(:,:,1) = img2(:,:,1) ./ 360;

            % circular mask
            [h,w,~] = size(img1);
            [X,Y] = meshgrid(1:w,1:h);
            c = fix([w h]/2); r = 50;
            mask = ((X-c(1)).^2 + (Y-c(2)).^2) < r^2;

            out = cv.bitwise_xor(img1, img2);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_xor(img1, img2);
            assert(isequaln(out, expected));

            out = cv.bitwise_xor(img1, img2, 'Mask',mask);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_xor(img1, img2, mask);
            assert(isequaln(out, expected));

            out = cv.bitwise_xor(img1, img2, 'Mask',mask, 'Dest',img1);
            validateattributes(out, {class(img1)}, {'size',size(img1)});
            expected = my_bitwise_xor(img1, img2, mask, img1);
            assert(isequaln(out, expected));
        end

        function test_vectors
            A = uint8([240 255 85]);
            B = uint8([15 15 170]);
            mask = [true false true];

            % uint8([255 240 255])
            out = cv.bitwise_xor(A, B);
            validateattributes(out, {class(A)}, {'size',size(A)});
            expected = bitxor(A, B);
            assert(isequal(out, expected));

            % uint8([255 0 255])
            out = cv.bitwise_xor(A, B, 'Mask',mask);
            validateattributes(out, {class(A)}, {'size',size(A)});
            expected = bitxor(A, B);
            expected(~mask) = 0;
            assert(isequal(out, expected));

            % uint8([255 255 255])
            out = cv.bitwise_xor(A, B, 'Mask',mask, 'Dest',A);
            validateattributes(out, {class(A)}, {'size',size(A)});
            expected = bitxor(A, B);
            expected(~mask) = A(~mask);
            assert(isequal(out, expected));
        end

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

end

function out = my_bitwise_xor(src1, src2, mask, dst)
    %MY_BITWISE_XOR  Similar to cv.bitwise_xor using core MATLAB functions

    if nargin < 3, mask = true(size(src1,1), size(src1,2)); end
    if nargin < 4, dst = zeros(size(src1), class(src1)); end

    % calculate bitwise XOR
    if isinteger(src1)
        out = bitxor(src1, src2);
    elseif isfloat(src1)
        if isa(src1, 'double')
            klass = 'uint64';
        elseif isa(src1, 'single')
            klass = 'uint32';
        end
        out = bitxor(typecast(src1(:), klass), typecast(src2(:), klass));
        out = reshape(typecast(out, class(src1)), size(src1));
    end

    % apply masking
    for k=1:size(out,3)
        out_slice = out(:,:,k);
        dst_slice = dst(:,:,k);
        out_slice(~mask) = dst_slice(~mask);
        out(:,:,k) = out_slice;
    end
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function sol = gsp_regression_tik(G ,M, y , tau, param )
%GSP_REGRESSION_TIK Regression using graph and Tikhonov
%   Usage: sol = gsp_regression_tik(G ,M, y , tau );
%          sol = gsp_regression_tik(G ,M, y , tau, param );
%
%   Input parameters:
%       G   : Graph
%       M   : Mask (to determine with label is known)
%       y   : label (total size of the problem)
%       tau : regularization parameter (weight for tv)
%       param : optional structure of parameters
%
%   Output parameters:
%       sol : Solution of the problem
%
%   This function solve the following problem
%
%   .. argmin_x  || M x - y ||_2^2 + tau || nabla_G x ||_2^2
%
%   If tau is set to zero, then the following problem is solved
%
%   ..  argmin_x   || nabla_G x ||_2^2    s. t.  M x - y = 0
%
%   For the las problem, this function can compute an exact solution if
%   *param.exact* is activated. It will be efficient if the number of
%   unlabelled points is low.
%
%   Additional parameters
%   ---------------------
%
%   * *param.verbose* : Verbosity of the iterative algorithm
%   * *param.maxit* : maximum number of iteration for PCG
%
%   This function uses the UNLocBoX.
%
%   See also: gsp_classification_tik gsp_regression_tv
%

% Author: Nathanael Perraudin
% Date  : 24 July 2015
% Testing: test_graph_ml

%% Optional parameters

if nargin<5
    param = struct;
end

if nargin<4
    tau = 0;
end

if ~isfield(param,'verbose'), param.verbose = 1; end
if ~isfield(param,'tol'), param.tol = 1e-6; end
if ~isfield(param,'maxit'), param.maxit = 200; end
% if ~isfield(param,'direct'), param.direct = (tau==0); end
% if ~isfield(param,'exact'), param.exact = (numel(M)-nnz(M))<1000; end

if tau==0
    if param.verbose
        fprintf('Using direct solution \n')
    end
    if (numel(M) == size(M,1)) || (numel(M) == size(M,2))
        indl = find(M);
        indu = find(1-M);
    else
        error('I cannot handle this case yet');
    end
    Luu = (G.L(indu,indu));
    Wul = - G.L(indu,indl);
    tmp = (Wul * y(indl,:));
%     if ~param.exact
%         if ~isfield(param,'order'), param.order = 30; end
%         Gtemp.L = Luu;
%         Gtemp.N = size(Luu,1);
% %        Gtemp = gsp_estimate_lmax(Gtemp);
% %         cheb_coeffs = gsp_cheby_coeff(Gtemp, @(x) pinv_n(x,1e-8),...
% %         param.order, param.order +1);
% %         solt = gsp_cheby_op(Gtemp, cheb_coeffs, tmp);
%         paramt.method = 'lanczos';
%         paramt.order = param.order;
%         solt = gsp_filter_analysis(Gtemp, @(x) pinv_n(x,1e-8),tmp,paramt);
%
%     else
%        solt = pinv(full(Luu)) * tmp;
         solt = Luu \ tmp;
%     end
    sol = y;
    sol(indu,:) = solt;
    return
else

    %% prepare the graph
    if ~isfield(G,'lmax')
        G = gsp_estimate_lmax(G);
    end

    %% set the
    % setting the function f2 (see unlocbox for help)

    Mop =@(x) bsxfun(@times,M,x);

%     fg.grad = @(x) 2*Mop(Mop(x)-y);
%     fg.eval = @(x) norm(Mop(x)-y)^2;
%     fg.beta = 2;
% %     paramtik.verbose = param.verbose -1;
% %     ftik.prox = @(x,T) gsp_prox_tik(x,tau * T,G,paramtik);
%     ftik.eval = @(x) tau* sum(gsp_norm_tik(G,x));
%     ftik.grad = @(x) 2 * G.L * x;
%     ftik.beta = 2 * G.lmax;

%
%
% else
% %     param_b2.verbose = param.verbose -1;
% %     param_b2.y = y;
% %     param_b2.A = @(x) M.*x;
% %     param_b2.At = @(x) M.*x;
% %     param_b2.tight = param.tight;
% %     param_b2.epsilon = 0;
% %     fproj.prox = @(x,T) proj_b2(x,T,param_b2);
% %     fproj.eval = @(x) eps;
%
%     fproj.prox = @(x,T) x - Mop(x) + Mop(y);
%     fproj.eval = @(x) eps;
%     ftik.eval = @(x) sum(gsp_norm_tik(G,x));
%     Ltmp = G.L + G.L';
%     ftik.grad = @(x) Ltmp*x;
%     ftik.beta = 2*G.lmax;
% end

% %% solve the problem
%
% % setting different parameter for the simulation
% paramsolver = param;
%
% if tau > 0
%     sol = forward_backward(y,ftik,fg,paramsolver);
% else
%     sol = forward_backward(y,fproj,ftik,paramsolver);
% end

% sol = sol(logical(1-M));
%
% sol = reshape(sol,[],size(M,2));

    A = @(z)  vec ( Mop(reshape(z,[],size(y,2))) ...
        +  tau * G.L * reshape(z,[],size(y,2)));
    b =  vec(Mop(y));
    sol = pcg(A,b,param.tol,param.maxit,[],[],b);
    sol = reshape(sol,[],size(y,2));
end

end

% function r =  pinv_n(x,t)
%
% r = double(abs(x)>t) .* 1./x;
%
% end

%   * *param.direct* : Direct computation of the exact solution (only for
%     $\tau = 0$). (Default tau==0)
%   * *param.exact* : Exact computation of the exact solution (only for
%     $\tau = 0$ and param.direct = 0). (Default: (numel(M)-nnz(M))<1000 )

repo: gspbox


text:
%%
% Test function to compare nn_bnorm and its GPU/CPU implementation with
% using VLFEAT
%%

gpu = false;
gpu = true ;

T = 1 ;
x = randn(64,64,32,32,'single') ;
g = randn(32,1,'single') ;
b = randn(32,1,'single') ;

if gpu
  x = gpuArray(x) ;
  g = gpuArray(g) ;
  b = gpuArray(b) ;
end

a=vl_nnbnorm(x,g,b);
a_=vl_nnbnorm_old(x,g,b);

vl_testsim(a,a_)

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function cpv_test01 ( )

%*****************************************************************************80
%
%% CPV_TEST01 seeks the CPV of Integral ( -1 <= t <= 1 ) exp(t) / t dt
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 March 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CPV_TEST01:\n' );
  fprintf ( 1, '  CPV of Integral ( -1 <= t <= 1 ) exp(t) / t dt\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   N           Estimate             Error\n' );
  fprintf ( 1, '\n' );
  exact = 2.11450175075;
  a = -1.0;
  b = +1.0;
  for n = 2 : 2 : 8
    value = cpv ( @f01, a, b, n );
    fprintf ( 1, '  %2d  %24.16g  %14.6g\n', n, value, abs ( value - exact ) );
  end

  return
end
function value = f01 ( t )

%*****************************************************************************80
%
%% F01 evaluates the integrand of Integral ( -1 <= t <= 1 ) exp(t) / t dt
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T, the argument.
%
%    Output, real VALUE, the value of the integrand.
%
  value = exp ( t );

  return
end


text:
function [grad, elec] = mne2grad(hdr, dewar, coilaccuracy)

% MNE2GRAD converts a header from a fif file that was read using the MNE toolbox into
% a gradiometer structure that can be understood by the FieldTrip low-level forward
% and inverse routines.
%
% Use as
%   [grad, elec] = mne2grad(hdr, dewar, coilaccuracy)
% where
%   dewar        = boolean, whether to return it in dewar or head coordinates (default = false, i.e. head coordinates)
%   coilaccuracy = empty or a number (default = [])
%
% See also CTF2GRAD, BTI2GRAD

%%%%%%%%%%%%%%%%%% BEGIN REVISION HISTORY %%%%%%%%%%%%%%%%%%
%
% Robert Oostenveld 22/01/2016 reimplemented construction of integration
% points on basis of coil_def.dat file from MNE
%
% Laurence Hunt 03/12/2008 (with thanks to Joachim Gross's original script
% based on fiff_access). lhunt@fmrib.ox.ac.uk
%
% Teresa Cheung 09/24/2011 revision to Laurence Hunt's script. The coil
% geometry has been revised to better reflect the coil dimensions of the
% Vectorview (306) planar gradiometers and magnetometers. Coil geometry for
% the Neuromag-122 planar gradiometer is different and is not read with
% this script (coil_type should equal 2)
%
%%%%%%%%%%%%%%%%%% END REVISION HISTORY %%%%%%%%%%%%%%%%%%

% Copyrights (C) 2016, Robert Oostenveld
% Copyrights (C) 2011, Teresa Cheung
% Copyrights (C) 2008, Laurence Hunt
% Copyrights (C) 2005, Joachim Gross
%
% 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 nargin<2 || isempty(dewar)
  dewar = false;
end

if nargin<3 || isempty(coilaccuracy)
  % if empty it will use the original code
  % otherwise it will use the specified accuracy coil definition from the MNE coil_def.dat
  coilaccuracy = [];
end

% orig = fiff_read_meas_info(filename);
% orig = fiff_setup_read_raw(filename);

if isfield(hdr, 'orig')
  orig = hdr.orig; % use the original MNE header, not the FieldTrip header
else
  orig = hdr; % assume that it is the original header
end

% start with empty gradiometer definition
grad = [];

% start with empty electrode definition
elec = [];

% begin by transforming all channel locations into the desired coordinate system, if possible
if ~dewar
  if ~isempty(orig.dev_head_t)
    orig.chs = fiff_transform_meg_chs(orig.chs,orig.dev_head_t);
    orig.chs = fiff_transform_eeg_chs(orig.chs,orig.dev_head_t); % EEG channels are normally stored in head coordinates anyway, but what the heck
  else
    ft_warning('No device to head transform available in fif file');
    ft_warning('MEG channels will likely have coordinates in device frame, not head frame');
  end
else
  if ~isempty(orig.dev_head_t)
    % compute the transformation from head to device
    orig.head_dev_t.trans = inv(orig.dev_head_t.trans);
    orig.head_dev_t.from  = orig.dev_head_t.to;
    orig.head_dev_t.to    = orig.dev_head_t.from;
    orig.chs = fiff_transform_meg_chs(orig.chs,orig.head_dev_t); % MEG channels are normally stored in dewar coordinates anyway, but what the heck
    orig.chs = fiff_transform_eeg_chs(orig.chs,orig.head_dev_t);
  else
    ft_warning('No device to head transform available in fif file');
    ft_warning('EEG channels will likely have coordinates in head frame, not device frame');
  end
end

if ~isempty(coilaccuracy)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % use the coil definitions from the MNE coil_def.dat file
  % these allow for varying accuracy which is specified by
  % coilaccuracy = 0, 1 or 2
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  ft_hastoolbox('mne', 1);
  [ftver, ftpath] = ft_version;
  def = mne_load_coil_def(fullfile(ftpath, 'external', 'mne', 'coil_def.dat'));

  k = 1;

  grad.chanpos = nan(length(orig.chs),3);
  grad.chanori = nan(length(orig.chs),3);
  grad.label = cell(length(orig.chs),1);
  grad.tra      = nan(length(orig.chs),0);
  grad.chantype = repmat({'unknown'}, size(grad.label));
  grad.chanunit = repmat({'unknown'}, size(grad.label));

  for i=1:length(orig.chs)
    thisdef = def([def.id]==orig.chs(i).coil_type & [def.accuracy]==coilaccuracy);
    if isempty(thisdef)
      continue;
    end
    weight = thisdef.coildefs(:,1);
    pos = thisdef.coildefs(:,2:4);
    ori = thisdef.coildefs(:,5:7);
    H = orig.chs(i).coil_trans;
    pos = ft_warp_apply(H, pos);
    T = H;
    R = H;
    T(1:3,1:3) = 0; % remove the rotation, keep the translation
    R(1:3,4)   = 0; % remove the translation, keep the rotation
    ori = ft_warp_apply(R, ori); % only apply the rotation
    for j=1:thisdef.num_points
      grad.coilpos(k,:) = pos(j,:);
      grad.coilori(k,:) = ori(j,:);
      grad.tra(i,k)     = weight(j);
      k = k + 1;
    end
    grad.label{i} = orig.ch_names{i};
    grad.chanpos(i,:) = T(1:3,4);
    grad.chanori(i,:) = ft_warp_apply(R, [0 0 1]);
  end

  grad.label = grad.label(:);
  grad.unit  = 'm'; % the coil_def.dat file is in meter

  % all the chs.kinds and chs.coil_types are obtained from the MNE manual, p.210-211
  for sel=find([orig.chs.kind]==1 & [orig.chs.coil_type]==2)' % planar gradiometers
    grad.chantype(sel) = {'megplanar'}; %Neuromag-122 planar gradiometer
  end
  for sel=find([orig.chs.kind]==1 & [orig.chs.coil_type]==3012)' %planar gradiometers
    grad.chantype(sel) = {'megplanar'}; %Type T1 planar grad
  end
  for sel=find([orig.chs.kind]==1 & [orig.chs.coil_type]==3013)' %planar gradiometers
    grad.chantype(sel) = {'megplanar'}; %Type T2 planar grad
  end
  for sel=find([orig.chs.kind]==1 & [orig.chs.coil_type]==3014)' %planar gradiometers
    grad.chantype(sel) = {'megplanar'}; %Type T3 planar grad
  end
  for sel=find([orig.chs.kind]==1 & [orig.chs.coil_type]==3022)' %magnetometers
    grad.chantype(sel) = {'megmag'};    %Type T1 magenetometer
  end
  for sel=find([orig.chs.kind]==1 & [orig.chs.coil_type]==3023)' %magnetometers
    grad.chantype(sel) = {'megmag'};    %Type T2 magenetometer
  end
  for sel=find([orig.chs.kind]==1 & [orig.chs.coil_type]==3024)' %magnetometers
    grad.chantype(sel) = {'megmag'};    %Type T3 magenetometer
  end
  for sel=find([orig.chs.kind]==1 & [orig.chs.coil_type]==7001)' %axial gradiometer
    grad.chantype(sel) = {'megaxial'};
  end
  for sel=find([orig.chs.kind]==301)' %MEG reference channel, located far from head
    grad.chantype(sel) = {'ref'};
  end
  for sel=find([orig.chs.kind]==2)'   %EEG channels
    grad.chantype(sel) = {'eeg'};
  end
  for sel=find([orig.chs.kind]==201)' %MCG channels
    grad.chantype(sel) = {'mcg'};
  end
  for sel=find([orig.chs.kind]==3)' %Stim channels
    if any(ismember([orig.chs(sel).logno], [101 102])) % new systems: 101 (and 102, if enabled) are digital; low numbers are 'pseudo-analog' (if enabled)
      grad.chantype(sel([orig.chs(sel).logno] == 101)) = {'digital trigger'};
      grad.chantype(sel([orig.chs(sel).logno] == 102)) = {'digital trigger'};
      grad.chantype(sel([orig.chs(sel).logno] <= 32))  = {'analog trigger'};
      others = [orig.chs(sel).logno] > 32 & [orig.chs(sel).logno] ~= 101 & ...
        [orig.chs(sel).logno] ~= 102;
      grad.chantype(sel(others)) = {'other trigger'};
    elseif any(ismember([orig.chs(sel).logno], [14 15 16])) % older systems: STI 014/015/016 are digital; lower numbers 'pseudo-analog'(if enabled)
      grad.chantype(sel([orig.chs(sel).logno] == 14)) = {'digital trigger'};
      grad.chantype(sel([orig.chs(sel).logno] == 15)) = {'digital trigger'};
      grad.chantype(sel([orig.chs(sel).logno] == 16)) = {'digital trigger'};
      grad.chantype(sel([orig.chs(sel).logno] <= 13)) = {'analog trigger'};
      others = [orig.chs(sel).logno] > 16;
      grad.chantype(sel(others)) = {'other trigger'};
    else
      ft_warning('There does not seem to be a suitable trigger channel.');
      grad.chantype(sel) = {'other trigger'};
    end
  end
  for sel=find([orig.chs.kind]==202)' %EOG
    grad.chantype(sel) = {'eog'};
  end
  for sel=find([orig.chs.kind]==302)' %EMG
    grad.chantype(sel) = {'emg'};
  end
  for sel=find([orig.chs.kind]==402)' %ECG
    grad.chantype(sel) = {'ecg'};
  end
  for sel=find([orig.chs.kind]==502)' %MISC
    grad.chantype(sel) = {'misc'};
  end
  for sel=find([orig.chs.kind]==602)' %Resp
    grad.chantype(sel) = {'respiration'};
  end

  % FIFF.FIFF_UNIT_HZ  = 101;
  % FIFF.FIFF_UNIT_N   = 102;
  % FIFF.FIFF_UNIT_PA  = 103;
  % FIFF.FIFF_UNIT_J   = 104;
  % FIFF.FIFF_UNIT_W   = 105;
  % FIFF.FIFF_UNIT_C   = 106;
  % FIFF.FIFF_UNIT_V   = 107;
  % FIFF.FIFF_UNIT_F   = 108;
  % FIFF.FIFF_UNIT_OHM = 109;
  % FIFF.FIFF_UNIT_MHO = 110;
  % FIFF.FIFF_UNIT_WB  = 111;
  % FIFF.FIFF_UNIT_T   = 112;
  % FIFF.FIFF_UNIT_H   = 113;
  % FIFF.FIFF_UNIT_CEL = 114;
  % FIFF.FIFF_UNIT_LM  = 115;
  % FIFF.FIFF_UNIT_LX  = 116;

  for sel=find([orig.chs.unit]==107)'
    grad.chanunit(sel) = {'V'};
  end
  for sel=find([orig.chs.unit]==112)'
    grad.chanunit(sel) = {'T'};
  end
  for sel=find([orig.chs.unit]==201)'
    grad.chanunit(sel) = {'T/m'};
  end

  remove = cellfun(@isempty, grad.label);
  grad.label    = grad.label(~remove);
  grad.chantype = grad.chantype(~remove);
  grad.chanunit = grad.chanunit(~remove);
  grad.tra      = grad.tra(~remove,:);
  grad.chanpos  = grad.chanpos(~remove,:);
  grad.chanori  = grad.chanori(~remove,:);

else
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % use the original implementation
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % how many Planar gradiometers?
  nPlaGrad = 0;
  for i = 1:orig.nchan
    nPlaGrad = nPlaGrad +(orig.chs(i).coil_type==3012 | orig.chs(i).coil_type==3013 | orig.chs(i).coil_type==3014 | orig.chs(i).coil_type==2);
  end

  % how many Magnetometers?
  nMag = 0;
  for i = 1:orig.nchan
    nMag = nMag +(orig.chs(i).coil_type==3022 | orig.chs(i).coil_type==3023 | orig.chs(i).coil_type==3024);
  end

  % how many Axial gradiometers?
  nAxGrad = 0;
  for i = 1:orig.nchan
    nAxGrad = nAxGrad +(orig.chs(i).coil_type==7001); % babySQUID
  end

  % % how many IAS (internal active shielding) channels?
  % nIAS = 0;
  % for i = 1:orig.nchan;
  %   nIAS = nIAS + ~isempty(strmatch('IAS', orig.chs(i).ch_name));
  % end

  % how many magnetic field sensors in total?
  nSensors = nPlaGrad + nMag + nAxGrad;

  % how many coils in total?
  nCoils = nPlaGrad*2 + nMag + nAxGrad*2;

  % intialise grad structure
  grad.coilpos  = zeros(nCoils,3);
  grad.coilori  = zeros(nCoils,3);
  grad.tra      = zeros(nSensors,nCoils);
  grad.unit     = 'cm'; % see below for the conversion, the original fif units are in meter
  grad.label    = cell(nSensors,1);

  if dewar
    grad.coordsys = 'dewar';
  else
    grad.coordsys = 'neuromag';
  end

  % initialise elec structure, this can remain empty
  elec = [];

  % define coils
  kCoil = 1;
  kChan = 1;
  % cf. Joachim's original script - I've implemented it this way in case MEG
  % and EEG channels are not listed first in the .fif file; this shouldn't
  % ever be the case but acts as a safety net...

  for n = 1:orig.nchan
    if (orig.chs(n).coil_type==3022 || orig.chs(n).coil_type==3023 || orig.chs(n).coil_type==3024) % magnetometer
      t = orig.chs(n).coil_trans;

      % TC 2011 09 24 I have changed the coil definition, the original was
      % grad.coilpos(kCoil,:) = t(1:3,4);

      grad.coilpos(kCoil,:) = t(1:3,4)+0.0003*t(1:3,3);
      grad.coilori(kCoil,:) = t(1:3,3);
      grad.tra(kChan,kCoil) = 1;
      kCoil = kCoil+1;
      grad.label{kChan} = deblank(orig.ch_names{n});
      grad.chantype{kChan,1}='megmag';
      kChan = kChan+1;

    elseif (orig.chs(n).coil_type==3012 || orig.chs(n).coil_type==3013 || orig.chs(n).coil_type==3014 || orig.chs(n).coil_type==2) % planar gradiometer
      t = orig.chs(n).coil_trans;

      % TC 2011 09 24 I have changed the coil definition, the original was
      % grad.coilpos(kCoil,:) = 0.0000*t(1:3,3)+t(1:3,4)-0.0084*t(1:3,1)); % for the 1st coil
      % grad.coilpos(kCoil,:) = 0.0000*t(1:3,3)+t(1:3,4)+0.0084*t(1:3,1)); % for the 2nd coil

      grad.coilpos(kCoil,:) = 0.0003*t(1:3,3)+t(1:3,4)-0.0084*t(1:3,1); % for the 1st coil
      grad.coilori(kCoil,:) = t(1:3,3);
      grad.tra(kChan,kCoil) = -1;
      kCoil = kCoil+1;

      grad.coilpos(kCoil,:) = 0.0003*t(1:3,3)+t(1:3,4)+0.0084*t(1:3,1); % for the 2nd coil
      grad.coilori(kCoil,:) = t(1:3,3);
      grad.tra(kChan,kCoil) = 1;
      kCoil = kCoil+1;

      grad.label{kChan} = deblank(orig.ch_names{n});
      grad.chantype{kChan,1}='megplanar';
      kChan = kChan+1;

    elseif (orig.chs(n).coil_type==7001)  % babySQUID axial gradiometer bottom coils
      t = orig.chs(n).coil_trans;

      grad.coilpos(kCoil,:)=t(1:3,4);  % for the 1st coil
      grad.coilori(kCoil,:)=t(1:3,3);
      grad.tra(kChan,kCoil)=1;
      kCoil=kCoil+1;

      grad.coilpos(kCoil,:)=t(1:3,4)+0.050*t(1:3,3);  % for the 2nd coil
      grad.coilori(kCoil,:)=t(1:3,3);
      grad.tra(kChan,kCoil)=-1;
      kCoil=kCoil+1;

      grad.label{kChan}=deblank(orig.ch_names{n});
      grad.chantype{kChan,1}='megaxial';
      kChan=kChan+1;

    else
      % do nothing - either an EEG channel or something else such as a stim channel
    end
  end

  % check we've got all the MEG channels:
  kChan = kChan-1;
  if kChan ~= (nPlaGrad + nMag +nAxGrad)
    ft_error('Number of MEG channels identified does not match number of channels in grad structure');
  end

  % determine the type of acquisition system
  if nAxGrad>0
    grad.type = 'babysquid74';
  elseif nPlaGrad>122 && nMag~=0
    grad.type = 'neuromag306';
  elseif nPlaGrad<=122 && nMag==0
    grad.type = 'neuromag122';
  else
    % do not specify type of acquisition system
  end

  % multiply by 100 to get cm
  grad.coilpos = 100*grad.coilpos;
end

% how many EEG channels?
nEEG = 0;
for i = 1:orig.nchan
  nEEG = nEEG +(orig.chs(i).kind==2);
end

% define EEG channels
if nEEG>0
  elec.elecpos  = zeros(nEEG,3);
  elec.unit     = 'cm';
  elec.label    = cell(nEEG,1);

  % Amendments to overcome problem with fiff_read_meas_info.m when >60 channels (thanks to Rik Henson)

  dig_eeg = find([orig.dig.kind]==3); % Find EEG digitisations

  if ~isempty(dig_eeg)
    dig_eeg(1) = []; % Remove reference
  end

  chn_eeg = find([orig.chs.kind]==2); % Find EEG channels

  kChan = 0;
  for n = chn_eeg
    kChan = kChan+1;
    if nEEG<=60
      elec.elecpos(kChan,1:3) = orig.chs(n).eeg_loc(1:3);
    else
      if kChan<=numel(dig_eeg)
        elec.elecpos(kChan,1:3) = orig.dig(dig_eeg(kChan)).r;
      else
        ft_warning('not all EEG channel positions have been digitized');
        elec.elecpos(kChan,1:3) = nan;
      end
    end
    elec.label{kChan} = deblank(orig.ch_names{n});
  end

  % check we've got all the EEG channels:
  if kChan~=(nEEG)
    ft_error('Number of EEG channels identified does not match number of channels in elec structure!!!!!');
  end

  % multiply by 100 to get cm
  elec.elecpos = 100*elec.elecpos;
end

% construct SSP projectors, these can be applied later to the data and grad
% structure using FT_DENOISE_SSP
grad.balance = [];
grad.balance.current = 'none';

for i = 1:length(hdr.projs)
  hdr.projs(i).active = 'true';
  montage = [];
  montage.tra = mne_make_projector(hdr.projs(i), hdr.ch_names, []);
  montage.labelorg = hdr.ch_names;
  montage.labelnew = hdr.ch_names;
  grad.balance.(fixname(hdr.projs(i).desc)) = montage;
end

% remove grad if completely empty
if size(grad.label,1) == 0
  grad = [];
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [align] = planarmove(in,config)
% planarmove ... alignment under assumption of planar motion
%
% [align] = planarmove(in,config)
% in, cam, config ... see the main GOCAL script
%
% align ... structures aligned wit the specified world frame
%
% $Id: planarmove.m,v 1.2 2005/05/20 15:31:31 svoboda Exp $

% fit a plane to the reconstructed points and estimate normal

plane.n = planefit(in.Xe(1:3,:)');

new.n = [0,0,1]'; % align the xy plane horizontally

rotaxis = cross(plane.n,new.n);
rotangle = acos( (plane.n'*new.n)/norm(plane.n)*norm(new.n) );

R = nfi2r(rotaxis,rotangle);
s = 3;
t = [0,0,1]' - s*R*mean(in.Xe(1:3,:)')';

align.simT.s = s;
align.simT.R = R;
align.simT.t = t;

[align.P, align.X]							= align3d(in.Pe,in.Xe,align.simT);
% save aligned data
if 1 % SAVE_STEPHI | SAVE_PGUHA
	[align.Cst,align.Rot] = savecalpar(align.P,config);
end
drawscene(align.X,align.Cst',align.Rot,61,'cloud','Graphical Output Validation: Aligned data, TopView',config.cal.cams2use);

set(gca,'CameraTarget',[0,0,1]);
set(gca,'CameraPosition',[0,0,2]);

figure(61),
% print -depsc graphevalaligned.eps
eval(['print -depsc ', config.paths.data, 'topview.eps'])

drawscene(align.X,align.Cst',align.Rot,62,'cloud','Graphical Output Validation: Aligned data, SideView',config.cal.cams2use);

set(gca,'CameraTarget',[0,0,0.9]);
set(gca,'CameraPosition',[2,0,0.9]);

figure(62),
% print -depsc graphevalaligned.eps
eval(['print -depsc ', config.paths.data, 'sideview.eps'])

return

author: strawlab

repo: MultiCamSelfCal


text:
function M1 = spm_eeg_inv_datareg(S)
% Co-registration of two sets of fiducials according to sets of
% corresponding points and (optionally) headshapes.
% rigid co-registration
%           1: fiducials based (3 landmarks: nasion, left ear, right ear)
%           2: surface matching between sensor mesh and headshape
%           (starts with a type 1 registration)
%
% FORMAT M1 = spm_eeg_inv_datareg(S)
%
% Input:
%
% S  - input struct
% fields of S:
%
% S.sourcefid  - EEG fiducials (struct)
% S.targetfid = MRI fiducials
% S.template  - 1 - input is a template (for EEG)
%               0 - input is an individual head model
%               2 - input is a template (for MEG) - enforce uniform scaling
%
% S.useheadshape - 1 use headshape matching 0 - don't
%
%
% Output:
% M1 = homogenous transformation matrix
%
% If a template is used, the senor locations are transformed using an
% affine (rigid body) mapping.  If headshape locations are supplied
% this is generalized to a full twelve parameter affine mapping (n.b.
% this might not be appropriate for MEG data).
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

% Jeremie Mattout
% $Id: spm_eeg_inv_datareg.m 7544 2019-03-15 16:20:16Z vladimir $

if ~isfield(S, 'targetfid')
    error('Target fiducials are missing');
else
    targetfid = ft_convert_units(S.targetfid, 'mm');
end

if ~isfield(S, 'sourcefid')
    error('Source are missing');
else
    sourcefid = ft_convert_units(S.sourcefid, 'mm');
    [sel1, sel2] = spm_match_str(targetfid.fid.label, sourcefid.fid.label);
    sourcefid.fid.pnt = sourcefid.fid.pnt(sel2, :);
    sourcefid.fid.label = sourcefid.fid.label(sel2);

    targetfid.fid.pnt = targetfid.fid.pnt(sel1, :);
    targetfid.fid.label = targetfid.fid.label(sel1);
end

if ~isfield(S, 'template')
    S.template = 0;
end

% Estimate-apply rigid body transform to sensor space
%--------------------------------------------------------------------------
M1 = spm_eeg_inv_rigidreg(targetfid.fid.pnt', sourcefid.fid.pnt');

sourcefid = ft_transform_geometry(M1, sourcefid);

if S.template

    % constrained affine transform
    %--------------------------------------------------------------------------
    aff   = S.template;
    for i = 1:64

        % scale
        %----------------------------------------------------------------------
        M       = pinv(sourcefid.fid.pnt(:))*targetfid.fid.pnt(:);
        M       = sparse(1:4,1:4,[M M M 1]);

        sourcefid = ft_transform_geometry(M, sourcefid);

        M1      = M*M1;

        % and move
        %----------------------------------------------------------------------
        M       = spm_eeg_inv_rigidreg(targetfid.fid.pnt', sourcefid.fid.pnt');

        sourcefid = ft_transform_geometry(M, sourcefid);

        M1      = M*M1;

        if (norm(M)-1)< eps
            break;
        end
    end
else
    aff = 0;
end

% Surface matching between the scalp vertices in MRI space and the
% headshape positions in data space
%--------------------------------------------------------------------------
if ~isempty(sourcefid.pnt) && S.useheadshape

    headshape = sourcefid.pnt;
    scalpvert = targetfid.pnt;

    % load surface locations from sMRI
    %----------------------------------------------------------------------
    if size(headshape,2) > size(headshape,1)
        headshape = headshape';
    end
    if size(scalpvert,2) > size(scalpvert,1)
        scalpvert = scalpvert';
    end

    % intialise plot
    %----------------------------------------------------------------------
    if ~spm('CmdLine')
        h    = spm_figure('GetWin','Graphics');
        spm_figure('Select',h);
        spm_clf(h);
        Fmri = plot3(scalpvert(:,1),scalpvert(:,2),scalpvert(:,3),'ro','MarkerFaceColor','r');
        hold on;
        Fhsp = plot3(headshape(:,1),headshape(:,2),headshape(:,3),'bs','MarkerFaceColor','b');
        axis off image
        drawnow
    else
        Fmri = [];
        Fhsp = [];
    end

    % nearest point registration
    %----------------------------------------------------------------------
    M = spm_eeg_inv_icp(scalpvert',headshape',targetfid.fid.pnt',sourcefid.fid.pnt',Fmri,Fhsp,aff);

    % transform headshape and eeg fiducials
    %----------------------------------------------------------------------
    sourcefid = ft_transform_geometry(M, sourcefid);
    M1        = M*M1;
end

author: spm


text:
function [X, R] = imstack2vectors(S, MASK)
[M, N, n] = size(S);
if nargin == 1
   MASK = true(M, N);
else
   MASK = MASK ~= 0;
end

[I, J] = find(MASK);
R = [I, J];

Q = M*N;
X = reshape(S, Q, n);

MASK = reshape(MASK, Q, 1);

X = X(MASK, :);


text:
function [ts] = timestamp_neuralynx(tsl, tsh)

% TIMESTAMP_NEURALYNX merge the low and high part of Neuralynx timestamps
% into a single uint64 value

% Copyright (C) 2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: timestamp_neuralynx.m 2885 2011-02-16 09:41:58Z roboos $

if ~isa(tsl, 'uint32') && ~isa(tsl, 'int32')
  error('invalid input');
elseif ~isa(tsh, 'uint32') && ~isa(tsl, 'int32')
  error('invalid input');
end

% convert the 32 bit low and 32 bit high timestamp into a 64 bit integer
dum = zeros(2, length(tsh), 'uint32');
if littleendian
  dum(1,:) = tsl;
  dum(2,:) = tsh;
else
  dum(1,:) = tsh;
  dum(2,:) = tsl;
end

ts = typecast(dum(:), 'uint64');

author: PatternRecognition

repo: OpenBMI


text:
function cdf = logistic_cdf ( x, a, b )

%*****************************************************************************80
%
%% LOGISTIC_CDF evaluates the Logistic CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real CDF, the value of the CDF.
%
  cdf = 1.0 / ( 1.0 + exp ( ( a - x ) / b ) );

  return
end


text:
%Test script for the ARGUS method

clear; %close all;
%path added for regression methods
%Set seed for repeatability
%s = RandStream('mt19937ar','Seed',1);
%RandStream.setGlobalStream(s);
rand('seed',1)

%Set initial conditions
intialSize =5;
dim = 2;
%Set stop conditions
maxPoints = 1000;
pntsAdd = 10;

%Set function handle for loop (not needed for ARGUS)
stdr = .5;
%fun = @(x)sin(6*x)+stdr*randn(length(x),1);
fun = @(x)sin(6*sum(x,2))+stdr*randn(length(x(:,1)),1);
range = 0:.01:1;
[x1 x2] = meshgrid(range);
y = sin(6*(x1(:)+x2(:)));
subplot(2,2,2)
hold off;
surf(x1,x2,reshape(y,size(x1)),'EdgeColor',[.5 .5 .5],'FaceAlpha',0);
hold on;
subplot(2,2,4)
hold off;
surf(x1,x2,stdr*ones(size(x2)),'EdgeColor',[.5 .5 .5],'FaceAlpha',0);
hold on;
%plot(0:.01:1,sin(6*(0:.01:1)'))

%Initial conditions
x = lhsdesign(intialSize,dim);
y= feval(fun,x);
XYdata = add2XYdata([],x,y);
subplot(2,2,4)
scatter3(x(:,1),x(:,2),zeros(size(x(:,1))),'rs')
subplot(2,2,2)
scatter3(x(:,1),x(:,2),y,'rs')
subplot(2,2,1)
hold off
scatter(x(:,1),x(:,2),'rs')
hold on
%plot(x,y,'rs')

%Work horse
i=intialSize;
[newpoints hist] = ARGUS(XYdata,pntsAdd);
ynew = feval(fun,newpoints);
XYdata = add2XYdata(XYdata,newpoints,ynew);
%plot(newpoints,ynew,'ks')
%text(newpoints,ynew,['\leftarrow' sprintf('%i',i)]);
brk =floor(pntsAdd*hist.ratio(end));
subplot(2,2,1)
scatter(newpoints(1:brk,1),newpoints(1:brk,2),'ko')
scatter(newpoints(brk+1:end,1),newpoints(brk+1:end,2),'bx')
%legend('Initial Deisgn','Replication','Exploration','Location','Best')
%legend('boxoff')
h_oldstErr=surf(x1,x2,reshape(feval(hist.stErrfun,[x1(:),x2(:)]),size(x1)),'Edgecolor','none','FaceAlpha',.2);
colorbar
xlabel('X1')
ylabel('X2')
title('Standard Error')
subplot(2,2,3)
plot(1:hist.itt,floor(pntsAdd*hist.ratio)/pntsAdd,'bs-')
xlabel('Iteration')
ylabel('Replications/Exploration')
subplot(2,2,2)
title('Mean')
scatter3(newpoints(1:brk,1),newpoints(1:brk,2),ynew(1:brk),'ko')
scatter3(newpoints(brk+1:end,1),newpoints(brk+1:end,2),ynew(brk+1:end),'bx')
h_oldmu = surf(x1,x2,reshape(feval(hist.mufun,[x1(:),x2(:)]),size(x1)),'Edgecolor','none','FaceAlpha',.5);
%text(newpoints(:,1),newpoints(:,2),ynew,['\leftarrow' sprintf('%i',i)]);
subplot(2,2,4)
title('Standard Deviation')
scatter3(newpoints(1:brk,1),newpoints(1:brk,2),zeros(size(newpoints(1:brk,1))),'ko')
scatter3(newpoints(brk+1:end,1),newpoints(brk+1:end,2),zeros(size(newpoints(brk+1:end,1))),'bx')
h_oldsig = surf(x1,x2,reshape(feval(hist.stdfun,[x1(:),x2(:)]),size(x1)),'Edgecolor','none','FaceAlpha',.5);
%text(newpoints(:,1),newpoints(:,2),zeros(size(newpoints(:,1))),['\leftarrow' sprintf('%i',i)]);
i=intialSize+pntsAdd;
while maxPoints>i
    [newpoints hist] = ARGUS(XYdata,pntsAdd,hist);
    ynew = feval(fun,newpoints);
    XYdata = add2XYdata(XYdata,newpoints,ynew);
    %plot(newpoints,ynew,'ks')
    brk =floor(pntsAdd*hist.ratio(end));
    subplot(2,2,1)
    scatter(newpoints(1:brk,1),newpoints(1:brk,2),'ko')
    scatter(newpoints(brk+1:end,1),newpoints(brk+1:end,2),'bx')
    h_newstErr=surf(x1,x2,reshape(feval(hist.stErrfun,[x1(:),x2(:)]),size(x1)),'Edgecolor','none','FaceAlpha',.2);
    %colorbar
    delete(h_oldstErr);
    h_oldstErr=h_newstErr;
    subplot(2,2,3)
    plot(1:hist.itt,floor(pntsAdd*hist.ratio)/pntsAdd,'bs-')
    xlabel('Iteration')
    ylabel('Replications/Exploration')
    subplot(2,2,2)
    scatter3(newpoints(1:brk,1),newpoints(1:brk,2),ynew(1:brk),'ko')
    scatter3(newpoints(brk+1:end,1),newpoints(brk+1:end,2),ynew(brk+1:end),'bx')
    h_newmu = surf(x1,x2,reshape(feval(hist.mufun,[x1(:),x2(:)]),size(x1)),'Edgecolor','none','FaceAlpha',.5);
    %text(newpoints(:,1),newpoints(:,2),ynew,['\leftarrow' sprintf('%i',i)]);
    delete(h_oldmu);
    h_oldmu=h_newmu;
    subplot(2,2,4)
    scatter3(newpoints(1:brk,1),newpoints(1:brk,2),zeros(size(newpoints(1:brk,1))),'ko')
    scatter3(newpoints(brk+1:end,1),newpoints(brk+1:end,2),zeros(size(newpoints(brk+1:end,1))),'bx')
    h_newsig = surf(x1,x2,reshape(feval(hist.stdfun,[x1(:),x2(:)]),size(x1)),'Edgecolor','none','FaceAlpha',.5);
    %text(newpoints(:,1),newpoints(:,2),zeros(size(newpoints(:,1))),['\leftarrow' sprintf('%i',i)]);
    delete(h_oldsig);
    h_oldsig=h_newsig;
    %text(newpoints,ynew,['\leftarrow' sprintf('%i',i)]);
    i=i+pntsAdd;
    if mod(i-intialSize,3*pntsAdd)==0
        fprintf('Done with %i\n',i);end
end


text:
classdef IMMOEA_F9 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing IM-MOEA

%------------------------------- Reference --------------------------------
% R. Cheng, Y. Jin, K. Narukawa, and B. Sendhoff, A multiobjective
% evolutionary algorithm using Gaussian process-based inverse modeling,
% IEEE Transactions on Evolutionary Computation, 2015, 19(6): 838-856.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = [1,zeros(1,obj.D-1)+10];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            t = X(:,2:obj.D).^(1./(1+3*repmat(2:obj.D,size(X,1),1)/obj.D)) - repmat(X(:,1),1,obj.D-1);
            g = sum(t.^2/4000,2) - prod(cos(t./repmat(sqrt(1:obj.D-1),size(X,1),1)),2) + 2;
            PopObj(:,1) = X(:,1);
            PopObj(:,2) = g.*(1-sqrt(PopObj(:,1)./g));
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
classdef CellNodesDescriptor < handle

    properties (SetAccess = private)
        nodesInXmin
        nodesInXmax
        nodesInYmin
        nodesInYmax
        cornerNodes
    end

    properties (Access = private)
        x
        y
        xmin
        xmax
        ymin
        ymax
        allNodes
    end

    methods (Access = public)

        function obj = CellNodesDescriptor(coord)
            obj.init(coord);
            obj.obtainCellLimits();
            obj.obtainNodesInCellFaces();
            obj.obtainCornerNodes();
        end

    end

    methods (Access = private)

        function init(obj,coord)
            obj.x = coord(:,1);
            obj.y = coord(:,2);
            obj.allNodes(:,1) = 1:size(obj.x,1);
        end

        function obtainCellLimits(obj)
            obj.xmin = min(obj.x);
            obj.xmax = max(obj.x);
            obj.ymin = min(obj.y);
            obj.ymax = max(obj.y);
        end

        function obtainNodesInCellFaces(obj)
            obj.nodesInXmin = obj.getNodesInCellFace(obj.x,obj.xmin);
            obj.nodesInXmax = obj.getNodesInCellFace(obj.x,obj.xmax);
            obj.nodesInYmin = obj.getNodesInCellFace(obj.y,obj.ymin);
            obj.nodesInYmax = obj.getNodesInCellFace(obj.y,obj.ymax);
        end

        function obtainCornerNodes(obj)
            XminYmin = intersect(obj.nodesInXmin,obj.nodesInYmin);
            XminYmax = intersect(obj.nodesInXmin,obj.nodesInYmax);
            XmaxYmin = intersect(obj.nodesInXmax,obj.nodesInYmin);
            XmaxYmax = intersect(obj.nodesInXmax,obj.nodesInYmax);
            obj.cornerNodes = [XminYmin XminYmax XmaxYmin XmaxYmax]';
        end

        function nodesInFace = getNodesInCellFace(obj,x,xLim)
            dist2face = abs(x - xLim);
            isInFace = dist2face < 1e-13;
            nodesInFace = obj.allNodes(isInFace);
        end
    end

end

author: SwanLab

repo: Swan


text:
function [img,hdr] = makeAnalyzeFromIfiles(ifileDir, outFileName, options,startIndex)
% [img,hdr] = makeAnalyzeFromIfiles(ifileDir, [outFileName], [options],[startIndex])
%
% ifileDir shall include the filename prefix- eg. 'spgr1/I' for I-files
% that look like 'I.001', 'I.002.dcm', etc...
%
% outFileName should not include '.hdr' or '.img'- just the base name.
% It defaults to a reasonable file name in the current working directory.
%
% options is a cell array specifying various little options. Currently
% supported options:
%   'verbose' : displays various bits of information and progress (default).
%   'silent'  : opposite of verbose.
%
% We will grab all suitable I-Files in ifileDir, the startIndex input is
% removed. If you did acquire many anatomies using the same prescription,
% you shall put them in separate folders respectively. -- Junjie
%
% The data in the resulting Analyze file will be oriented according to
% the default Analyze convention- 'transverse unflipped' (orient code '0' in
% the .hdr file). That is, the image will be reoriented so that left-right is
% along the x-axis, anterior-posterior is along the y-axis, superior-inferior
% is along the z-axis, and the leftmost, anterior-most, superior-most point is
% at 0,0,0 (which, for Analyze, is the lower left-hand corner of the last slice).
%
% To help you find the 0,0,0 point, a 4-pixel rectange is drawn there with a pixel
% value equal to the maximum image intensity.
%
% The reorientation involves only cannonical rotations and mirror flips- no
% interpolation is performed. Also, note that the reorientation depends on info
% in the GE I-file header- if this info is wrong, the reorientation will be wrong.
%
% REQUIRES:
%   * Stanford Anatomy and filter functions
%
% SEE ALSO: loadAnalyze, saveAnalyze, analyze2mrGray, mrGray2Analyze,
% makeAnalyzeFromRaw (that does all the hard work).
%
% HISTORY:
%   2002.05.31 RFD (bob@white.stanford.edu): wrote it.
%   2002.06.05 RFD: got it to actually work reliably.
%   2002.06.14 RFD: cleaned up the code and added more comments.
%   2003.06.19 RFD: Moved logic for computing the transform to the
%   cannonical pixel order from makeAnalyzeFromRaw to a separate function
%   (computeCannonicalXformFromIfile). What does that have to do with this
%   function? Well, in doing so, things got simple enough that we can just
%   do everything in here.
%   2004.01.24 Junjie: disable the startIndex input, with explanation.
%   2004.06.17 RFD: optionally reurns img and hdr.

if(~exist('ifileDir','var') | isempty(ifileDir))
    [f,p] = uigetfile({'*.dcm','DICOM';'*.*','All files'}, 'Select the first raw image file...');
    if(isnumeric(f))
        disp('User canceled.'); return;
    end
    ifileDir = fullfile(p,f);
    if(~exist('outFileName','var') | isempty(outFileName))
        [p,f] = fileparts(ifileDir);
        outFileName = fullfile(pwd, f);
    end
end
if(~ischar(ifileDir))
    help(mfilename);
    return;
end

if(exist('outFileName','file'))
    error([outFileName 'exists!']);
end
if(~exist('options','var'))
    options = {};
end
% if(~exist('startIndex','var'))
%     startIndex=1;
% end
% indSuffix = sprintf('%03d',startIndex);
if exist('startIndex','var');
    warning('startIndex no longer a valid input. Each folder shall contain only one series of anatomies, so that the code auto detect first and last slices');
end

[img2std, ifileBaseName, mmPerVox, imDim, notes] = computeCannonicalXformFromIfile(ifileDir);

% Load the image data
disp(['Loading data from ', ifileDir, '...']);
img = int16(makeCubeIfiles(ifileDir, imDim([1,2]), 1:imDim(3))); %[startIndex:startIndex+imDim(3)-1]));

% The following does all the hard work.
%makeAnalyzeFromRaw(img, [ifileDir,'.001'], sprintf('%s.%03d',ifileDir,nSlices), nSlices, outFileName, options);
[img, mmPerVoxNew] = applyCannonicalXform(img, img2std, mmPerVox);

disp(['Saving ',outFileName,'.hdr and ',outFileName,'.img ...']);
hdr = saveAnalyze(img, outFileName, mmPerVoxNew, notes);

if(nargout<1)
    clear img hdr;
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%--- help for decompose_parameter_name ---
%
%  INTERNAL FUNCTION
%
%

author: jmaih