text:
function main_image = PTKGetMainRegionExcludingPaddingBorder(original_image, threshold_image, minimum_region_volume_mm3, include_interior_regions, reporting)
    % PTKGetMainRegionExcludingPaddingBorder. Finds the largest connected region in a
    % binary 3D volume, excluding any regions which touch an inner border
    % created by extending an inner padding ring from the top and bottom of the image and growing horizontally.
    %
    % This is a heuristic that find the lung region while excluding regions exterior to the body even
    % if they are connected to the exterior through the airways or by the
    % lung extending beyond the image limits.
    %
    % Syntax:
    %     main_image = PTKGetMainRegionExcludingPaddingBorder(threshold_image, reporting)
    %
    % Inputs:
    %     threshold_image - a binary 3D volume as a PTKImage.
    %
    %     minimum_region_volume_mm3 (optional) - ignore any regions below this
    %         threshold
    %
    %     include_interior_regions - if true, the segmentation will include
    %         interior disconnected regions within the ROI even though they
    %         might not be part of the lungs and airways. Generally this
    %         might be set to true when finding a mask for airway
    %         segmentation, but set to false when finding a mask for lung
    %         segmentation
    %
    %     reporting (optional) - an object implementing the CoreReporting
    %         interface for reporting progress and warnings
    %
    % Outputs:
    %     main_image - a PTKImage binary volume of the found region
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    if ~isa(threshold_image, 'PTKImage')
        error('Requires a PTKImage as input');
    end

    if nargin < 2 || isempty(minimum_region_volume_mm3)
        minimum_region_volume_mm3 = 0;
    end

    if nargin < 3
        reporting = CoreReportingDefault;
    end

    reporting.ShowProgress('Finding lung region');

    if ~isempty(threshold_image.PaddingValue)
        padding_image = original_image.RawImage == threshold_image.PaddingValue;
    else
        padding_image = false(original_image.ImageSize);
    end
    padding_image(1, :, :) = true;
    padding_image(end, :, :) = true;
    padding_image(:, 1, :) = true;
    padding_image(:, end, :) = true;
    threshold_excluding_padding = threshold_image.RawImage > 0;
    threshold_excluding_padding(padding_image) = false;

    % Find an inner border of the padding at the top of the image
    top_slice_padding = padding_image(:, :, 1);
    top_slice_padding = imdilate(top_slice_padding, strel([0 1 0; 1 1 1; 0 1 0]), 'same');
    top_slice_connect_points = top_slice_padding & threshold_excluding_padding(:, :, 1);

    % Find an inner border of the padding at the bottom of the image
    bottom_slice_padding = padding_image(:, :, 1);
    bottom_slice_padding = imdilate(bottom_slice_padding, strel([0 1 0; 1 1 1; 0 1 0]));
    bottom_slice_connect_points = bottom_slice_padding & threshold_excluding_padding(:, :, end);

    vert_con = zeros([3,3,3]);
    vert_con(:,:,1) = [0 0 0; 0 1 0; 0 0 0];
    vert_con(:,:,2) = [0 0 0; 0 1 0; 0 0 0];
    vert_con(:,:,3) = [0 0 0; 0 1 0; 0 0 0];
    vert_cc = bwconncomp(threshold_excluding_padding, vert_con);
    lm_vert = labelmatrix(vert_cc);
    lm_vert_top = lm_vert(:, :, 1);
    lm_vert_bottom = lm_vert(:, :, end);
    top_border_indices = lm_vert_top(top_slice_connect_points(:));
    bottom_border_indices = lm_vert_bottom(bottom_slice_connect_points(:));
    vert_border_indices = unique([top_border_indices; bottom_border_indices]);

    % We extend each point of the inner border from the top and bottom of
    % the image vertically until they encounter non-air voxels. The
    % resulting set of voxels we define as border points.
    is_border_point = ismember(lm_vert, vert_border_indices);

    % Now get horizontal connected components
    hor_con = zeros([3,3,3]);
    hor_con(:,:,2) = [0 1 0; 1 1 1; 0 1 0];
    cc_hor = bwconncomp(threshold_excluding_padding, hor_con);
    lm_hor = labelmatrix(cc_hor);

    % Remove horizontal components that are connected to the border points
    % we determined earlier
    hor_border_indices = unique(lm_hor(is_border_point(:)));
    remove_points = ismember(lm_hor, hor_border_indices);

    threshold_excluding_padding(remove_points) = false;

    bordered_image = true(threshold_image.ImageSize + [2 2 2]);
    bordered_image(2:end-1, 2:end-1, 2:end-1) = threshold_excluding_padding;
    bordered_image(2:end-1, 2:end-1, 1) = false;
    bordered_image(2:end-1, 2:end-1, end) = false;

    image_opening_params_mm = [0, 2, 4, 6];
    image_opening_index = 1;
    still_searching = true;
    if max(threshold_image.VoxelSize) >= 2
        connectivity = 6;
    else
        connectivity = 26;
    end

    last_ball_element = -1;

    while still_searching

        opening_mm = image_opening_params_mm(image_opening_index);
        if opening_mm > 0
            ball_element = CoreImageUtilities.CreateBallStructuralElement(threshold_image.VoxelSize, opening_mm);
        else
            ball_element = [];
        end
        if isequal(ball_element, last_ball_element)
            % We skip the next opening if the structural element is
            % the same as last time (which could happen with large voxel
            % sizes)
            results = [];
        else
            if include_interior_regions
                max_num_regions = 10000;
            else
                max_num_regions = 2;
            end
            [results, CC] = OpenAndGetRegions(bordered_image, ball_element, connectivity, threshold_image.VoxelSize, minimum_region_volume_mm3, max_num_regions, ~include_interior_regions, reporting);
        end
        if numel(results) > 0 || image_opening_index == numel(image_opening_params_mm)
            still_searching = false;
        else
            image_opening_index = image_opening_index + 1;
            last_ball_element = ball_element;
        end
    end

    if numel(results) == 0
        reporting.Error('PTKGetMainRegionExcludingPaddingBorder:Failed', 'Unable to extract the main image region');
    end

    if ~include_interior_regions
        % Unless we include all interior regions, then determine whether to
        % combine the two results or only choose one of them
        if numel(results) > 1
            % If more than one region was found (after excluding boundary-touching
            % regions), then check if they are disconnected left and right lungs
            bounding_boxes = regionprops(CC, 'BoundingBox');
            bb_1 = bounding_boxes(results(1)).BoundingBox;
            bb_1 = bb_1([2,1,3,5,4,6]); % Bounding box is [xyz] but Matlab indices are [yxz]
            bb_2 = bounding_boxes(results(2)).BoundingBox;
            bb_2 = bb_2([2,1,3,5,4,6]); % Bounding box is [xyz] but Matlab indices are [yxz]
            image_centre = round(size(bordered_image)/2);

            image_centre_j = image_centre(2);
            centre_offset = round(image_centre_j/2);

            use_both_regions = false;
            if (bb_1(2) >= image_centre_j - centre_offset) && (bb_2(5) + bb_2(2) < image_centre_j + centre_offset)
                reporting.LogVerbose('I appear to have found 2 disconnected lungs. I am connecting them.');
                use_both_regions = true;
            end

            if (bb_2(2) >= image_centre_j - centre_offset) && (bb_1(5) + bb_1(2) < image_centre_j + centre_offset)
                reporting.LogVerbose('I appear to have found 2 disconnected lungs. I am connecting them.');
                use_both_regions = true;
            end
        else
            use_both_regions = false;
        end
    end

    bordered_image(:) = false;
    if include_interior_regions
        % Include ALL regions
        for index = 1 : numel(results)
            bordered_image(CC.PixelIdxList{results(index)}) = true;
        end
    else
        % Include the main or the main two regions
        bordered_image(CC.PixelIdxList{results(1)}) = true;
        if use_both_regions
            bordered_image(CC.PixelIdxList{results(2)}) = true;
        end
    end

    main_image = threshold_image.BlankCopy;
    main_image.ChangeRawImage(bordered_image(2:end-1, 2:end-1, 2:end-1));
end

function [results, CC] = OpenAndGetRegions(bordered_image_input, ball_element, connectivity, voxel_size, minimum_region_volume_mm3, max_num_regions, exclude_out_of_roi, reporting)

    bordered_image = bordered_image_input;
    if ~isempty(ball_element)
        bordered_image = imopen(bordered_image_input > 0, ball_element);
    end

    % Borders are marked as segmented - this ensures all components
    % touching the border are connected and allows us to eliminate them
    % when extracting the lung

    % Obtain connected component matrix
    CC = bwconncomp(bordered_image, connectivity);
    stats = regionprops(CC, 'BoundingBox');

    % Find largest region
    num_pixels = cellfun(@numel, CC.PixelIdxList);
    [sorted_largest_areas, sorted_largest_areas_indices] = sort(num_pixels, 'descend');

    % Remove regions that are below the volume threshold - NOTE we must
    % reduce the size due to possible gaps due to noise
    pixel_volume = prod(voxel_size);

    small_area_mask = pixel_volume*sorted_largest_areas >= minimum_region_volume_mm3/100;
    sorted_largest_areas_indices = sorted_largest_areas_indices(small_area_mask);
    sorted_largest_areas = sorted_largest_areas(small_area_mask);

    result_index = 1;
    results = [];
    index_in_sorted_array = 1;
    while (result_index <= max_num_regions && index_in_sorted_array <= length(sorted_largest_areas_indices))
        current_region_being_checked = sorted_largest_areas_indices(index_in_sorted_array);
        bordered_image(:) = false;
        bordered_image(CC.PixelIdxList{current_region_being_checked}) = true;
        if (bordered_image(1) || bordered_image(end))
            % This region is connected to the edge
%             reporting.ShowMessage('PTKGetMainRegionExcludingBorder:LargestROIConnectedToExterior', 'The largest region connected with the edge of the volume. I''m assuming this region is outside the body so choosing the next largest region');
        else
            bb = stats(current_region_being_checked).BoundingBox;
            if exclude_out_of_roi && (bb(2) > (5/9)*size(bordered_image_input, 2))
                disp('removed');
                % ToDo
            else
                results(result_index) = current_region_being_checked;
                result_index = result_index + 1;
            end
        end
        index_in_sorted_array = index_in_sorted_array + 1;
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function p = prior_loggaussian(varargin)
%PRIOR_LOGGAUSSIAN  Log-Gaussian prior structure
%
%  Description
%    P = PRIOR_LOGGAUSSIAN('PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    creates Log-Gaussian prior structure in which the named
%    parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    P = PRIOR_LOGGAUSSIAN(P, 'PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    modify a prior structure with the named parameters altered
%    with the specified values.
%
%    Parameters for Log-Gaussian prior [default]
%      mu       - location [0]
%      s2       - scale squared (variance) [1]
%      mu_prior - prior for mu [prior_fixed]
%      s2_prior - prior for s2 [prior_fixed]
%
%  See also
%    PRIOR_*
%
% Copyright (c) 2000-2001,2010 Aki Vehtari
% Copyright (c) 2010 Jaakko Riihimäki

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'PRIOR_LOGGAUSSIAN';
  ip.addOptional('p', [], @isstruct);
  ip.addParamValue('mu',0, @(x) isscalar(x));
  ip.addParamValue('mu_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('s2',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('s2_prior',[], @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'Log-Gaussian';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'Log-Gaussian')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('mu',ip.UsingDefaults)
    p.mu = ip.Results.mu;
  end
  if init || ~ismember('s2',ip.UsingDefaults)
    p.s2 = ip.Results.s2;
  end
  % Initialize prior structure
  if init
    p.p=[];
  end
  if init || ~ismember('mu_prior',ip.UsingDefaults)
    p.p.mu=ip.Results.mu_prior;
  end
  if init || ~ismember('s2_prior',ip.UsingDefaults)
    p.p.s2=ip.Results.s2_prior;
  end

  if init
    % set functions
    p.fh.pak = @prior_loggaussian_pak;
    p.fh.unpak = @prior_loggaussian_unpak;
    p.fh.lp = @prior_loggaussian_lp;
    p.fh.lpg = @prior_loggaussian_lpg;
    p.fh.recappend = @prior_loggaussian_recappend;
  end

end

function [w, s, h] = prior_loggaussian_pak(p)

  w=[];
  s={};
  h=[];
  if ~isempty(p.p.mu)
    w = p.mu;
    s=[s; 'Log-Gaussian.mu'];
    h = 1;
  end
  if ~isempty(p.p.s2)
    w = [w log(p.s2)];
    s=[s; 'log(Log-Gaussian.s2)'];
    h = [h 1];
  end
end

function [p, w] = prior_loggaussian_unpak(p, w)

  if ~isempty(p.p.mu)
    i1=1;
    p.mu = w(i1);
    w = w(i1+1:end);
  end
  if ~isempty(p.p.s2)
    i1=1;
    p.s2 = exp(w(i1));
    w = w(i1+1:end);
  end
end

function lp = prior_loggaussian_lp(x, p)

  lJ = -log(x);    % =log(1/x)=log(|J|) of transformation
  xt  = log(x);    % transformed x
  lp = 0.5*sum(-log(2*pi) -log(p.s2)- 1./p.s2 .* sum((xt-p.mu).^2,1)) +sum(lJ);

  if ~isempty(p.p.mu)
    lp = lp + p.p.mu.fh.lp(p.mu, p.p.mu);
  end
  if ~isempty(p.p.s2)
    lp = lp + p.p.s2.fh.lp(p.s2, p.p.s2) + log(p.s2);
  end
end

function lpg = prior_loggaussian_lpg(x, p)

  lJg = -1./x;     % gradient of log(|J|) of transformation
  xt   = log(x);   % transformed x
  xtg  = 1./x;     % derivative of transformation
  lpg = xtg.*(1./p.s2).*(p.mu-xt) + lJg;

  if ~isempty(p.p.mu)
    lpgmu = sum((1./p.s2).*(xt-p.mu)) + p.p.mu.fh.lpg(p.mu, p.p.mu);
    lpg = [lpg lpgmu];
  end
  if ~isempty(p.p.s2)
    lpgs2 = (sum(-0.5*(1./p.s2-1./p.s2.^2.*(xt-p.mu).^2 )) + p.p.s2.fh.lpg(p.s2, p.p.s2)).*p.s2 + 1;
    lpg = [lpg lpgs2];
  end
end

function rec = prior_loggaussian_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
  if ~isempty(p.p.mu)
    rec.mu(ri,:) = p.mu;
  end
  if ~isempty(p.p.s2)
    rec.s2(ri,:) = p.s2;
  end
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function ys = gaussian_likelihood(xs,Priors,Mu,Sigmas )
%GAUSSIAN_LIKELIHOOD Summary of this function goes here
%   Detailed explanation goes here
%
%   input ----------------------------------------------------------------
%
%       xs: (D x N)
%

K = size(Priors,2);
N = size(xs,2);
D = size(xs,1);

ys = zeros(1,N);

for i=1:N
    for k=1:K
        ys(i) = ys(i) + Priors(k).*gaussPDF(xs(:,i),Mu(:,k),Sigmas(:,:,k));
    end
end

end

author: epfl-lasa


text:
function [ complexdata ] = read_complex( fid, datasize, offset, ...
    datatype, complexbool, varargin )
%READ_COMPLEX Generic function for reading binary complex data.
%
%   Assumes i and q components interleaved (adjacent for a single pixel).
%   Uses regular fread functions to read in data.
%
%   complexdata = read_complex( fid, datasize, offset, datatype,...
%      complexbool, dim1range, dim2range, subsample, decimationfun )
%
%   INPUTS:
%      FID:           MATLAB file identifier.  Must refer to a file that is
%         open for reading.
%      DATASIZE:      1x2 array [number of elements in first dimension,
%         number of elements in the second dimension].  Total number of
%         elements in each dimension of the data, with dimension ordered as
%         there are written in the file.
%      OFFSET:        Index (in bytes) from the beginning of the file to
%         the beginning of the data.  Default is 0 (beginning of file).
%      DATATYPE:      MATLAB string for specifying binary data precision
%         ('uint8','float32',etc.) See FREAD documentation for full list.
%         Default is 'float32'.
%      COMPLEXBOOL:   Boolean whether data is complex (true) or real
%      DIM1RANGE:     1x2 array [first element, last element].  If reading
%         a subimage, this is the range of data to read, in the first
%         dimension (as written in the file).  Default is entire image
%         range.  Can also use an empty array ([]) to specify entire range.
%      DIM2RANGE:     1x2 array [first element, last element].   Default is
%         entire image range.  Can also use an empty array ([]) to specify
%         entire range.
%      SUBSAMPLE:     1x2 array [subsample rate in first dimension,
%         subsample rate in second dimension].  Default is [1 1] (no
%         subsampling).
%
%   OUTPUTS:
%      COMPLEXDATA:   Array of complex data values, of data type single
%         complex.  Is oriented so that the first dimension in the file
%         order matches the first dimension in MATLAB (vertical when
%         display in imshow, imagesc, etc.), so this may need to be
%         transposed for display for many image types.
%
% Author: Wade Schwartzkopf (NGA/IDT)
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if complexbool
    complexdata=read_bip(fid, datasize, offset, datatype, 2, true, varargin{:});
    complexdata=complex(complexdata(:,:,1:2:end),complexdata(:,:,2:2:end));
else
    complexdata=read_bip(fid, datasize, offset, datatype, 1, true, varargin{:});
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
classdef StrainDescriptor < FieldNameDescriptor

    methods (Access = protected)

        function loadFieldVariable(obj)
            obj.fieldName = 'strain';
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [R,theta] = round_sra(X,opt)
%% Given the moment matrix X, round R and theta from X
if nargin < 2
    opt = 1;
end
X     = X{1};
n     = size(X,1);
N     = round((n-10)/10);
[V,D] = eig(X);
[~,I] = sort(diag(D),'descend');
V     = V(:,I);
%% take the opt-th eigenvector
nropt = length(opt);
x     = V(:,opt);
R     = zeros(3,3,nropt);
theta = zeros(N,nropt);
for k = 1:nropt
    xk    = x(:,k);
    xk    = xk/xk(1);
    %% round that eigenvector
    Rk    = reshape( xk(1+blkIndices(1,9)),3,3 );
    Rk    = project2SO3(Rk);

    thetak= sign( xk(10+1:10+N) );
    thetak(thetak==0) = 1;

    R(:,:,k)     = Rk;
    theta(:,k)   = thetak;
end
end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function centroids = faceCentroids(nodes, faces)
%FACECENTROIDS Compute centroids of a mesh faces.
%
%   NORMALS = faceCentroids(VERTICES, FACES)
%   VERTICES is a set of 3D points  (as a N-by-3 array), and FACES is
%   either a N-by-3 index array or a cell array of indices. The function
%   computes the centroid of each face, and returns a Nf-by-3 array
%   containing their coordinates.
%
%   Example
%     [v e f] = createIcosahedron;
%     normals1 = faceNormal(v, f);
%     centros1 = faceCentroids(v, f);
%     figure; drawMesh(v, f);
%     hold on; axis equal; view(3);
%     drawVector3d(centros1, normals1);
%
%
%   See also:
%   meshes3d, drawMesh, faceNormal, convhull, convhulln
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2006-07-05
% Copyright 2006 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)

% deprecation warning
warning('geom3d:deprecated', ...
    [mfilename ' is deprecated, use ''meshFaceCentroids'' instead']);

if isnumeric(faces)
    % trimesh or quadmesh
    nf = size(faces, 1);
    centroids = zeros(nf, size(nodes, 2));
    if size(nodes, 2) == 2
        % planar case
        for f = 1:nf
            centroids(f,:) = polygonCentroid(nodes(faces(f,:), :));
        end
    else
        % 3D case
        for f = 1:nf
            centroids(f,:) = polygonCentroid3d(nodes(faces(f,:), :));
        end
    end
else
    % mesh with faces stored as cell array
    nf = length(faces);
    centroids = zeros(nf, size(nodes, 2));
    if size(nodes, 2) == 2
        % planar case
        for f = 1:nf
            centroids(f,:) = polygonCentroid(nodes(faces{f}, :));
        end
    else
        % 3D case
        for f = 1:nf
            centroids(f,:) = polygonCentroid3d(nodes(faces{f}, :));
        end
    end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ddydphi, d2gdx2, logL, dy, dy2, vy]= VBA_get_dL (gx, dG_dPhi, y, type, Qy, sigmaHat)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [ddydphi, d2gdx2, logL, dy, dy2, vy]= VBA_get_dL (gx, dG_dPhi, y, type, Qy, sigmaHat)
% Compute usefull intermediate values describing the misfit between a model
% prediction and an observation.
%
% IN:
%   - gx: model prediction about the observation y (1st order moment)
%   - dG_dPhi: derivative of gx wrt observation parameters
%   - y: actual observation
%   - type: flag defining the distribution of the observation (0: gaussian,
%     1: Bernouilli, 2: categorical)
%   - Qy: scaling matrix of the 2nd order moment of the prediction (only
%     gaussian observations)
%   - sigmaHat: scaling factor of the 2nd order moment of the prediction
%     (only gaussian observations)
%
% OUT:
%   - ddydphi: gradient of the prediction error wrt observation parameters
%   - d2gdx2: hessian of the prediction
%   - logL: log-likelihood of the observation given the prediction
%   - dy: prediction error
%   - dy2: normalized squared deviation
%   - vy: prediction variance
%
% /////////////////////////////////////////////////////////////////////////

if nargin<5
end

% prediction error
dy = y - gx;

switch type

    case 0 %--- normal
        vy=(1./sigmaHat).*diag(VBA_inv(Qy));
        ddydphi = sigmaHat.*(dG_dPhi*Qy*dy);
        d2gdx2 = sigmaHat.*(dG_dPhi*Qy*dG_dPhi');
        dy2=dy'*Qy*dy ;
        logL = - 0.5*sigmaHat.*dy2 ;
        logL = logL + 0.5*VBA_logDet(Qy*sigmaHat) - 0.5*numel(dy2)*log(2*pi) ;

    case 1 %--- binomial
        gx = VBA_finiteBinomial (gx);
        vy = gx.*(1-gx) ;
        ddydphi = dG_dPhi*(dy./vy);
        temp = y./(gx).^2 - (y-1)./(1-gx).^2;
        d2gdx2 = dG_dPhi*diag(temp)*dG_dPhi';
        logL = y'*log(gx) + (1-y)'*log(1-gx);
        dy2 = sum(temp);

   case 2   %--- multinomial
        gx = VBA_finiteBinomial (gx);
        vy = gx.*(1-gx) ;
        ddydphi = dG_dPhi*(y./gx);
        d2gdx2 = dG_dPhi*diag(y./gx.^2)*dG_dPhi';
        dy2 = sum(y./(gx).^2);
        logL = log(gx)'*y;

end

end

author: MBB-team

repo: VBA-toolbox


text:
function imgRec = HDRJPEG2000Dec(name)
%
%
%       imgRec = HDRJPEG2000Dec(name)
%
%
%       Input:
%           -name: the prefix of the compressed HDR images using JPEG HDR
%
%       Output:
%           -imgRec: the reconstructed HDR image
%
%     Copyright (C) 2011-2015  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

info = imfinfo(name); %read metadata
decoded = sscanf(cell2mat(info.Comments), '%g', 7);

nBit = decoded(end);
imgRec = double(imread(name)) / (2^nBit - 1);
for i=1:size(imgRec, 3)
    xMax = decoded((i - 1) * 2 + 1);
    xMin = decoded((i - 1) * 2 + 2);
    %range expansion
    imgRec(:,:,i) = exp(imgRec(:,:,i) * (xMax - xMin) + xMin) - 1e-6;
end

imgRec(imgRec < 0.0) = 0;

end

author: banterle


text:
function center = polylineCentroid(varargin)
%POLYLINECENTROID Compute centroid of a curve defined by a series of points
%
%   PT = polylineCentroid(POINTS);
%   Computes center of mass of a polyline defined by POINTS. POINTS is a
%   [NxD] array of double, representing a set of N points in a
%   D-dimensional space.
%
%   PT = polylineCentroid(PTX, PTY);
%   PT = polylineCentroid(PTX, PTY, PTZ);
%   Specifies points as separate column vectors
%
%   PT = polylineCentroid(..., TYPE);
%   Specifies if the last point is connected to the first one. TYPE can be
%   either 'closed' or 'open'.
%
%   Example
%   poly = [0 0;10 0;10 10;20 10];
%   polylineCentroid(poly)
%   ans =
%       [10 5]
%
%   See also:
%   polygons2d, centroid, polygonCentroid, polylineLength
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 22/05/2006.
%

%% process input arguments

% check whether the curve is closed
closed = false;
var = varargin{end};
if ischar(var)
    if strcmpi(var, 'closed')
        closed = true;
    end
    varargin = varargin(1:end-1);
end

% extract point coordinates
if length(varargin)==1
    points = varargin{1};
elseif length(varargin)==2
    points = [varargin{1} varargin{2}];
end

% compute centers and lengths composing the curve
if closed
    centers = (points + points([2:end 1],:))/2;
    lengths = sqrt(sum(diff(points([1:end 1],:)).^2, 2));
else
    centers = (points(1:end-1,:) + points(2:end,:))/2;
    lengths = sqrt(sum(diff(points).^2, 2));
end

% centroid of edge centers weighted by edge length
%weigths = repmat(lengths/sum(lengths), [1 size(points, 2)]);
center = sum(centers.*repmat(lengths, [1 size(points, 2)]), 1)/sum(lengths);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function S = inv(C)
% stiffness to compliance tensor
%
% Input
%  C - @complianceTensor
%
% Output
%  S - @complianceTensor
%

S = complianceTensor(inv@tensor(C));

author: mtex-toolbox

repo: mtex

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


text:
function C = slcovlarge(X, w, vmean, cachesize)
%SLCOVLARGE Computes large covariance matrix using memory-efficient way
%
% $ Syntax $
%   - C = slcovlarge(X, w, vmean, cachesize)
%
% $ Arguments $
%   - X:            the sample matrix
%   - w:            the weights of samples (default = [])
%                   if it is specifed, it should be a 1 x n vector
%   - vmean:        the pre-computed mean (default = [])
%   - cachesize:    the size of working cache (in the unit of Mbytes)
%
% $ Description $
%   - C = slcovlarge(X, w, vmean, cachesize) computes the large covariance
%     matrix within a memory-limited context. Compared to slcov, which is
%     faster at the expense of using more memory, it is slower, however
%     the memory used is under control. Thus it can effectively prevent
%     from the situation of out of memory.
%
% $ Remarks $
%   - The memory increased in the function will not exceed
%     the size of the covariance matrix plus the cachesize.
%
%   - For vmean
%       - if it is empty, then a mean vector will be calculated
%       - if it is zero, then the data is assumed to be centralized
%       - if it is specified, then its size should be d x 1.
%
%   - The memory estimated to be used for each section is
%     if has no weight
%           2 x sd x sn + 2 x (sd + sn).
%     else
%           3 x sd x sn + 2 x (sd + sn).
%     end
%     Here sd and sn are the sub-dimension
%     and sub-number of that section.
%
% $ History $
%   - Created by Dahua Lin, on Aug 17, 2006
%

%% parse and verify input arguments

if nargin < 4
    raise_lackinput('slcovlarge', 4);
end

if ndims(X) ~= 2 || ~isnumeric(X)
    error('sltoolbox:invalidarg', ...
        'X should be a 2D numeric matrix');
end

[d, n] = size(X);

if ~isempty(w)
    if ~isequal(size(w), [1, n])
        error('sltoolbox:sizmismatch', ...
            'The w should be an 1 x n vector');
    end
end

% convert cache size
cacheelems = floor(cachesize * 1e6 / 8);

if isempty(vmean)
    if cacheelems < d
        error('sltoolbox:notenoughmem', ...
            'The cache size is not enough to hold even a mean vector');
    end
    cacheelems = cacheelems - d;
    vmean = slmean(X, w);
else
    if ~isequal(vmean, 0)
        if ~isequal(size(vmean), [d, 1])
            error('sltoolbox:sizmismatch', ...
                'The vmean should be a d x 1 vector');
        end
    end
end

%% decide partition struct

if cacheelems > 2*d*n + 2*(d+n)   % can compute without partitioning
    ps = [];
else
    if isempty(w)
        cblks = 2;
    else
        cblks = 3;
    end

    sd = floor( (cacheelems - 2 * n) / (cblks * n + 2) );
    divide_row = false;
    if sd <= 0
        sd = 1;
        if cacheelems < cblks + 4
            error('sltoolbox:notenoughmem', ...
                'The cache is not large enough');
        end
        sn = floor((cacheelems - 2) / (cblks + 2));
        divide_row = true;
    end

    ps = slpartition(d, 'maxblksize', sd);
    if divide_row
        rps = slpartition(d, 'maxblksize', sn);
    end
end

%% compute

if isempty(ps)

    if ~isempty(w);
        X = weight_x(X, w);
    end
    X = shift_x(X, vmean);
    Xt = X';
    C = X * Xt;

else

    if ~divide_row      % each row is taken as integral

        C = zeros(d, d);
        nsecs = length(ps.sinds);

        for i = 1 : nsecs
            for j = 1 : nsecs
                si = ps.sinds(i);
                ei = ps.einds(i);
                sj = ps.sinds(j);
                ej = ps.einds(j);

                if isempty(w)
                    curXj = X(sj:ej, :);
                    curXj = shift_x(curXj, vmean, sj, ej);
                    curXjt = curXj';
                    clear curXj;

                    curXi = X(si:ei, :);
                    curXi = shift_x(curXi, vmean, si, ei);
                else
                    curXj = X(sj:ej, :);
                    curXj = shift_x(curXj, vmean, sj, ej);
                    curXj = weight_x(curXj, w);
                    clear curwj;
                    curXjt = curXj';
                    clear curXj;

                    curXi = X(si:ei, :);
                    curXi = shift_x(curXi, vmean, si, ei);
                    curXi = weight_x(curXi, w);
                    clear curwi;
                end

                C(si:ei, sj:ej) = curXi * curXjt;

            end
        end

    else                % even each row need to be divided

        slignorevars(rps);

        error('sltoolbox:rterror', ...
            'In current implementation, row-division is not implemeted yet');
    end

end

%% scale down

if isempty(w)
    C = C / n;
else
    tw = sum(w);
    C = C / tw;
end

%% Core computing function

function sx = weight_x(sx, w)

sn = size(sx, 2);
for i = 1 : sn
    sx(:, i) = sx(:, i) * sqrt(w(i));
end

function sy = shift_x(sx, vmean, sidx, eidx)

[sd, sn] = size(sx);
if ~isequal(vmean, 0)

    if nargin == 2
        curmean = vmean;
    else
        curmean = vmean(sidx:eidx);
    end

    sy = zeros(sd, sn);
    for i = 1 : sn
        sy(:,i) = sx(:,i) - curmean;
    end
else
    sy = sx;
end

author: lmthang


text:
function [errKN,errKI] = getErrInfIFE3D(uh, fun, mesh, femI, fem, dind)

%% USAGE: generate u-uh on all Gaussian points
%
% INPUTS:
% fun --- coefficient function
% mesh --- a struct data contains very rich mesh information.
% fem --- global DoF for test function space
% femI --- global DoF for test function space
% dind --- derivative info for test function
%            dind = [0,0,0]: function value
%            dind = [1,0,0]: Dx value
%            dind = [0,1,0]: Dy value
%            dind = [0,0,1]: Dz value
%
% OUTPUTS:
% matrix --- global (mass, stiffness, ...) matrix.

% Last Modified: 07/11/2020 by Xu Zhang
%% 0. Initializaiton
if strcmp(fem.type,'P1')||strcmp(fem.type,'DGP1')||strcmp(fem.type,'CR')
    feEvalBas = @evalP1Bas3D;
elseif strcmp(fem1.type,'P2')||strcmp(fem.type,'DGP2')
    feEvalBas = @evalP2Bas3D;
end
%% 1. Error on noninterface elements
dof = fem.ldof;
gx = fem.gx; gy = fem.gy; gz = fem.gz;
ntID = find(mesh.tLoc > 0);
gxN = gx(ntID,:); gyN = gy(ntID,:); gzN = gz(ntID,:);
uhK = uh(fem.t); uhKN = uhK(ntID,:);
tuN = feval(fun, gxN, gyN, gzN);
uKN = 0;
for i = 1:dof
    BAS = fem.bas(ntID,:,i);
    uKN = uKN + uhKN(:,i).*feEvalBas(BAS, gxN, gyN, gzN, dind);
end
errKN = norm(tuN-uKN,'inf');

%% 2. Error on interface elements
gxI = femI.gx; gyI = femI.gy; gzI = femI.gz;
itID = find(mesh.tLoc < 0); ntI = length(itID);
uhK = uh(fem.t); uhKitmp = uhK(itID,:); uhKI = zeros(length(femI.t),4);
id = 0;
for i = 1:ntI
    tmp = uhKitmp(i,femI.locID(i,:));
    uhKI(id+1:id+femI.quadT(i),:) = repmat(tmp,femI.quadT(i),1);
    id = id + femI.quadT(i);
end
tuI = feval(fun,gxI,gyI,gzI);
uKI = 0;
for i = 1:dof
    BAS = femI.bas(:,:,i);
    uKI = uKI + uhKI(:,i).*feEvalBas(BAS, gxI, gyI, gzI, dind);
end
errKI = norm(tuI-uKI,'inf');

repo: ifem


text:
t1 = t1/1000;
sl = 181;

figure;imagesc(ni.data(:,:,sl));axis image;colormap gray;colorbar;title('bob');
figure;imagesc(t1(:,:,sl));axis image;colormap gray;colorbar;title('deoni');

wm = roipoly;
[x,y] = ind2sub(ni.dim(1:2),find(wm));
wmInds = sub2ind(ni.dim(1:3),x,y,ones(size(x))*sl);
gm = roipoly;
[x,y] = ind2sub(ni.dim(1:2),find(gm));
gmInds = sub2ind(ni.dim(1:3),x,y,ones(size(x))*sl);

dm = [mean(t1(wmInds)) std(t1(wmInds)); mean(t1(gmInds)) std(t1(gmInds))];
km = [mean(ni.data(wmInds)),std(ni.data(wmInds)); mean(ni.data(gmInds)) std(ni.data(gmInds))];

(dm(2,1)-dm(1,1))/mean(dm(:,2))

(km(2,1)-km(1,1))/mean(km(:,2))

% Save a bunch of slices
sl = 115;

im = flipud(t1(:,:,sl)');
im = im./max(im(:));
imwrite(im,sprintf('images/deoni_T1_%03d.png',sl));

im = flipud(ni.data(:,:,sl)');
im = im./max(im(:));
imwrite(im,sprintf('images/bob_T1_%03d.png',sl));

allData = cat(4,s(:).imData);
[img, clipVals] = mrAnatHistogramClip(allData(:),.4,.99);
clear img allData;
for(ii=1:length(s))
   im = flipud(s(ii).imData(:,:,sl)');
   im(im>clipVals(2)) = clipVals(2);
   im(im<clipVals(1)) = clipVals(1);
   im = im-clipVals(1);
   im = im./max(im(:));
   if(ii==find(tiInds))
      imwrite(im,sprintf('images/raw_IR_%03d.png',sl));
   else
      imwrite(im,sprintf('images/raw_%02ddeg_%03d.png',s(ii).flipAngle,sl));
   end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% test for front propagation on 3D meshes
%
%   Copyright (c) 2007 Gabriel Peyre

path(path, 'toolbox/');
path(path, '../toolbox_graph_data/');
path(path, '../toolbox_graph_data/off/');

rep = ['results/geodesic-mesh/'];
if not(exist(rep))
    mkdir(rep);
end

disp('Loading mesh.');
if not(exist('name'))
    name = 'bunny';
    name = 'elephant-50kv';
    name = 'david50kf';
    name = 'david_head';
    name = 'hand';
end
options.name = name;
[vertex,faces] = read_mesh(name);
clf;
plot_mesh(vertex, faces,options);
saveas(gcf, [rep name '-mesh.png'], 'png');

nverts = max(size(vertex));
if not(exist('nstart'))
    nstart = 1;
end
start_points = [1 round(nverts/2)];
start_points = floor(rand(nstart,1)*nverts)+1;
start_points = start_points(:);
options.end_points = [];
% enforce a nice starting position for the first point
switch name
    case 'bunny'
        start_points(1) = 8900;
    case 'elephant-50kv'
        start_points(1) = 24575;
    case 'david50kf'
        start_points(1) = 20361;
    case 'david_head'
        start_points(1) = 18080;
    case 'hand'
        start_points(1) = 29719;
end

disp('Performing propagation.');
[D,S,Q] = perform_fast_marching_mesh(vertex, faces, start_points, options);

% compute geodesics
npaths = max(20,4*nstart);
% npaths = 1;
[tmp,I] = sort( D(:) ); I = I(end:-1:1); I = I(1:round(nverts*1));
end_points = floor( rand(npaths,1)*(length(I)-1) )+1;
end_points = I(end_points);
% [tmp,I] = sort( D(:) ); end_points(1) = I(end);

options.v2v = compute_vertex_ring(faces);
options.e2f = compute_edge_face_ring(faces);

disp('Extracting geodesics');
options.method = 'discrete';
options.method = 'continuous';
paths = compute_geodesic_mesh(D,vertex,faces, end_points, options);

options.colorfx = 'equalize';
plot_fast_marching_mesh(vertex,faces, D, paths, options);
saveas(gcf, [rep name '-geodesics-' num2str(nstart) '.png'], 'png');

author: alecjacobson

repo: gptoolbox


function [CVAR] = cvar_(y,lags,options)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'cvar_' estimates the paramters of VAR model with classical methods
% Impulse Response Function CI are computed using bootstrap methods

% Core Inputs:
% - y, data columns variables
% - lags, lag order of the VAR

% Additonal Inputs collected options:
% - options are not mandatory.
% See the Hitchhiker's guide for more details.
% https://github.com/naffe15/BVAR_/blob/master/HitchhikerGuide_.pdf

% Output:

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

if nargin < 2
    error('the cvar_ funtion needs at least two inputs: data and number of lags');
end
if lags < 1
    error('lags cannot be zero or negative');
end
% number of observable variables
ny                  = size(y, 2);

%********************************************************
%* DEFAULT SETTINGS
%********************************************************
% Control random number generator
if isOctave == 0
    isMatlab = 1;
    rng('default');
    rng(999);
else
    isMatlab = 0;
    % pkg load optim
    randn('state',999);
    rand('state',999);
end

% Default Settings (they can all be changed in 'options' see below)
K                   = 5000;         % number of draws from the posterior
hor                 = 24;           % horizon for the IRF
fhor                = 12;           % horizon for the forecasts
nethor              = 12;           % network horizon/ splioover and connectedness
firstobs            = lags+1;       % first observation
presample           = 0;            % using a presample for setting the hyper-parameter of the Minnesosta prior
noconstant          = 0;            % when 0, includes a constatn in the VAR
timetrend           = 0;            % when 1, includes a time trend in the VAR
long_run_irf        = 0;            % when 0, it does not compute long run IRF
irf_1STD            = 1;            % when 1, IRF are computed as 1SD increase. Else, IRF are compued as unitary increase in the shock
cfrcst_yes          = 0;            % no conditional forecast unless defined in options
% non_explosive_      = 0;            %
heterosked          = 0;
replacement         = 1;           % with replacement when bootstrapping
bootstrap           = 1;

signs_irf           = 0;
narrative_signs_irf = 0;
zeros_signs_irf     = 0;
proxy_irf           = 0;
heterosked_irf      = 0;
hmoments_signs_irf  = 0;
nexogenous          = 0;
exogenous           = [];
cnnctdnss_          = 0;
Ridge_              = 0;
Lasso_              = 0;
ElasticNet_         = 0;
set_irf             = 0;

if any(any(isnan(y))) %== 1
    error('cvar_ cannot handle missing observation; use bvar_ instead');

end

% declaring the names for the observable variables
for v = 1 : ny
    eval(['varnames{'   num2str(v) '} =  ''Var' num2str(v) ''';'])
end

%********************************************************
%* CUSTOMIZED SETTINGS
%********************************************************
if nargin > 2
    if isfield(options,'vnames')==1
        varnames = options.vnames;
    end
    %======================================================================
    % Inference options
    %======================================================================
    if isfield(options,'K')==1
        K = options.K;
    end
    if isfield(options,'firstobs')==1
        firstobs = options.firstobs;
        if firstobs < lags + 1
            error('firstobs need to be larger than lags +1')
        end
    end
    if isfield(options,'presample')==1
        presample = options.presample;
    end
    if isfield(options,'noconstant')==1
        noconstant = options.noconstant;
    end
    if isfield(options,'timetrend')==1
        timetrend  = options.timetrend;
        noconstant = 0;
    end
%     if isfield(options,'non_explosive_')==1
%         non_explosive_  = options.non_explosive_;
%     end
    if isfield(options,'heterosked_weights')==1
        ww  = options.heterosked_weights;
        heterosked = 1;
    end
    %======================================================================
    % Bootstrap options
    %======================================================================
    if isfield(options,'bootstrap')==1
        bootstrap = options.bootstrap;
    end
    if isfield(options,'replacement')==1
        replacement = options.replacement;
    end
    %======================================================================
    % Exogenous Variables options
    %======================================================================
    if isfield(options,'exogenous')==1 || isfield(options,'controls')==1
        if isfield(options,'controls')==1
            exogenous = options.controls;
        else
            exogenous = options.exogenous;
        end
        nexogenous = size(exogenous,2);
        if size(exogenous,1) ~= size(y,1)+ fhor && size(exogenous,1) ~= size(y,1)
            error('Size Mismatch between endogenous and exogenos variables; exo must be either T or T+fhor');
        end
        if any(isnan(exogenous(lags+1:end,:)))
            error('Exogenous variables cannnot be ''nan'' from lags+1 onward.');
        end
    end
    %======================================================================
    % IRF options
    %======================================================================
    if isfield(options,'hor') ==1
        hor = options.hor;
    end
    if isfield(options,'set_irf') ==1
        % set_irf is the # of rotations for set identified systems in wihch
        % Phi and Sigma are point estimates (e.g. OLS, Ridge, Lasso or
        % Elastic Nets); default 0
        set_irf = options.set_irf;
    end
    if isfield(options,'long_run_irf')==1
        % Activating Long run IRF
        long_run_irf = options.long_run_irf;
    end
    if isfield(options,'heterosked_regimes')==1
        % Activating identification via heteroskedasticity
        if length(options.heterosked_regimes) ~= size(y,1)
            error('heterosked_regimes must have the same time dimension of y')
        end
        heterosked_irf     = 1;
        heterosked_regimes = options.heterosked_regimes(lags+1:end);
        if any(heterosked_regimes>1) || any(heterosked_regimes<0)
            error('heterosked_regimes must contain zero (first regime) and one (second regime)')
        end
    end
    if isfield(options,'signs')==1
        % Activating IRF with sign restrictions (mulitple horizons allowed)
        signs_irf       = 1;
        signs           = options.signs;
        if iscellstr(signs) == 0
            error(['options.signs should be a cell array. Each cell must contain a string with the format'...
                '\n''y(a,b,c)<0'' or ''y(a,b,c)>0'' where a, b and c are integers.',...
                '\na = index of the variable',...
                '\nb = horizon',...
                '\nc = index of the shock'],class(zeros))
        end
    end
    if isfield(options,'narrative')==1
        if signs_irf  == 0
            warning('You did not provide any sign restrictions.')
            signs{1} = 'isempty(y(1,1,1))==0';
        end
        if signs_irf  == 1
            signs_irf       = 0;  % disactivating signs
        end
        narrative_signs_irf = 1;
        narrative           = options.narrative ;
    end
    if isfield(options,'hmoments')==1
        if signs_irf  == 0
            warning('You did not provide any sign restrictions.')
            signs{1} = 'isempty(y(1,1,1))==0';
        end
        if signs_irf  == 1
            signs_irf       = 0;  % disactivating signs
        end
        hmoments_signs_irf = 1;
        hmoments           = options.hmoments;
        [f]                = hmoments2matrix(hmoments,ny);
    end
    if isfield(options,'zeros_signs')==1
        if signs_irf  == 1
            signs_irf       = 0;  % disactivating signs
        end
        if narrative_signs_irf  == 1
            narrative_signs_irf       = 0;  % disactivating narrative
        end
        % Activating IRF with zeros and sign restrictions (mulitple horizons NOT allowed)
        zeros_signs_irf  = 1;
        zeros_signs      = options.zeros_signs;
        if iscellstr(zeros_signs) == 0
            error(['options.zeros_signs should be a cell array.'...
                '\nEach cell must contain a string with the following format'...
                '\nFor sign restrictions ''y(a,b)=1'' or ''y(a,b)=-1'',',...
                '\nFor short run zero restriction ''ys(a,b)=0'',',...
                '\nFor long run restriction ''yr(a,1,b)=0'' where a and b are integers.',...
                '\na = index of the variable',...
                '\nb = index of the shock'],class(zeros_signs))
        end
        [f,sr] = sign2matrix(zeros_signs,ny);
        if isfield(options,'var_pos')==0
            var_pos = ones(1,ny);
        else
            var_pos = options.var_pos;
        end
    end
    if isfield(options,'proxy')==1
        % Activating IRF with provy
        proxy_irf   = 1;
        in.proxies  = options.proxy;
        in.vars     = y;
        in.p        = lags;
        in.compute_F_stat = 0;
        if isfield(options,'proxy_end') == 1
            in.T_m_end  = options.proxy_end;
        else
            in.T_m_end  = 0;  %if the times series of the instrument ends when VAR data ends
        end
        in.irhor    = hor;
        if isnumeric(in.proxies) == 0
            error(['options.proxy should be a numeric array (nans or inf not allowed)'],class(in.proxies))
        end
        inols =in;
    end
    if isfield(options,'irf_1STD')==1
        % Activating of unitary IRF, i.e. a unitary increase in the shocks
        % (instead of 1 STD)
        irf_1STD = options.irf_1STD;
    end
    %======================================================================
    % (Un)Conditional Forecasts options
    %======================================================================
    if isfield(options,'fhor')==1
        fhor = options.fhor;
        if fhor < 1
            error('Forecast horizon must be positive')
        end
    end
    if isfield(options,'endo_index')==1
        % Forecast conditional on the path of an endogenous var
        cfrcst_yes      = 1;
        if isfield(options,'endo_path')== 0
            error('You need to provide the p[ath for the endogenou variable')
        end
        % rows forecasts, column variables
        endo_path       = options.endo_path;
        endo_path_index = options.endo_index;
        if length(endo_path_index) ~= size(endo_path)
            error(['Mismatch beween the number of endogenous paths and the number of conditioned variables'...
                '\nE.g. the # of conditioned variables must coincide with the # of column in ''options.endo_path'''],class(endo_path_index));
        end
        if isfield(options,'exo_index')==1
            % Forecast conditional on the path of an endo var using only a
            % subset of shocks. notice that the # of endo and # exo must coincide
            cfrcst_yes  = 2;
            exo_index   = options.exo_index;
            Omega       = eye(ny);
            %             if isfield(options,'Omegaf')==1
            %                 Omegaf = options.Omegaf;
            %             end
            if length(exo_index) ~= length(endo_path_index)
                error('the # of conditioned endogenous and # exogenous shocks used must coincide');
            end
        end

    end
    %======================================================================
    % Regularization options
    %======================================================================
    if isfield(options,'Ridge')==1
        Ridge_     = 1;%options.Ridge;
        if isfield(options.Ridge,'est') == 1
            Ridge_ = options.Ridge.est;
        end
        %cnnctdnss_ = 1;
        if isfield(options.Ridge,'lambda') == 1
            Ridge_lambda = options.Ridge.lambda;
        else
            warning('You did not specify a value for the penalization parameter (options.Ridge.lambda)')
            warning('I am using lambda = 0.02')
            Ridge_lambda = 0.02;
        end
    end
    if isfield(options,'Lasso')==1
        % (matlab stat toolbox needed)
        if exist('lasso') ~= 2
            error('Cannot estimate VAR with Lasso: matlab stat toolbox needed')
        end
        Lasso_     = 1;%options.Lasso;
        if isfield(options.Lasso,'est') == 1
            Lasso_ = options.Lasso.est;
        end
        %cnnctdnss_ = 1;
        if isfield(options.Lasso,'lambda') == 1
            Lasso_lambda = options.Lasso.lambda;
        else
            warning('You did not specify a value for the penalization parameter (options.Lasso.lambda)')
            %warning('Use the largest value of Lambda that gives a nonnull model')
            warning('I am using lambda = 0.05')
            Lasso_lambda = 0.05;
        end
    end
    if isfield(options,'ElasticNet')==1
        % (matlab stat toolbox needed)
        if  exist('lasso') ~= 2
            error('Cannot estimate VAR with ElasticNet: matlab stat toolbox needed')
        end
        ElasticNet_ = 1;%options.ElasticNet;
        if isfield(options.ElasticNet,'est') == 1
            ElasticNet_ = options.ElasticNet.est;
        end
        %cnnctdnss_  = 1;
        if isfield(options.ElasticNet,'lambda') == 1
            ElasticNet_lambda = options.ElasticNet.lambda;
        else
            warning('You did not specify a value for the penalization parameter (options.ElasticNet.lambda)')
            warning('I am using lambda = 0.05')
            ElasticNet_lambda = 0.05;
        end
        if isfield(options.ElasticNet,'alpha') == 1
            ElasticNet_alpha = options.ElasticNet.alpha;
        else
            warning('You did not specify a value for the relative penalization parameter (options.ElasticNet.alpha)')
            warning('I am using 0.5')
            ElasticNet_alpha = 0.5;
        end
    end
	%======================================================================
    % Network-Spillover-Connectedness options
    %======================================================================
    if isfield(options,'nethor')==1
        nethor = options.nethor;
    end
    if isfield(options,'connectedness')==1
        cnnctdnss_ = options.connectedness;
        if cnnctdnss_ == 2 && set_irf == 0
           set_irf = 100;
        end
    end
end

%********************************************************
%* Consistency Checks
%********************************************************
nobs = size(y,1)-firstobs+1;
if (firstobs+ nobs-1)> size(y,1)
    fprintf('Incorrect or missing specification of the number of observations. nobs can be at most %4u\n',size(y,1)-firstobs+1);
    error('Inconsistent number of observations.')
end
if firstobs + presample + lags >= nobs
    error('presample too large')
end
if firstobs + presample  <= lags
    error('firstobs+presample should be > # lags (for initializating the VAR)')
end
if cfrcst_yes ~= 0 && nexogenous > 0
    warning('I will not use exogenous variables with conditional forecasts');
    nexogenous = 0;
end
if size(exogenous,1) == size(y,1) && nexogenous > 0
    warning('For forecast purposes, I will assume that exo are zero out-of sample.')
    warning('To change this, include the exogenous forecasts in options.exogenous.')
    %fprintf('To change this, include the exogenous forecasts in options.exogenous.\n')
    exogenous = [exogenous; zeros(fhor,nexogenous)];
end
if (ElasticNet_==1 && Lasso_==1 && Ridge_==1) || ...
        (ElasticNet_==1 && Lasso_==1) || ...
        (Lasso_==1 && Ridge_==1) || ...
        (ElasticNet_==1 && Ridge_==1)
%    fprintf('Incorrect or missing specification of the number of observations. nobs can be at most %4u\n',size(y,1)-firstobs+1);
    error('You chose more than one regularization estimator, Ridge Lasso or ElasticNet. Please chose only one.')
end

%********************************************************
%* Estimate VAR
%********************************************************

idx = firstobs+presample-lags:firstobs+nobs-1;
nx  = 1;
if noconstant
    nx = 0;
end

% organize data as  yy = XX B + E
[yy,XX] = YXB_(y(idx, :),lags,[nx timetrend]);

% preparing the data for LS estimation
ydata   = y(idx, :);
T       = size(ydata, 1);
if T-lags < lags*ny + nx %+ flat*(ny+1)
    error('Less observations than regressors: increase the # of obs or decrease the # of lags.')
end
xdata   = ones(T,nx);
if timetrend == 1
    % xdata = [xdata [1:T]'];
    xdata = [xdata [1-lags : T-lags]'];
end
if nexogenous > 0
    xdata = [xdata exogenous(idx,:)];
    XX    = [XX exogenous(idx(1)+lags : idx(end),:)];
end

% OLS estimate:
if heterosked == 0
    varols  = rfvar3(ydata, lags, xdata, [T; T], 0, 0);
else
    varols  = rfvar3(ydata, lags, xdata, [T; T], 0, 0, ww);
end

if Ridge_ == 1
    varRidge   = varols;
    varRidge.B = (eye(size(varRidge.B,1)) + Ridge_lambda*varols.xxi)\varRidge.B;
    % varRidge.B = inv(eye(size(varRidge.B,1)) + Ridge_lambda*varols.xxi)*varRidge.B;
    varRidge.u = varRidge.y - varRidge.X*varRidge.B;
end
if Lasso_ == 1
    varLasso   = varols;
    for vv = 1 : ny
        varLasso.B(:,vv) = lasso(varLasso.X,varLasso.y(:,vv),'Lambda',Lasso_lambda);
    end
    varLasso.u = varLasso.y - varLasso.X*varLasso.B;
end
if ElasticNet_ == 1
     varElasticNet   = varols;
     for vv = 1 : ny
         varElasticNet.B(:,vv) = ...
             lasso(varElasticNet.X,varElasticNet.y(:,vv),'Lambda',ElasticNet_lambda,'Alpha',ElasticNet_alpha);
     end
     varElasticNet.u = varElasticNet.y - varElasticNet.X*varElasticNet.B;
end

%**************************************************
%* Generating draws form the Posterior Distribution
%**************************************************

% Preallocation of memory
% Matrices for collecting draws from Posterior Density
% Last dimension corresponds to a specific draw
Phi_boots     = zeros(ny*lags+nx+timetrend + nexogenous, ny, K);   % Autoregressive Parameters
Sigma_boots   = zeros(ny,ny,K);                     % Shocks Covariance
ir_boots      = zeros(ny,hor,ny,K);                 % variable, horizon, shock and draws - Cholesky IRF
irlr_boots    = zeros(ny,hor,ny,K);                 % variable, horizon, shock and draws - Long Run IRF
Qlr_boots     = zeros(ny,ny,K);                     % long run impact matrix
e_boots       = zeros(size(yy,1), ny,K);                  % residuals
yhatfut_no_shocks         = NaN(fhor, ny, K);   % forecasts with shocks
yhatfut_with_shocks       = NaN(fhor, ny, K);   % forecast without the shocks
yhatfut_cfrcst            = NaN(fhor, ny, K);   % forecast conditional on endogenous path
if signs_irf == 1
    irsign_boots = ir_boots;
    Omega_boots    = Sigma_boots;
end
if nexogenous > 0
    irx_boots = zeros(ny,hor,nexogenous,K);
    % Ox_boots  = Sigma_boots;
end
if narrative_signs_irf == 1
    irnarrsign_boots = ir_boots;
    Omegan_boots     = Sigma_boots;
end
if hmoments_signs_irf == 1
    irhmomsign_boots = ir_boots;
    Omegam_boots     = Sigma_boots;
end
if zeros_signs_irf == 1
    irzerosign_boots   = ir_boots;
    Omegaz_boots       = Sigma_boots;
end
if proxy_irf == 1
    irproxy_boots = ir_boots;
end
if heterosked_irf == 1
   irheterosked_boots = ir_boots;
   Omegah_boots       = Sigma_boots;
end
if cnnctdnss_ == 1
    CnndtnssIndex         = nan(K,1);
    CnndtnssFromAlltoUnit = nan(ny,K);
    CnndtnssFromUnitToAll = nan(ny,K);
    Ctheta                = nan(ny,ny,K);
end

% Settings for the forecasts
forecast_data.xdata       = ones(fhor, nx);
if timetrend
    forecast_data.xdata = [forecast_data.xdata (T-lags+1 : T-lags+fhor)'];
end
if nexogenous>0
    forecast_data.xdata = [forecast_data.xdata exogenous(T-lags+1 : T-lags+fhor,:)];
end
forecast_data.initval     = ydata(end-lags+1:end, :);

% # of RHS variables
nk               = ny*lags+nx+timetrend + nexogenous;

waitbar_yes = 0;
if K > 99
    waitbar_yes = 1;
    wb = waitbar(0, 'Bootstrapping');
end

% LS quantities for bootstrapping
Phi_    = varols.B;
e_      = varols.u;
Sigma_  = 1/(nobs-nk) * (e_' * e_);

% use penalized LS if activated in options
if Ridge_ == 1
    Phi_    = varRidge.B;
    e_      = varRidge.u;
    Sigma_  = 1/(nobs-nk) * (e_' * e_);
end
if Lasso_ == 1
    Phi_    = varLasso.B;
    e_      = varLasso.u;
    Sigma_  = 1/(nobs-nk) * (e_' * e_);
end
if ElasticNet_ == 1
    Phi_    = varElasticNet.B;
    e_      = varElasticNet.u;
    Sigma_  = 1/(nobs-nk) * (e_' * e_);
end

% setting of the forecast
deterministic_data.xdata = ones(nobs, nx);
if timetrend
    deterministic_data.xdata = [deterministic_data.xdata [1-lags : T-lags]'];
end
if nexogenous>0
    deterministic_data.xdata = [deterministic_data.xdata exogenous(idx,:)];
end

% Start bootstrap procedure
for  d =  1 : K
    % Bootstrap residuals
    switch bootstrap
        case 1
            % sampling with replacement
            % generating K new samples by bootstrapping the VAR reduced form errors
            e_boots(:,:,d) = bootstrap_(e_, 1, replacement);
            % generating K new initial conditions
            initval     =  bootstrap_(ydata(1 : firstobs-1,:), 1, replacement);
        case 2
            % Wild bootstrap based on simple distribution (~Rademacher)
            rr = 1-2*(rand(nobs,1)>0.5);
            %T x n randomly choose the sign of the time T shocks (all)
            e_boots(:,:,d) = (e_ .* (rr*ones(1,ny)));
            initval        = ydata(1 : firstobs-1,:);
    end
    % With bootstrapped errors / intial conditions construct new data
    deterministic_data.initval = initval;
    [~,ystar0]                 = forecasts(deterministic_data,Phi_,Sigma_,nobs,lags,e_boots(:,:,d),1);
    ystar                      = [initval; ystar0];
    % Construct left and right ahnd side variables for LS methods
    [yys,XXs] = YXB_(ystar(idx, :),lags,[nx timetrend]);
    % correct for heteroskedasticity of known form (ww) if any
    if heterosked == 1
        wws = [ww; ones(length(yys)-length(ww),1) ];
        yys = yys ./ repmat(wws,1,size(y,2)) ;
        XXs = XXs ./ repmat(wws,1,size(X,2));
    end
    % Compute LS estimator
    [vl,d_,vr] = svd(XXs,0);
    di         = 1./diag(d_);
    Phi        = (vr .* repmat(di', nk, 1)) * vl' * yys;
    % Penalized the LS estimator if activated in options
    if Ridge_ == 1
        xxi        = vr .* repmat(di',nk,1);
        iXX        = xxi * xxi';
        Phi = (eye(nk) + Ridge_lambda*iXX)\Phi;
    end
    if Lasso_ == 1
        LassoPhi = nan(size(Phi));
        for vv = 1 : ny
            LassoPhi(:,vv) = lasso(XXs,yys(:,vv),'Lambda',Lasso_lambda);
        end
        Phi = LassoPhi;
    end
    if ElasticNet_ == 1
        ElasticNetPhi = nan(size(Phi));
        for vv = 1 : ny
            ElasticNetPhi(:,vv) = ...
                lasso(XXs,yys(:,vv),'Lambda',ElasticNet_lambda,'Alpha',ElasticNet_alpha);
        end
        Phi = ElasticNetPhi;
    end
    % compute residual and covariance matrix
    errors  = yys - XXs*Phi;
    Sigma   = 1/(nobs-nk)*(errors' * errors);
    % store the draw
    Phi_boots(:,:,d)   = Phi;
    Sigma_boots(:,:,d) = Sigma;

    %======================================================================
    % IRF
    % Compute the impulse response functions
    % with cholesky
    if irf_1STD == 1
        % one STD increase
        ir_boots(:,:,:,d)      = iresponse(Phi(1 : ny*lags, 1 : ny),Sigma,hor,eye(ny));
    else
        % one percent increase
        ir_boots(:,:,:,d)      = iresponse(Phi(1 : ny*lags, 1 : ny),Sigma,hor,eye(ny),0);
    end
    if nexogenous > 0
        Phi1         = Phi(1 : ny*lags + nx + timetrend, 1 : ny);
        Exo1         = zeros(ny,nexogenous);
        Exo1(:,1:nexogenous) = Phi(ny*lags + nx + timetrend +1 ....
                            : ny*lags + nx + timetrend + nexogenous, 1 : ny)';
        Sig1         = eye(ny);
        % unitary increase
        irx_boots(:,:,:,d) = iresponse(Phi1,Sig1,hor,Exo1);
    end
    % define the identity rotation
    Omega = eye(ny);

    % with long run restrictions
    if long_run_irf == 1
        [irlr,Qlr]             = iresponse_longrun(Phi(1 : ny*lags, 1 : ny),Sigma,hor,lags);
        irlr_boots(:,:,:,d)    = irlr;
        Qlr_boots(:,:,d)       = Qlr;
        Omega                  = Qlr;
    end
%     % with sign restrictions
%     if signs_irf == 1
%         [irsign,Omega]         = iresponse_sign(Phi(1 : ny*lags, 1 : ny),Sigma,hor,signs);
%         irsign_boots(:,:,:,d)  = irsign;
%         Omega_boots(:,:,d)     = Omega;
%     end
%     % with narrative and sign restrictions
%     if narrative_signs_irf == 1
%         [irnarrsign,Omega]         = iresponse_sign_narrative(errors,Phi(1 : ny*lags, 1 : ny),Sigma,hor,signs,narrative);
%         irnarrsign_boots(:,:,:,d)  = irnarrsign;
%         Omegan_boots(:,:,d)        = Omega;
%     end
%     % with higher-moments and sign restrictions
%     if hmoments_signs_irf == 1
%         [irhmomsign,Omega]         = iresponse_sign_hmoments(errors,Phi(1 : ny*lags, 1 : ny),Sigma,hor,signs,hmoments,f);
%         irhmomsign_boots(:,:,:,d)  = irhmomsign;
%         Omegam_boots(:,:,d)        = Omega;
%     end
    % with zeros and sign restrictions
    if zeros_signs_irf == 1         %= iresponse_zeros_signs( Phi,Sigma,bvar1.hor,lags,var_pos,f,sr);
        [irzerosign,Omega]          = iresponse_zeros_signs(Phi,Sigma,hor,lags,var_pos,f,sr);
        irzerosign_boots(:,:,:,d)   = irzerosign;
        Omegaz_boots(:,:,d)         = Omega;
    end
    % with proxy
    if proxy_irf == 1
        in.res                  = e_boots(:,:,d);
        in.Phi                  = Phi_boots(:,:,d)  ;
        in.Sigma                = Sigma;
        tmp_                    = iresponse_proxy(in);
        irproxy_boots(:,:,1,d)  = tmp_.irs';
        clear tmp_
    end
    % with heteroskedasticity
    if heterosked_irf == 1
        [irheterosked,Omegah]       = iresponse_heterosked(Phi(1 : ny*lags, 1 : ny),errors,hor,heterosked_regimes);
        irheterosked_boots(:,:,:,d) = irheterosked;
        Omegah_boots(:,:,d)         = Omegah;
    end

    %======================================================================
    % Forecasts
    % compute the out of sample forecast (unconditional)
    [frcst_no_shock,frcsts_with_shocks] = forecasts(forecast_data,Phi,Sigma,fhor,lags);
    yhatfut_no_shocks(:,:,d)            = frcst_no_shock;
    yhatfut_with_shocks(:,:,d)          = frcsts_with_shocks;

    if cfrcst_yes == 1
        % Forecast conditional on the path of an endo var using all shocks
        [sims_with_endopath,EPS(:,:,d)] = ...
            cforecasts(endo_path,endo_path_index,forecast_data,Phi,Sigma);
        yhatfut_cfrcst(:,:,d) = sims_with_endopath;

    elseif cfrcst_yes == 2
        % Forecast conditional on the path of an endo var using only a
        % subset of shocks. notice that the # of endo and # exo must coincide
        % Omega is the structural orthonormal matrix
        [sims_with_endopath,EPS(:,:,d)] = ...
            cforecasts2(endo_path,endo_path_index,exo_index,forecast_data,Phi,Sigma,Omega);
        yhatfut_cfrcst(:,:,d) = sims_with_endopath;
    end

    %======================================================================
    % Connectedness
    if cnnctdnss_ > 0
%         if use_omega == 1
%             [C] = connectedness(Phi,Sigma,nethor,Omega);
%         else
%         end
        if cnnctdnss_ == 2
            Sigma_lower_chol = chol(Sigma)';
            [C] = connectedness(Phi(1 : ny*lags, 1 : ny),Sigma,nethor, Sigma_lower_chol * Omega );
        elseif cnnctdnss_ == 1
            [C] = connectedness(Phi(1 : ny*lags, 1 : ny),Sigma,nethor);
        end
        CnndtnssIndex(d,1)         = C.Index;
        CnndtnssFromAlltoUnit(:,d) = C.FromAllToUnit;
        CnndtnssFromUnitToAll(:,d) = C.FromUnitToAll;
        Ctheta(:,:,d)              = C.theta;
    end

    if waitbar_yes, waitbar(d/K, wb); end
%     dd = 0; % reset
end
if waitbar_yes, close(wb); end

%********************************************************
%* Storing the resutls
%*******************************************************

%==========================================================================
% classical inference: LS estimator
CVAR.Phi_ols    = varols.B;
CVAR.e_ols      = varols.u;
CVAR.Sigma_ols  = 1/(nobs-nk)*varols.u'*varols.u;
[CVAR.InfoCrit.AIC, CVAR.InfoCrit.HQIC, CVAR.InfoCrit.BIC] = IC(CVAR.Sigma_ols, CVAR.e_ols, nobs, nk);
% the model with the lowest IC is preferred

% OLS irf
% with cholesky
CVAR.ir_ols      = iresponse(CVAR.Phi_ols(1 : ny*lags, 1 : ny),CVAR.Sigma_ols,hor,eye(ny));
% with long run
if long_run_irf == 1
    [irlr,Qlr]              = iresponse_longrun(CVAR.Phi_ols(1 : ny*lags, 1 : ny),CVAR.Sigma_ols,hor,lags);
    CVAR.irlr_ols           = irlr;
    CVAR.Qlr_ols(:,:)       = Qlr;
end
% set identified IRF
if set_irf > 0
    wb = waitbar(0, ['Generating rotations for set-identification - LS Estimator']);
    % with sign restrictions
    if signs_irf == 1
        for d1 = 1 : set_irf
            [CVAR.irsign_ols(:,:,:,d1),CVAR.Omega_ols(:,:,d1)] = ...
                iresponse_sign(CVAR.Phi_ols(1 : ny*lags, 1 : ny),CVAR.Sigma_ols,hor,signs);
            waitbar(d1/set_irf, wb);
        end
    end
    % with narrative and sign restrictions
    if narrative_signs_irf == 1
        for d1 = 1 : set_irf
            [CVAR.irnarrsign_ols(:,:,:,d1),CVAR.Omegan_ols(:,:,d1)] = ...
                iresponse_sign_narrative(CVAR.e_ols,CVAR.Phi_ols(1 : ny*lags, 1 : ny),CVAR.Sigma_ols,hor,signs,narrative);
            waitbar(d1/set_irf, wb);
        end
    end
    % with zeros and sign restrictions
    if zeros_signs_irf == 1
        for d1 = 1 : set_irf
            [CVAR.irzerosign_ols(:,:,:,d1),CVAR.Omegaz_ols(:,:,d1)] = ...
                iresponse_zeros_signs(CVAR.Phi_ols,CVAR.Sigma_ols,hor,lags,var_pos,f,sr);
            waitbar(d1/set_irf, wb);
        end
    end
    close(wb)
end
% proxy
if proxy_irf == 1
    inols.res               = CVAR.e_ols;
    inols.Phi               = CVAR.Phi_ols;
    inols.Sigma             = CVAR.Sigma_ols;
    inols.compute_F_stat    = 1;
    tmp_                    = iresponse_proxy(inols);
    CVAR.irproxy_ols(:,:,1) = tmp_.irs';
    CVAR.proxy.F_m          = tmp_.F_m;
    CVAR.proxy.F_m_rob      = tmp_.F_m_rob;
    CVAR.proxy.R2adj_m      = tmp_.R2adj_m;
    CVAR.proxy.data         = options.proxy;
    clear tmp_
end
% with heteroskedasticity
if heterosked_irf == 1
    [CVAR.irheterosked_ols,CVAR.Omegah_ols] = ...
        iresponse_heterosked(CVAR.Phi_ols(1 : ny*lags, 1 : ny),CVAR.e_ols,hor,heterosked_regimes);
end
% test the normality of the ols VAR residuals (matlab stat toolbox needed)
if  exist('kstest') ==2
    for gg = 1 : ny
        [H,Pv] = kstest(CVAR.e_ols(:,gg)/sqrt(CVAR.Sigma_ols(gg,gg)));
        CVAR.HP(gg,:) = [H,Pv];
        %      H = 0 => Do not reject the null hypothesis at the 5% significance
        %      level.
    end
else
    CVAR.HP = [];
end

%==========================================================================
% Penalized Approaches (Regularization)
penalizationstrn = {'Ridge','Lasso','ElasticNet'};
% loop across penalization approaches
for pp = 1 : 3
    Phi_ = []; Sigma_ = []; u_ = []; Omega_(:,:,1) = eye(ny);  InfoCrit_ = [];
    if eval([ penalizationstrn{pp} '_ == 1'])
        % AR param
        eval(['CVAR.' penalizationstrn{pp} '.Phi    = var' penalizationstrn{pp} '.B;']);
        eval(['Phi_   = CVAR.' penalizationstrn{pp} '.Phi;'])
        % Error term
        eval(['CVAR.' penalizationstrn{pp} '.e      = var' penalizationstrn{pp} '.u;']);
        eval(['u_ = CVAR.' penalizationstrn{pp} '.e;'])
        % Covariance Matrix
        eval(['CVAR.' penalizationstrn{pp} '.Sigma  = 1/(nobs-nk) * var' penalizationstrn{pp} '.u'' * var' penalizationstrn{pp} '.u;'])
        eval(['Sigma_ = CVAR.' penalizationstrn{pp} '.Sigma;'])
        % Recursive IRFs
        eval(['CVAR.' penalizationstrn{pp} '.ir      = iresponse(Phi_(1 : ny*lags, 1 : ny),Sigma_,hor,eye(ny));']);
        % info crit
        [InfoCrit_.AIC, InfoCrit_.HQIC, InfoCrit_.BIC] ...
            = IC(Sigma_, u_, nobs, nk);
        eval(['CVAR.' penalizationstrn{pp} '.InfoCrit      = InfoCrit_;']);
        % IRFs with different identification schemes
        % point identification: LR
        if long_run_irf == 1
            eval(['[CVAR.' penalizationstrn{pp} '.irlr,Omega_] = iresponse_longrun(Phi_(1 : ny*lags, 1 : ny),Sigma_,hor,lags);'])
        end
        % point identification: with heteroskedasticity
        if heterosked_irf == 1
            eval(['[CVAR.' penalizationstrn{pp} '.irheterosked,Omega_] = iresponse_heterosked(Phi_(1 : ny*lags, 1 : ny),u_,hor,heterosked_regimes);'])
        end
        % set identification: signs
        if set_irf > 0
            wb = waitbar(0, ['Generating rotations for set-identification - ' penalizationstrn{pp} ' Estimator']);
        end
        if signs_irf == 1
            for d1 = 1 : set_irf
                eval(['[CVAR.' penalizationstrn{pp} '.irsign_boots(:,:,:,' num2str(d1) ...
                    '),Omega_(:,:,' num2str(d1) ')] = iresponse_sign(Phi_(1 : ny*lags, 1 : ny),Sigma_,hor,signs);'])
                waitbar(d1/set_irf, wb);
            end
        end
        % set identification: narrative and sign restrictions
        if narrative_signs_irf == 1
            for d1 = 1 : set_irf
                eval(['[CVAR.' penalizationstrn{pp} '.irnarrsign_boots(:,:,:,' num2str(d1) ...
                    '),Omega_(:,:,' num2str(d1) ')] = iresponse_sign_narrative(Phi_(1 : ny*lags, 1 : ny),Sigma_,hor,signs,narrative);'])
                waitbar(d1/set_irf, wb);
            end
        end
        % set identification: zeros and sign restrictions
        if zeros_signs_irf == 1
            for d1 = 1 : set_irf
                eval(['[CVAR.' penalizationstrn{pp} 'irzerosign_boots.(:,:,:,' num2str(d1) ...
                    '),Omega_(:,:,' num2str(d1) ')] = iresponse_zeros_signs(Phi_(1 : ny*lags, 1 : ny),Sigma_,hor,lags,var_pos,f,sr);'])
                waitbar(d1/set_irf, wb);
            end
        end
        if set_irf>0, close(wb); end
    end
    if cnnctdnss_ == 1 % default identification (Pesaran and Shin)
        eval(['CVAR.' penalizationstrn{pp} '.Connectedness = connectedness(Phi_(1 : ny*lags, 1 : ny),Sigma_,nethor);']);
    elseif cnnctdnss_ == 2 % customized identification
        Sigma_lower_chol = chol(Sigma_)';
        % this part is not correct as the median orthogonal rotation is not
        % orthogonal
        Omegam           = median(Omega_,3);
        eval(['CVAR.' penalizationstrn{pp} '.Connectedness = connectedness(Phi_(1 : ny*lags, 1 : ny), Sigma_, nethor, Sigma_lower_chol * Omegam);']);
    end

end

%==========================================================================
% Store the output
% the last dimension of these objects corresponds to a bootstrapped sample

% inference and IRFs
CVAR.Phi_boots    = Phi_boots;          % draws from the autoregressive part
CVAR.Sigma_boots  = Sigma_boots;        % draws from the covarance matrix
CVAR.ir_boots     = ir_boots;           % draws from the IRF with cholesky
CVAR.irlr_boots   = irlr_boots;         % draws from the IRF with Long Run
CVAR.Qlr_boots    = Qlr_boots;          % Long Run Rotation matrix
CVAR.lags         = lags;               % lags
CVAR.N            = ny;                 % number of variables
CVAR.e_boots      = e_boots;            % residuals
CVAR.XX           = XX;                 % regressors (no dummy)
CVAR.yy           = yy;                 % dependent  (no dummy)

% prediction
CVAR.fhor         = fhor;               % forecast horizon
CVAR.hor          = hor;                % IRF horizon
CVAR.forecasts.no_shocks      = yhatfut_no_shocks;         % trajectories of forecasts without shocks
CVAR.forecasts.with_shocks    = yhatfut_with_shocks;       % trajectories of forecasts with shocks
CVAR.forecasts.conditional    = [];            % trajectories of conditional forecasts
CVAR.forecasts.EPScond        = [];            % shocks of conditional forecasts
if cfrcst_yes ~= 0
    CVAR.forecasts.conditional    = yhatfut_cfrcst;            % trajectories of forecasts
    CVAR.forecasts.EPScond        = EPS;                       % shocks of forecasts
end
CVAR.forecast_data            = forecast_data;

%
CVAR.varnames     = varnames;
CVAR.ndraws       = K;

if signs_irf == 1 && narrative_signs_irf == 0
    CVAR.irsign_boots = irsign_boots;
    CVAR.Omegas       = Omega_boots;
else
    CVAR.irsign_boots = [];
    CVAR.Omegas       = [];
end
if narrative_signs_irf == 1
    CVAR.irnarrsign_boots = irnarrsign_boots;
    CVAR.Omegan           = Omegan_boots;
else
    CVAR.irnarrsign_boots = [];
    CVAR.Omegan           = [];
end
if hmoments_signs_irf == 1
    CVAR.irhmomsign_boots = irhmomsign_boots;
    CVAR.Omegam           = Omegam_boots;
else
    CVAR.irhmomsign_boots = [];
    CVAR.Omegam           = [];
end
if zeros_signs_irf == 1
    CVAR.irzerosign_boots   = irzerosign_boots;
    CVAR.Omegaz             = Omegaz_boots;
else
    CVAR.irzerosign_boots = [];
    CVAR.Omegaz           = [];
end
if heterosked_irf == 1
    CVAR.irheterosked_boots   = irheterosked_boots;
    CVAR.Omegah_boots         = Omegah_boots;
    CVAR.Omegah               = CVAR.Omegah_boots;
else
    CVAR.irheterosked_boots   = [];
    CVAR.Omegah_boots         = [];
    CVAR.Omegah               = [];
end
if proxy_irf == 1
    CVAR.irproxy_boots = irproxy_boots;
else
    CVAR.irproxy_boots= [];
end
if nexogenous > 0
    CVAR.irx_boots = irx_boots;
else
    CVAR.irx_boots = [];
end

% connetedness computed using LS (or penalized LS if activated in options)
if cnnctdnss_
    CVAR.Connectedness.Index         = CnndtnssIndex;
    CVAR.Connectedness.FromAlltoUnit = CnndtnssFromAlltoUnit;
    CVAR.Connectedness.FromUnitToAll = CnndtnssFromUnitToAll;
    CVAR.Connectedness.theta         = Ctheta;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of cvar_.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function [event] = read_yokogawa_event(filename, varargin)

% READ_YOKOGAWA_EVENT reads event information from continuous,
% epoched or averaged MEG data that has been generated by the Yokogawa
% MEG system and software and allows those events to be used in
% combination with FieldTrip.
%
% Use as
%   [event] = read_yokogawa_event(filename)
%
% See also READ_YOKOGAWA_HEADER, READ_YOKOGAWA_DATA

% Copyright (C) 2005, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

event   = [];
handles = definehandles;

% get the options, the default is set below
chanindx    = ft_getopt(varargin, 'chanindx');
threshold   = ft_getopt(varargin, 'threshold');
detectflank = ft_getopt(varargin, 'detectflank');

% ensure that the required toolbox is on the path
if ft_hastoolbox('yokogawa_meg_reader')
  % read the dataset header
  hdr = read_yokogawa_header_new(filename);
  ch_info = hdr.orig.channel_info.channel;
  type = [ch_info.type];

  % determine the trigger channels (if not specified by the user)
  if isempty(chanindx)
    chanindx = find(type==handles.TriggerChannel);
  end

  % Use the MEG Reader documentation if more detailed support is required.
  if hdr.orig.acq_type==handles.AcqTypeEvokedRaw
    % read the trigger id from all trials
    event = getYkgwHdrEvent(filename);
    % use the standard FieldTrip header for trial events
    % make an event for each trial as defined in the header
    for i=1:hdr.nTrials
      event(end+1).type     = 'trial';
      event(end  ).sample   = (i-1)*hdr.nSamples + 1;
      event(end  ).offset   = -hdr.nSamplesPre;
      event(end  ).duration =  hdr.nSamples;
      if ~isempty(value)
        event(end  ).value    =  event(i).code;
      end
    end

  % Use the MEG Reader documentation if more detailed support is required.
  elseif hdr.orig.acq_type==handles.AcqTypeEvokedAve
    % make an event for the average
    event(1).type     = 'average';
    event(1).sample   = 1;
    event(1).offset   = -hdr.nSamplesPre;
    event(1).duration =  hdr.nSamples;

  elseif hdr.orig.acq_type==handles.AcqTypeContinuousRaw
    % Events annotated by users during measurements
    bookmark_tmp = [];
    bookmark_tmp = getYkgwHdrBookmark(filename);
    if ~isempty(bookmark_tmp)
      for i = 1:length(bookmark_tmp)
        event(end+1).sample = bookmark_tmp(i).sample_no + 1;
        event(end  ).type   =  'bookmarks';
        if ~isempty(bookmark_tmp(i).label)
          event(end  ).value = bookmark_tmp(i).label;
          %% 0 in default
        else
          event(end  ).value = 99;
        end
      end
    end
    clear bookmark_tmp;
  end

elseif ft_hastoolbox('yokogawa');

  % read the dataset header
  hdr = read_yokogawa_header(filename);

  % determine the trigger channels (if not specified by the user)
  if isempty(chanindx)
    chanindx = find(hdr.orig.channel_info(:,2)==handles.TriggerChannel);
  end

  if hdr.orig.acq_type==handles.AcqTypeEvokedRaw
    % read the trigger id from all trials
    fid   = fopen_or_error(filename, 'r');
    value = GetMeg160TriggerEventM(fid);
    fclose(fid);
    % use the standard FieldTrip header for trial events
    % make an event for each trial as defined in the header
    for i=1:hdr.nTrials
      event(end+1).type     = 'trial';
      event(end  ).sample   = (i-1)*hdr.nSamples + 1;
      event(end  ).offset   = -hdr.nSamplesPre;
      event(end  ).duration =  hdr.nSamples;
      if ~isempty(value)
        event(end  ).value    =  value(i);
      end
    end

  elseif hdr.orig.acq_type==handles.AcqTypeEvokedAve
    % make an event for the average
    event(1).type     = 'average';
    event(1).sample   = 1;
    event(1).offset   = -hdr.nSamplesPre;
    event(1).duration =  hdr.nSamples;

  elseif hdr.orig.acq_type==handles.AcqTypeContinuousRaw
    % the data structure does not contain events, but flank detection on the trigger channel might reveal them
    % this is done below
  end

else
    ft_error('cannot determine, whether Yokogawa toolbox is present');
end

% read the trigger channels and detect the flanks
if ~isempty(chanindx)
  trigger = read_trigger(filename, 'header', hdr, 'denoise', false, 'chanindx', chanindx, 'detectflank', detectflank, 'threshold', threshold);
  % combine the triggers and the other events
  event = appendstruct(event, trigger);
end

if isempty(event)
  ft_warning('no triggers were detected, please specify the "trigindx" option');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this defines some usefull constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handles = definehandles
handles.output = [];
handles.sqd_load_flag = false;
handles.mri_load_flag = false;
handles.NullChannel         = 0;
handles.MagnetoMeter        = 1;
handles.AxialGradioMeter    = 2;
handles.PlannerGradioMeter  = 3;
handles.RefferenceChannelMark = hex2dec('0100');
handles.RefferenceMagnetoMeter       = bitor( handles.RefferenceChannelMark, handles.MagnetoMeter );
handles.RefferenceAxialGradioMeter   = bitor( handles.RefferenceChannelMark, handles.AxialGradioMeter );
handles.RefferencePlannerGradioMeter = bitor( handles.RefferenceChannelMark, handles.PlannerGradioMeter );
handles.TriggerChannel      = -1;
handles.EegChannel          = -2;
handles.EcgChannel          = -3;
handles.EtcChannel          = -4;
handles.NonMegChannelNameLength = 32;
handles.DefaultMagnetometerSize       = (4.0/1000.0);       % Square of 4.0mm in length
handles.DefaultAxialGradioMeterSize   = (15.5/1000.0);      % Circle of 15.5mm in diameter
handles.DefaultPlannerGradioMeterSize = (12.0/1000.0);      % Square of 12.0mm in length
handles.AcqTypeContinuousRaw = 1;
handles.AcqTypeEvokedAve     = 2;
handles.AcqTypeEvokedRaw     = 3;
handles.sqd = [];
handles.sqd.selected_start  = [];
handles.sqd.selected_end    = [];
handles.sqd.axialgradiometer_ch_no      = [];
handles.sqd.axialgradiometer_ch_info    = [];
handles.sqd.axialgradiometer_data       = [];
handles.sqd.plannergradiometer_ch_no    = [];
handles.sqd.plannergradiometer_ch_info  = [];
handles.sqd.plannergradiometer_data     = [];
handles.sqd.eegchannel_ch_no   = [];
handles.sqd.eegchannel_data    = [];
handles.sqd.nullchannel_ch_no   = [];
handles.sqd.nullchannel_data    = [];
handles.sqd.selected_time       = [];
handles.sqd.sample_rate         = [];
handles.sqd.sample_count        = [];
handles.sqd.pretrigger_length   = [];
handles.sqd.matching_info   = [];
handles.sqd.source_info     = [];
handles.sqd.mri_info        = [];
handles.mri                 = [];

author: spm


text:
function MISE=GaussKernelGaussMixMISE(n,H,w,mu,P)
%%GAUSSKERNELGAUSSMIXMISE Given a Gaussian mixture probability density
%            function (PDF) that one is attempting to approximate using a
%            kernel density estimator with Gaussian kernels all having
%            square root covariance matrix H, determine the mean integrated
%            squared error (MISE) of the estimate. The MISE is the expected
%            value of int_{R}(\hat{f}(x)-f(x))^2dx where \hat{f}(x) is the
%            kernel estimator based on a set of random samples of the true
%            density f(x) (since it is based on random samples, the mean
%            part is the expected value being taken over all possible
%            samples). This function is useful for assessing different
%            kernel bandwidth estimators (estimators of H) when using
%            Gaussian mixture sample problems.
%
%INPUTS: n The number of samples that the kernel-density estimator will
%          use.
%        H The dXd bandwidth matrix used in the kernel density estimator
%          (as one would use in the function kernelApprox).
%        w A numCompX1 or 1XnumComp vector of the weights associated with
%          each of the numComp components of the Gaussian mixture.
%       mu The dXnumComp matrix of the mean of each component of the
%          Gaussian mixture.
%        P The dXdXnumComp set of covariance matrices of the components of
%          the Gaussian mixture.
%
%OUTPUTS: MISE The MISE of thekernel estimator with the given number of
%              samples and bandwidth matrix.
%
%This function implements the formula for the MISE in an arbitrary number
%of dimensions in Section 4 of [1].
%
%REFERENCES:
%[1] M. P. Wand and M. C. Jones, "Comparison of smoothing parameterizations
%    in bivariate kernel density estimation," Journal of the American
%    Statistical Association, vol. 88, no. 422, pp. 520-528, Jun. 1993.
%
%December 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The number of components in the
k=size(mu,2);
d=size(mu,1);

%The definition of the bandwidth used in the paper is the square of what is
%used here.
H=H*H';

Omega0=zeros(k,k);
Omega1=zeros(k,k);
Omega2=zeros(k,k);

for l=1:k
    for lp=l:k
       Sigma=P(:,:,l)+P(:,:,lp);

       Omega0(l,lp)=GaussianD.PDF(mu(:,l),mu(:,lp),Sigma);
       Omega1(l,lp)=GaussianD.PDF(mu(:,l),mu(:,lp),Sigma+H);
       Omega2(l,lp)=GaussianD.PDF(mu(:,l),mu(:,lp),Sigma+2*H);
    end
end

MISE=(1/n)*(4*pi)^(-d/2)*1/sqrt(det(H))+w(:)'*((1-1/n)*Omega2-2*Omega1+Omega0)*w(:);

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:
close all;
clear all;
clc;

data_file_path = 'bin/mmv_phase_transition_snr_30db_s_16.mat';
options.export = true;
options.export_dir = 'bin';
options.export_name = 'mmv_snr_30_db_s_16';
options.chosen_ks = [2, 4, 8, 16, 32, 64];
options.subtitle = 'MMV, SNR=30dB s=16';
spx.pursuit.PhaseTransitionAnalysis.print_results(data_file_path, ...
    'CoSaMP', options);

author: indigits

repo: sparse-plex


text:
function c = intersect(a,b)
%INTERSECT    Intersection of intervals; empty components set to NaN
%
%   c = intersect(a,b)
%
%Result is an ***outer*** inclusion:
%   c includes the true intersection of a and b
%   components of c are set to NaN, if the corresponding components of
%     a and b are NaN, or if the intersection is empty
%   Thus a component NaN of the result c does NOT mean an empty intersection,
%     but rather that no information is available.
%
%===> To check empty intersection please use the INTLAB function
%===>   "EmptyIntersect"
%
%Input a and b must be both real or both complex
%

% written  10/16/98     S.M. Rump
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 12/06/05     S.M. Rump  correction of result (thanks to Andreas Rauh, Ulm, for hint)
% modified 06/11/06     S.M. Rump  correction of result (thanks to Andreas Rauh, Ulm, for hint)
% modified 09/10/07     S.M. Rump  redesign
% modified 04/22/09     S.M. Rump  comment changed, based on emptyintersect
%

  [empty,c] = emptyintersect(a,b);

repo: OpenMAS


text:
function trx = ReadMWTBlobFile(filename,trx,id)

% multiple files input, then call recursively
if iscell(filename),
  if nargin < 2,
    trx = [];
  end
  if nargin < 3,
    id = [];
  end
  for i1 = 1:numel(filename),
    trx = ReadMWTBlobFile(filename{i1},trx,id);
  end
  return;
end

nparams = 10;
nspinepts = 11;
contouroffset = '0';

if nargin < 2 || isempty(trx),
  trx = [];
  ids = [];
  havereadcontour = false;
  havereadspine = false;
else
  ids = [trx.id];
  havereadcontour = isfield(trx,'xcontour');
  havereadspine = isfield(trx,'xspine');
end

blobsext = '.blobs';
[~,~,ext] = fileparts(filename);
isblobs = strcmp(blobsext,ext);

fid = fopen(filename,'r');
if fid < 1,
  error('Could not open file %s for reading',filename);
end

% read in the id
if isblobs,
  while true,
    id = nan;
    isdone = false;
    while true,
      l1 = fgetl(fid);
      if ~ischar(l1),
        isdone = true;
        break;
      end
      if isempty(l1),
        continue;
      end
      if l1(1) == '%',
        id = str2double(l1(2:end));
      end
      break;
    end
    if isdone,
      break;
    end
    if isnan(id),
      warning('Could not parse id for line >%s<',l1);
      continue;
    end
    ReadMWTBlobFile1();
  end
else
  if nargin < 3 || isempty(id),
    id = 1;
  end
  ReadMWTBlobFile1();
end
%
% % set nframes
% for i1 = 1:numel(trx),
%   trx(i1).nframes = trx(i1).endframe - trx(i1).firstframe + 1;
% end

% sort by ids
[ids,order] = sort(ids);
trx = trx(order);

% check for missing data
for i1 = 1:numel(trx),
  if any(isnan(trx(i1).timestamps)),
    warning('%d frames of data not set for target %d',nnz(isnan(trx(i1).timestamps)),i1);
  end
end

fclose(fid);

  function ReadMWTBlobFile1()

    idi = find(ids == id,1);
    if ~isempty(idi),
      warning('Reading data for id %d more than once',id);
    end

    while true,
      filepos = ftell(fid);
      l = fgetl(fid);
      if ~ischar(l),
        break;
      end
      if isempty(l),
        continue;
      end
      % new id
      if isblobs && l(1) == '%',
        fseek(fid,filepos,'bof');
        break;
      end

      % parse the parameters in the data line
      [data_curr,count,errmsg,nextindex] = sscanf(l,'%f',nparams);
      if count < nparams || ~isempty(errmsg),
        warning('Could not parse data for id %d, line >%s<, error message: %s',id,l,errmsg);
        continue;
      end
      l = l(nextindex:end);

      % read in the contour, spine if it exists
      readcontour = false;
      readspine = false;
      while ~isempty(l),
        [w,count,~,nextindex] = sscanf(l,'%s',1);
        if count == 1,
          l = l(nextindex:end);
          switch w,
            case '%',
              [tmp,count,~,nextindex] = sscanf(l,'%f',nspinepts*2);
              if count == nspinepts*2,
                readspine = true;
                xspine = tmp(1:2:end-1);
                yspine = tmp(2:2:end);
              else
                warning('Error parsing spine points');
              end
              l = l(nextindex:end);
            case '%%',
              [tmp,count,~,nextindex] = sscanf(l,'%f',3);
              l = l(nextindex:end);
              if count == 3,
                xstart = tmp(1);
                ystart = tmp(2);
                ncontourpts = tmp(3);
                [a,count,~,nextindex] = sscanf(l,'%s',1);
                l = l(nextindex:end);
                if count == 1,
                  % convert to binary
                  b = dec2bin(a - contouroffset,8) == '1';
                  b = b(:,3:end)';
                  % first bit is 1 when y is changing
                  isy = b(1:2:end-1);
                  % second bit is 1 when positive change
                  ispos = b(2:2:end);
                  isy = isy(1:ncontourpts-1);
                  ispos = ispos(1:ncontourpts-1);
                  dx = double(~isy);
                  dx(~ispos) = -dx(~ispos);
                  dy = double(isy);
                  dy(~ispos) = -dy(~ispos);
                  xcontour = cumsum([xstart,dx]);
                  ycontour = cumsum([ystart,dy]);
                  readcontour = true;
                end
              end
            otherwise,
              warning('Unknown delimiter %s',w);
          end
        end
      end

      if ~readcontour && havereadcontour,
        xcontour = {};
        ycontour = {};
      end
      if ~readspine && havereadspine,
        xspine = nan(nspinepts,1);
        yspine = nan(nspinepts,1);
      end

      f = data_curr(1);
      timestamps = data_curr(2);
      x = data_curr(3);
      y = data_curr(4);
      area = data_curr(5);
      vx = data_curr(6);
      vy = data_curr(7);
      a = sqrt(vx^2+vy^2);
      theta = atan2(vy,vx);
      b = data_curr(8);
      length = data_curr(9);
      width = data_curr(10);

      if readspine,
        xspine = xspine + x;
        yspine = yspine + y;
      end

      % first time reading a contour
      if readcontour && ~havereadcontour,
        for tmpi = 1:numel(trx),
          nframescurr = trx(tmpi).endframe - trx(tmpi).firstframe + 1;
          trx(tmpi).xcontour = cell(1,nframescurr);
          trx(tmpi).ycontour = cell(1,nframescurr);
        end
      end
      % first time reading a spine
      if readspine && ~havereadspine,
        for tmpi = 1:numel(trx),
          nframescurr = trx(tmpi).endframe - trx(tmpi).firstframe + 1;
          trx(tmpi).xspine = nan(nspinepts,nframescurr);
          trx(tmpi).yspine = nan(nspinepts,nframescurr);
        end
      end

      if isempty(idi),
        idi = numel(ids) + 1;
        newtrx = struct('x',x,'y',y,'theta',theta,'a',a,'b',b,'id',id,...
          'firstframe',f,'nframes',nan,'endframe',f,'off',1-f,'timestamps',timestamps,...
          'area',area,'length',length,'width',width);
        if readcontour,
          newtrx.xcontour = {xcontour(:)};
          newtrx.ycontour = {ycontour(:)};
        elseif havereadcontour,
          newtrx.xcontour = {};
          newtrx.ycontour = {};
        end
        if readspine,
          newtrx.xspine = xspine(:);
          newtrx.yspine = yspine(:);
        elseif havereadspine,
          newtrx.xspine = nan(nspinepts,1);
          newtrx.yspine = nan(nspinepts,1);
        end

        if isempty(ids),
          trx = newtrx;
        else
          trx = structarrayset(trx,idi,newtrx);
        end
        ids(idi) = id;

      else

        % assume everything comes in order of frames, but some frames may
        % be skipped
        npad = f-trx(idi).endframe-1;
        if npad > 0,
          trx(idi).x(end+1:end+npad) = nan;
          trx(idi).y(end+1:end+npad) = nan;
          trx(idi).a(end+1:end+npad) = nan;
          trx(idi).b(end+1:end+npad) = nan;
          trx(idi).theta(end+1:end+npad) = nan;
          trx(idi).length(end+1:end+npad) = nan;
          trx(idi).width(end+1:end+npad) = nan;
          trx(idi).area(end+1:end+npad) = nan;
          trx(idi).timestamps(end+1:end+npad) = nan;

          if readcontour || havereadcontour,
            tmp = cell(1,npad);
            trx(idi).xcontour = [trx(idi).xcontour,tmp];
            trx(idi).ycontour = [trx(idi).ycontour,tmp];
          end

          if readspine || havereadspine,
            trx(idi).xspine(:,end+1:end+npad) = nan;
            trx(idi).yspine(:,end+1:end+npad) = nan;
          end
        end

        i = f - trx(idi).firstframe + 1;
        trx(idi).x(i) = x;
        trx(idi).y(i) = y;
        trx(idi).a(i) = a;
        trx(idi).b(i) = b;
        trx(idi).theta(i) = theta;
        trx(idi).length(i) = length;
        trx(idi).width(i) = width;
        trx(idi).area(i) = area;
        trx(idi).timestamps(i) = timestamps;

        if readcontour || havereadcontour,
          trx(idi).xcontour{i} = xcontour;
          trx(idi).ycontour{i} = ycontour;
        end

        if readspine || havereadspine,
          trx(idi).xspine(:,i) = xspine;
          trx(idi).yspine(:,i) = yspine;
        end
        trx(idi).endframe = max(f,trx(idi).endframe);

      end

      havereadcontour = havereadcontour || readcontour;
      havereadspine = havereadspine || readspine;

    end

    if isempty(idi),
      %fprintf('Read 0 frames for id %d\n',id);
    else
      trx(idi).nframes = trx(idi).endframe-trx(idi).firstframe+1;
      %fprintf('Read %d frames for id %d\n',trx(idi).nframes,id);
    end

    if havereadcontour,
      for tmpi = 1:numel(trx),
        nframescurr = trx(tmpi).nframes;
        nadd = nframescurr - numel(trx(tmpi).xcontour);
        if nadd > 0,
          trx(tmpi).xcontour(end+1:end+nadd) = cell(1,nadd);
          trx(tmpi).ycontour(end+1:end+nadd) = cell(1,nadd);
        end
      end
    end
    if havereadspine,
      for tmpi = 1:numel(trx),
        nframescurr = trx(tmpi).nframes;
        nadd = nframescurr - size(trx(tmpi).xspine,2);
        if nadd > 0,
          trx(tmpi).xspine(:,end+1:end+nadd) = nan;
          trx(tmpi).yspine(:,end+1:end+nadd) = nan;
        end
      end
    end

  end

end

author: kristinbranson

repo: JAABA


text:
%% Copyright (C) 2014-2016, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym asind (@var{x})
%% Symbolic inverse sin function with output in degrees.
%%
%% Example:
%% @example
%% @group
%% asind (sqrt (sym (2))/2)
%%   @result{} (sym) 45
%%
%% syms x
%% y = asind (x)
%%   @result{} y = (sym)
%%       180⋅asin(x)
%%       ───────────
%%            π
%% @end group
%% @end example
%%
%% @seealso{@@sym/sind, @@sym/asin}
%% @end defmethod

function y = asind(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('lambda a: deg(asin(a))', x);
end

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

%!test
%! f1 = asind (sym(1)/2);
%! f2 = asind (1/2);
%! assert (double (f1), f2, -eps)

%!test
%! D = [1 2; 3 4]/4;
%! A = sym([1 2; 3 4])/4;
%! f1 = asind (A);
%! f2 = asind (D);
%! assert (double (f1), f2, -eps)

repo: octsympy