function Spath = Simulate_Jump_Diffusion_func( N_sim, M, T, S_0, r, q, sigma, jumpModel, jumpParams)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Simulates Paths of Jump Diffusion Models with jumps (including simple Black-Scholes, no jumps)
%        Uses log-Euler scheme
% Returns: paths of dimension (N_sim, M+1), since they include S_0
%          ... Simulates N_sim paths, each row is a full path starting from S_0, ending with S_M (M+1 points in path)
% Author: Justin Lars Kirkby
%
% -----------------
% Params
% -----------------
% N_sim = # paths
% M = #time steps on [0,T], time step is dt=T/M, so each path has M+1 points
% T = time to maturity, ie path is on [0,T]
% S_0 = initial underlying value (e.g. S_0=100)
% r = interst rate (e.g. r = 0.05)
% q = dividend yield (e.g. q = 0.05)
% sigma = diffusion parameter (e.g. sigma = 0.2)
%
%===================================
% jumpModel: 0 = NoJumps, 1 = NormalJumps, 2 = DEJumps, 3 = MixedNormalJumps
%===================================
% jumpParams = paramters container containing all necessary params for models,
%            : if jumpModel = 0, no jump params are needed
%            : if jumpModel > 0, jumpParams must contain lambda, kappa, and any other model specific params (see below)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin < 8   % Then default is standard diffusion, NO JUMPS
    jumpModel = 0;  jumpParams = {};
end

dt = T/M;

%==============================
% Initialize Jump Model Params and JumpFunc (function handle)
%==============================
%%% NOTE:  Jump Model is of the form in LOG space
%%% X(m+1) = X(m) + drift + Brownian Component + sum(Jumps on [m,m+1])
%%% By Jump we mean log(Y), e.g. in Merton Model, Jump ~ Normal (since we are in log space )

if jumpModel > 0 %ie if there are jumps in the model
    lambda = jumpParams.lambda;
    kappa  = jumpParams.kappa;

    Zeta = r - q - lambda*kappa;  %NOTE: we are redefining r to include compensation for jump component
    lamdt = lambda*dt;

    if jumpModel == 1 %Normal Jumps, e.g. Merton
        muJ  = jumpParams.muJ;
        sigJ = jumpParams.sigJ;
        JumpFunc = @(n) sum(muJ +sigJ*randn(n,1)); %Generates n independent jumps and sums them

    elseif jumpModel == 2 %Double Exponenial Jumps
        p_up = jumpParams.p_up;
        eta1 = jumpParams.eta1;
        eta2 = jumpParams.eta2;
        JumpFunc = @(n) sum(DoubleExpoRnd(n,p_up, eta1,eta2));

    elseif jumpModel == 3 %Mixed normal Jumps
        p_up = jumpParams.p_up;
        a1 = jumpParams.a1;  b1 = jumpParams.b1;
        a2 = jumpParams.a2;  b2 = jumpParams.b2;
        JumpFunc = @(n) sum(MixedNormalRnd(n, p_up, a1,b1,a2,b2));
    end
else
    Zeta = r - q;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Spath = zeros(N_sim,M+1);
Spath(:,1) = S_0;
Sigsqdt = sigma*sqrt(dt);
drift = (Zeta - .5*sigma^2)*dt;

if jumpModel == 0
    for m = 1:M
        W1 = randn(N_sim,1);
        Spath(:,m+1) = Spath(:,m).*exp(drift + Sigsqdt*W1);  %log scheme
    end
else
    for m = 1:M
        Poi = PoissonRnd(N_sim, lamdt);  %Generate Poisson Column Vector of size N_Sim
        sumJumpsVec = zeros(N_sim,1);
        for n = 1:N_sim
            if Poi(n)>0
                sumJumpsVec(n) = JumpFunc(Poi(n));  %JumpFunc(Poi(n)) sums up Poi(n) many jumps from the jump distribution
            end
        end

        W1 = randn(N_sim,1);
        Spath(:,m+1) = Spath(:,m).*exp(drift + sumJumpsVec + Sigsqdt*W1);  %log scheme
    end
end

end


text:
%io_loadspec_data.m
%Jamie Near, McGill University 2014.
%
% USAGE:
% [out,out_w]=io_loadspec_data(filename,sw,Larmor,subspecs,te,tr);
%
% DESCRIPTION:
% Reads in philips MRS data (.data and .list files) using code adapted from
% PhilipsRead_data.m, provided as part of the Gannet software package by
% Richard Edden (gabamrs.blogspot.com).
%
% op_loadspec_data outputs the data in structure format, with fields
% corresponding to time scale, fids, frequency scale, spectra, and header
% fields containing information about the acquisition.  The resulting
% matlab structure can be operated on by the other functions in this MRS
% toolbox.  NOTE:  Since the Gannet code is geared towards edited GABA MRS
% data, this code may not be general enough to handle all types of MRS data.
% Suggestions are most welcome.
%
% INPUTS:
% filename   = filename of Philips .data file to be loaded.
% sw         = spectral width (Hz)
% Larmor     = Larmor frequency (Hz/ppm, ie.  127 for 3T)
% subspecs   = number of subspectra in the data (from spectral editing, ISIS, etc.)
% te         = echo time (ms).  Optional, default is [].
% tr         = repetition time (ms).  Optional, default is [].
%
% OUTPUTS:
% out        = Input water suppressed dataset in FID-A structure format.
% out_w      = Input water reference dataset in FID-A structure format.

function [out,out_w]=io_loadspec_data(filename,sw,Larmor,subspecs,te,tr);

if nargin<6
    if nargin<5
        te=[];
    end
    tr=[];
end

%read in the data using the philipsDataLoad.m (adapted from PhilipsRead_data.m)
[FullData,WaterData]=philipsDataLoad(filename);

%As far as I can tell, the data that comes out of the philipsDataLoad
%function is normally a N x Navgs matrix.  The coils have already been
%combined. The Navgs dimension contains all the subspectra, so we will
%split them now.  Note, that in the data-list format that I have seen, the
%edit-OFF subspectra appear in elements [1 2 5 6 9 10 13 14...] and the
%edit-ON subspectra appear in the elements [3 4 7 8 11 12 15 16...].  Other sequences
%may result in a different subspecs order, but for now we will separate the
%subspectra in this way.
%If the data has multiple subspectra
if subspecs>1
    %First make an vector that holds the indices of the ON subspectra:
    totalAvgs=size(FullData,2);
    OFFindices=[1:2:totalAvgs]-mod([0:(totalAvgs/2)-1],2);
    ONindices=[2:2:totalAvgs]+mod([1:totalAvgs/2],2);
    %Now split the subspectra out of the "averages" dimension:
    data(:,:,1)=FullData(:,OFFindices);
    data(:,:,2)=FullData(:,ONindices);
else
    data=FullData;
end

fids=squeeze(data);
fids_w=squeeze(WaterData)';

sz=size(fids);
sz_w=size(fids_w);

%Find the magnetic field strength:
Bo=Larmor/42.577;

%Find the number of averages:
Naverages=size(fids,2)*size(fids,3);
Naverages_w=size(fids_w,2)*size(fids_w,3);

%In Philips data/list format, coil channels have already been combined:
Ncoils=1;
Ncoils_w=1;

%Now create a record of the dimensions of the data array.
dims.t=1;
dims.coils=0;
dims.averages=2;
if subspecs>1
    dims.subSpecs=3;
else
    dims.subSpecs=0;
end

dims_w.t=1;
dims_w.coils=0;
dims_w.averages=2;
dims_w.subSpecs=0;

specs=fftshift(ifft(fids,[],dims.t),dims.t);
specs_w=fftshift(ifft(fids_w,[],dims_w.t),dims_w.t);

%Now get relevant scan parameters:*****************************

%Get Spectral width and Dwell Time
spectralwidth=sw;
dwelltime=1/spectralwidth;

%Get TxFrq
txfrq=Larmor*1e6;

%Leave date blank
date='';

%Find the number of averages.  'averages' will specify the current number
%of averages in the dataset as it is processed, which may be subject to
%change.  'rawAverages' will specify the original number of acquired
%averages in the dataset, which is unchangeable.
%FOR WATER SUPPRESSED DATA:
if dims.subSpecs ~=0
    if dims.averages~=0
        averages=sz(dims.averages)*sz(dims.subSpecs);
        rawAverages=averages;
    else
        averages=sz(dims.subSpecs);
        rawAverages=1;
    end
else
    if dims.averages~=0
        averages=sz(dims.averages);
        rawAverages=averages;
    else
        averages=1;
        rawAverages=1;
    end
end

%FOR WATER UNSUPPRESSED DATA:
if dims_w.subSpecs ~=0
    if dims_w.averages~=0
        averages_w=sz(dims_w.averages)*sz(dims_w.subSpecs);
        rawAverages_w=averages_w;
    else
        averages_w=sz(dims_w.subSpecs);
        rawAverages_w=1;
    end
else
    if dims_w.averages~=0
        averages_w=sz(dims_w.averages);
        rawAverages_w=averages_w;
    else
        averages_w=1;
        rawAverages_w=1;
    end
end

%Find the number of subspecs.  'subspecs' will specify the current number
%of subspectra in the dataset as it is processed, which may be subject to
%change.  'rawSubspecs' will specify the original number of acquired
%subspectra in the dataset, which is unchangeable.
%FOR WATER SUPPRESSED DATA:
if dims.subSpecs ~=0
    subspecs=sz(dims.subSpecs);
    rawSubspecs=subspecs;
else
    subspecs=1;
    rawSubspecs=subspecs;
end

%FOR WATER UNSUPPRESSED DATA:
if dims_w.subSpecs ~=0
    subspecs_w=sz(dims.subSpecs);
    rawSubspecs_w=subspecs_w;
else
    subspecs_w=1;
    rawSubspecs_w=subspecs_w;
end

%****************************************************************

%Calculate t and ppm arrays using the calculated parameters:
f=[(-spectralwidth/2)+(spectralwidth/(2*sz(1))):spectralwidth/(sz(1)):(spectralwidth/2)-(spectralwidth/(2*sz(1)))];
ppm=f/(Bo*42.577);
ppm=ppm+4.65;

t=[0:dwelltime:(sz(1)-1)*dwelltime];

%FOR WATER SUPPRESSED DATA
%FILLING IN DATA STRUCTURE
out.fids=fids;
out.specs=specs;
out.sz=sz;
out.ppm=ppm;
out.t=t;
out.spectralwidth=spectralwidth;
out.dwelltime=dwelltime;
out.txfrq=txfrq;
out.date=date;
out.dims=dims;
out.Bo=Bo;
out.averages=averages;
out.rawAverages=rawAverages;
out.subspecs=subspecs;
out.rawSubspecs=rawSubspecs;
out.seq='';
out.te=te;
out.tr=tr;
out.pointsToLeftshift=0;

%FILLING IN THE FLAGS
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;
out.flags.averaged=0;
out.flags.addedrcvrs=0;
out.flags.subtracted=0;
out.flags.writtentotext=0;
out.flags.downsampled=0;
if out.dims.subSpecs==0
    out.flags.isFourSteps=0;
else
    out.flags.isFourSteps=(out.sz(out.dims.subSpecs)==4);
end

%FOR WATER UNSUPPRESSED DATA
%FILLING IN DATA STRUCTURE
out_w.fids=fids_w;
out_w.specs=specs_w;
out_w.sz=sz_w;
out_w.ppm=ppm;
out_w.t=t;
out_w.spectralwidth=spectralwidth;
out_w.dwelltime=dwelltime;
out_w.txfrq=txfrq;
out_w.date=date;
out_w.dims=dims_w;
out_w.Bo=Bo;
out_w.averages=averages_w;
out_w.rawAverages=rawAverages_w;
out_w.subspecs=subspecs_w;
out_w.rawSubspecs=rawSubspecs_w;
out_w.seq='';
out_w.te=te;
out_w.tr=tr;
out_w.pointsToLeftshift=0;

%FILLING IN THE FLAGS
out_w.flags.writtentostruct=1;
out_w.flags.gotparams=1;
out_w.flags.leftshifted=0;
out_w.flags.filtered=0;
out_w.flags.zeropadded=0;
out_w.flags.freqcorrected=0;
out_w.flags.phasecorrected=0;
out_w.flags.averaged=0;
out_w.flags.addedrcvrs=0;
out_w.flags.subtracted=0;
out_w.flags.writtentotext=0;
out_w.flags.downsampled=0;
if out_w.dims.subSpecs==0
    out_w.flags.isFourSteps=0;
else
    out_w.flags.isFourSteps=(out.sz(out.dims.subSpecs)==4);
end

%DONE

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function uvPixValidInd = sr_get_uvpix(img, opt)

% SR_GET_UVPIX
%
%

if(ndims(img) == 3)
    img = rgb2gray(img);
end

% Gradient
[imgGx, imgGy] = gradient(img);
imgGx2 = imgGx.^2;
imgGy2 = imgGy.^2;
imgGxy = imgGx.*imgGy;

% Blur
h = fspecial('gaussian',[5 1], 1.5);
imgGx2 = conv2(h, h', imgGx2, 'same');
imgGy2 = conv2(h, h', imgGy2, 'same');

imgGradEng = imgGx2 + imgGy2;
% figure(1); imagesc(imgGradEng > opt.gradThres);
imgGradEng = imgGradEng(opt.pRad+1:end-opt.pRad ,opt.pRad+1:end-opt.pRad);

uvPixValidInd = imgGradEng(:) > opt.gradThres;

end

repo: SelfExSR


text:
function[tau, Is] = unimodal(I, N)

if( nargin < 2)
  N = 100;
end

[I_hist, pos] = hist(I(:),N);

[vmax k_max] = max(I_hist);

pmax = pos(k_max);

vend = I_hist(end);
pend = pos(end);

m = (vend - vmax) / (pend - pmax);
alpha = pi/2 - atan(m);

k = pmax:(pend-pmax)/(N-k_max):pend;

dk = sqrt( (vmax-I_hist(k_max:end)).^2 + (pmax - k).^2 );
mk = (I_hist(k_max:end) - vmax) ./ (k-pmax);
alphak = pi/2 - atan(mk);

dpk = dk.*abs(sin(alpha-alphak));

%  figure; plot(dpk);

[d1 p1] = max(dpk);

tau = pmax + pos(p1);
Is = I >= tau;

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [sFileOut, errMsg] = out_fopen(RawFile, FileFormat, sFileIn, ChannelMat, iChannels)
% OUT_FOPEN: Saves the header of a new empty binary file.
%
% INPUTS:
%    - RawFile    : Full path of the file to create
%    - FileFormat : String ('EEG-EDF', 'BST-BIN', ...)
%    - sFileIn    : Structure of the header of the file to create
%    - ChannelMat : Channel file associated with the file to save
%    - iChannels  : Subset of channels from input file to save in output file

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2014-2019

% Parse inputs
if (nargin < 5) || isempty(iChannels)
    iChannels = [];
end
% Output variables
sFileOut = [];
errMsg = [];
% Get default epoch size
EpochSize = bst_process('GetDefaultEpochSize', sFileIn);

% Select channels in channel file
if ~isempty(iChannels)
    sFileIn.channelflag = sFileIn.channelflag(iChannels);
end

% Trying to open the output file
try
    % Create file header
    switch (FileFormat)
        case 'EEG-EGI-RAW'
            sFileOut = out_fopen_egi(RawFile, sFileIn, ChannelMat);
        case 'EEG-BRAINAMP'
            sFileOut = out_fopen_brainamp(RawFile, sFileIn, ChannelMat);
        case 'BST-BIN'
            sFileOut = out_fopen_bst(RawFile, sFileIn, ChannelMat, EpochSize);
        case 'SPM-DAT'
            sFileOut = out_fopen_spm(RawFile, sFileIn, ChannelMat);
        case 'EEG-EDF'
            sFileOut = out_fopen_edf(RawFile, sFileIn, ChannelMat, EpochSize, iChannels);
        case 'FIF'
            error('copy input file');
        case 'CTF-CONTINUOUS'
            error('copy input file');
        otherwise
            error('Unsupported file format');
    end
catch
    % Get error message
    e = lasterror();
    if ~isempty(e)
        errMsg = str_striptag(e.message);
    end
end

author: brainstorm-tools


text:
% phd_test.m
% ====================================================>
% This is a test script which demonstrates the usage of the "SMC_PHD" class.
%
% SETUP:
%  * Before running the simulation, open "2_crossing_targets.mat" or "3_roaming_targets.mat" datasets, from the "datasets" folder
%  * The datasets have been extracted by simulating the motion of differential robots in a 2D-plane (x,y)
%  * The "gen_obs_cluttered_multi3" function takes as an input the ground truth data, including information about the measurement noise and clutter rate
%     and then produces 1xNk cell array of corrupted and cluttered measurements, Nk being the total number of timesteps

% Load the ground truth data
load('multiple-robot-tracking.mat');

% Plot settings
ShowPlots = 1;              % Set to 0 to hide plots
ShowPrediction = 0;         % Set to 0 to skip showing prediction
ShowUpdate = 1;             % Set to 0 to skip showing update

lambdaV = 0; % Expected number of clutter measurements over entire surveillance region
V = 10^2;     % Volume of surveillance region (10x10 2D-grid)
V_bounds = [0 10 0 10]; % [x_min x_max y_min y_max]
P_D = 1;

% Instantiate a Transitionamic model
transition_model = ConstantVelocityX('NumDims',2,'VelocityErrVariance',0.0001);

% Instantiate a Measurement model
measurement_model = LinearGaussianX('NumMeasDims',2,'NumStateDims',4,'MeasurementErrVariance',0.02,'Mapping',[1 3]);
%measurement_model = RangeBearing2CartesianX('NumStateDims',4,'MeasurementErrVariance',[0.001,0.02],'Mapping',[1 3]);

% Instantiate a clutter model
clutter_model = PoissonRateUniformPositionX('ClutterRate',lambdaV,'Limits',[V_bounds(1:2);V_bounds(3:4)]);

% Instantiate birth model
birth_model = DistributionBasedBirthModelX('Distribution', UniformDistributionX([V_bounds(1:2); ...
                                                                                [-0.1 0.1 ];...
                                                                                V_bounds(3:4);...
                                                                                [-0.1 0.1 ]]),...
                                           'BirthIntensity', 0.0001);

% Instantiate detection model
detection_model = ConstantDetectionProbabilityX('DetectionProbability',P_D);

% Compile the State-Space model
model = StateSpaceModelX(transition_model,...
                       measurement_model,...
                       'Clutter',clutter_model,...
                       'Birth', birth_model,...
                       'Detection', detection_model);

% Extract the ground truth data from the example workspace
NumIter = size(GroundTruth,2);

% Set BirthIntensity
NumTracks = 3;

% Generate DataList
meas_simulator = MultiTargetMeasurementSimulatorX('Model',model);
DataList = meas_simulator.simulate(GroundTruthStateSequence);

% Assign PHD parameter values
config.Model = model;
[priorParticles, priorWeights] = model.Birth.random(50000);
config.StatePrior = ParticleStateX(priorParticles,10*priorWeights);
config.BirthScheme = {'Expansion', 5000};
config.SurvivalProbability = 0.99;

% Instantiate PHD filter
myphd = SMC_PHDFilterX(config);

% Create figure windows
if(ShowPlots)
    img = imread('maze.png');

    % set the range of the axes
    % The image will be stretched to this.
    min_x = 0;
    max_x = 10;
    min_y = 0;
    max_y = 10;

    % make data to plot - just a line.
    x = min_x:max_x;
    y = (6/8)*x;

    figure('units','normalized','outerposition',[0 0 .5 1])
    ax(1) = gca;
    figure('units','normalized','outerposition',[.5 0 .5 1])
    ax(2) = gca;
end

% START OF SIMULATION
% ===================>
for k=1:NumIter
    fprintf('Iteration = %d/%d\n================>\n',k,NumIter);

    % Extract DataList at time k
    tempDataList = DataList(k);

    % Change PHD filter parameters
    myphd.MeasurementList = tempDataList; % New observations

    % Predict PHD filter
    myphd.predict();

    % Update PHD filter
    myphd.update();

    fprintf("Estimated number of targets: %f\n", sum(myphd.StatePosterior.Weights));
    % Plot update step results
    if(ShowPlots && ShowUpdate)
        % Plot data
        cla(ax(1));
         % Flip the image upside down before showing it
        imagesc(ax(1),[min_x max_x], [min_y max_y], flipud(img));

        % NOTE: if your image is RGB, you should use flipdim(img, 1) instead of flipud.
        hold on;
        h2 = plot(ax(1), DataList(k).Vectors(1,:),DataList(k).Vectors(2,:),'k*','MarkerSize', 10);
        for j=1:NumTracks
            states = [GroundTruthTracks(j).Trajectory.Vector];
            h2 = plot(ax(1), states(1,1:k), states(3,1:k),'b.-','LineWidth',1);
            h2 = plot(ax(1), states(1,k), states(3,k),'bo','MarkerSize', 10);
        end
        p_i = myphd.IntensityPerHypothesis(2:end);

        ValidHypothesisInds = find(p_i>0.8)+1;
        for j = 1:length(ValidHypothesisInds)
            hypInd = ValidHypothesisInds(j);
            dist = ParticleDistributionX(myphd.StatePrediction.Particles,myphd.weightsPerHypothesis_(hypInd,:));
            dist.resample(1000);
            plot(ax(1), dist.Particles(1,:), dist.Particles(3,:), '.');
            text(ax(1), dist.Mean(1,:), dist.Mean(3,:), num2str(p_i(hypInd-1)),'FontSize',20,'Color','k');
        end

        % set the y-axis back to normal.
        set(ax(1),'ydir','normal');
        str = sprintf('Robot positions (Update)');
        title(ax(1),str)
        xlabel('X position (m)')
        ylabel('Y position (m)')
        axis(ax(1),V_bounds)

        % Plot PHD
        cla(ax(2), 'reset');
        [bandwidth,density,X,Y]=kde2d(myphd.StatePosterior.Particles([1,3],:)');
        %contour3(X,Y,density,50);
        h = surf(ax(2),X,Y,density);
        shading interp
        colormap(ax(2), jet(3000))
        %set(h, 'edgecolor','none')
        hold on;
%         plot(ax(2), myphd.Particles(1,:), myphd.Particles(3,:), '.')
%         hold on;
        plot(ax(2), myphd.MeasurementList.Vectors(1,:), myphd.MeasurementList.Vectors(2,:), 'y*');
        axis(ax(2), [V_bounds]);
        str = sprintf('PHD intensity (Update)');
        xlabel(ax(2),'X position (m)')
        ylabel(ax(2),'Y position (m)')
        zlabel(ax(2),'Intensity')
        title(ax(2),str)
        pause(0.01)
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function quality=meshquality(node,elem,maxnode)
%
% quality=meshquality(node,elem)
%
% compute the Joe-Liu mesh quality measure of an N-D mesh (N<=3)
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2011/02/26
%
% input:
%    node:  node coordinates of the mesh (nn x 3)
%    elem:  element table of an N-D mesh (ne x (N+1))
%
% output:
%    quality: a vector of the same length as size(elem,1), with
%           each element being the Joe-Liu mesh quality metric (0-1) of
%           the corresponding element. A value close to 1 represents
%           higher mesh quality (1 means equilateral tetrahedron);
%           a value close to 0 means nearly degenerated element.
%
% reference:
%    A. Liu, B. Joe, Relationship between tetrahedron shape measures,
%                    BIT 34 (2) (1994) 268-287.
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin<3)
    maxnode=4;
end
if(size(elem,2)>maxnode)
    elem=elem(:,1:maxnode);
end
enum=size(elem,1);
vol=elemvolume(node,elem);
edges=meshedge(elem);
ed=node(edges(:,1),:)-node(edges(:,2),:);
ed=sum((ed.*ed)');
ed=sum(reshape(ed,[enum length(ed)/enum])')';
dim=size(elem,2)-1;

coeff=10/9; % for tetrahedral
if(dim==2)
    coeff=1;
end
quality=coeff*dim*2^(2*(1-1./dim))*3^((dim-1)/2)*vol.^(2/dim)./ed;
maxquality=max(quality);
if(maxquality>1)
    quality=quality./maxquality;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function noise = orderedNoiseParamInit(noise, y)

% ORDEREDNOISEPARAMINIT ORDERED noise parameter initialisation.
% The ordered categorical noise model is an ordinal regression noise
% model. The real line is divided into categories, which have some
% ordering (such as small, medium, large).
%
% FORMAT
% DESC initialises the ordered categorical
%  noise structure with some default parameters.
% ARG noise : the noise structure which requires initialisation.
% RETURN noise : the noise structure with the default parameters placed in.
%
% SEEALSO : noiseCreate, noiseParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

if nargin > 1
  noise.C = max(max(y))+1;
  noise.numProcess = size(y, 2);
  noise.bias = zeros(1, noise.numProcess);
  for i = 1:size(y, 2)
    noise.bias(i) = mean(y(find(~isnan(y(:, i))), i));
  end
else
  noise.bias = repmat(1/2, 1, noise.numProcess);
end
noise.nParams = noise.C-2 + noise.numProcess;

if noise.C > 2
  noise.widths = repmat(1/(noise.C-2), noise.C-2, 1);
  noise.transforms.index = [noise.numProcess+1:noise.nParams];
  noise.transforms.type = optimiDefaultConstraint('positive');
else
  noise.widths = [];
end
noise.variance = 0.1; % needs to be set a bit above zero for numerical reasons.

% Can handle missing values?
noise.missing = 1;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function triangulation_test213 ( )

%*****************************************************************************80
%
%% TEST213 tests TRIANGULATION_ORDER3_QUAD.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 January 2007
%
%  Author:
%
%    John Burkardt
%
  quad_num = 6;

  quad_w = [ ...
    0.1666666666666666, ...
    0.1666666666666666, ...
    0.1666666666666666, ...
    0.1666666666666666, ...
    0.1666666666666666, ...
    0.1666666666666666 ];
  quad_xy = [ ...
    0.659027622374092,  0.231933368553031; ...
    0.659027622374092,  0.109039009072877; ...
    0.231933368553031,  0.659027622374092; ...
    0.231933368553031,  0.109039009072877; ...
    0.109039009072877,  0.659027622374092; ...
    0.109039009072877,  0.231933368553031  ]';
  test_num = 4;
  triangle_order = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST213\n' );
  fprintf ( 1, '  TRIANGULATION_ORDER3_QUAD can apply a quadrature rule\n' );
  fprintf ( 1, '  to every triangle in a triangulated region,\n' );
  fprintf ( 1, '  and estimate the integral of a function over\n' );
  fprintf ( 1, '  that region.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NODE_NUM   TRI_NUM  Integral estim  Area of Region\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
%
%  Set up the grid.
%
    n = 2^( test - 1 );
    node_num = ( n + 1 ) * ( n + 1 );

    node_xy = zeros(2,node_num);

    k = 0;
    for j = 1 : n + 1
      y = ( j - 1 ) / ( n + 1 - 1 );
      for i = 1 : n + 1
        x = ( i - 1 ) / ( n + 1 - 1 );
        k = k + 1;
        node_xy(1:2,k) = [ x, y ]';
      end
    end
%
%  Set up the triangulation.
%
    triangle_num = 2 * n * n;

    triangle_node = zeros(triangle_order,triangle_num);

    k = 0;
    for j = 1 : n
      for i = 1 : n

        n11 = i     + ( j     - 1 ) * ( n + 1 );
        n12 = i     + ( j + 1 - 1 ) * ( n + 1 );
        n21 = i + 1 + ( j     - 1 ) * ( n + 1 );
        n22 = i + 1 + ( j + 1 - 1 ) * ( n + 1 );

        k = k + 1;
        triangle_node(1:triangle_order,k) = [ n11, n21, n12 ]';
        k = k + 1;
        triangle_node(1:triangle_order,k) = [ n22, n12, n21 ]';

      end
    end
%
%  Estimate the integral.
%
    [ quad_value, region_area ] = triangulation_order3_quad ( node_num, ...
      node_xy, triangle_order, triangle_num, triangle_node, @quad_fun, ...
      quad_num, quad_xy, quad_w );

    fprintf ( 1, '  %8d  %8d  %14f  %14f\n', ...
      node_num, triangle_num, quad_value, region_area );
%
%  Delete arrays.
%
    node_xy = [];
    triangle_node = [];

  end

  return
end

author: johannesgerer

repo: jburkardt-m