text:
function binfeatures = derivebinaryfeat(randf, Tr_Data, params, stage)
%DERIVEBINARYFEAT Summary of this function goes here
%   Function: Derive binary features for each sample given learned random forest
%   Detailed explanation goes here
%   Input:
%       lmarkID: the ID of landmark
%       randf:   learned random forest in current stage
%       Tr_Data: training data
%       params:  parameters for curent stage

% calculate the overall dimension of binary feature, concatenate the
% features for all landmarks, and the feature of one landmark is sum
% leafnodes of all random trees;

dims_binfeat = 0;

ind_bincode = zeros(size(randf));

for l = 1:size(randf, 1)
    for t = 1:size(randf, 2)
        ind_bincode(l, t) = randf{l, t}.num_leafnodes;
        dims_binfeat = dims_binfeat + randf{l, t}.num_leafnodes;
    end
end

% initilaize the memory for binfeatures
dbsize = length(Tr_Data);

% faster implementation
%{
tic;
binfeatures   = zeros(dbsize*(params.augnumber), dims_binfeat);
img_sizes     = zeros(2, dbsize*(params.augnumber));
shapes        = zeros([size(params.meanshape), dbsize*(params.augnumber)]);
tfs2meanshape = cell(1, dbsize*params.augnumber);

img_areas = zeros(1, dbsize);
parfor i = 1:dbsize
    img_areas(i) = Tr_Data{i}.width*Tr_Data{i}.height;
end

imgs_gray = zeros(dbsize, max(img_areas));

for i = 1:dbsize
    area = img_areas(i);
    imgs_gray(i, 1:area) = Tr_Data{i}.img_gray(:)';
    shapes(:, :, (i-1)*params.augnumber+1:i*params.augnumber) = Tr_Data{i}.intermediate_shapes{stage};

    for k = 1:params.augnumber
        img_sizes(:, (i-1)*params.augnumber+k) = [Tr_Data{i}.width; Tr_Data{i}.height];
        tfs2meanshape{(i-1)*params.augnumber+k} = Tr_Data{i}.tf2meanshape{k};
    end
end

binfeature_lmarks = cell(1, size(params.meanshape, 1));
num_leafnodes = zeros(1, size(params.meanshape, 1));
for l = 1:size(params.meanshape, 1)
    [binfeature_lmarks{l}, num_leafnodes(l)] = lbf_faster(randf{l}, imgs_gray, img_sizes, shapes(l, :, :), tfs2meanshape, params, stage);
end

cumnum_leafnodes = [0 cumsum(num_leafnodes)];
binfeatures_faster = binfeatures;
for l = 1:size(params.meanshape, 1)
    binfeatures_faster(:, cumnum_leafnodes(l)+1:cumnum_leafnodes(l+1)) = binfeature_lmarks{l};
end

toc;
%}
feats = cell(size(params.meanshape, 1), 1);
isleaf = cell(size(params.meanshape, 1), 1);
threshs = cell(size(params.meanshape, 1), 1);
cnodes = cell(size(params.meanshape, 1), 1);

% prepare for the derivation of local binary codes
for l = 1:size(params.meanshape, 1)
    % concatenate all nodes of all random trees
    rf = randf(l, :);
    num_rfnodes = 0;
    for t = 1:params.max_numtrees
        num_rfnodes = num_rfnodes + rf{t}.num_nodes;
    end

    % fast implementation
    feats{l} = zeros(num_rfnodes, 4);
    isleaf{l} = zeros(num_rfnodes, 1);
    threshs{l} = zeros(num_rfnodes, 1);
    cnodes{l} = zeros(num_rfnodes, 2);

    id_rfnode = 1;

    for t = 1:params.max_numtrees
        feats{l}(id_rfnode:(id_rfnode + rf{t}.num_nodes - 1), :)   = rf{t}.feat(1:rf{t}.num_nodes, :);
        isleaf{l}(id_rfnode:(id_rfnode + rf{t}.num_nodes - 1), :)  = rf{t}.isleafnode(1:rf{t}.num_nodes, :);
        threshs{l}(id_rfnode:(id_rfnode + rf{t}.num_nodes - 1), :) = rf{t}.thresh(1:rf{t}.num_nodes, :);
        cnodes{l}(id_rfnode:(id_rfnode + rf{t}.num_nodes - 1), :)  = rf{t}.cnodes(1:rf{t}.num_nodes, :);

        id_rfnode = id_rfnode + rf{t}.num_nodes;
    end
end

% extract feature for each samples
parfor i = 1:dbsize*(params.augnumber)
    k = floor((i-1)/(params.augnumber)) + 1;
    s = mod(i-1, (params.augnumber)) + 1;
    img_gray = Tr_Data{k}.img_gray;
    bbox     = Tr_Data{k}.intermediate_bboxes{stage}(s, :);
    shapes   = Tr_Data{k}.intermediate_shapes{stage};
    shape    = shapes(:, :, s);

    % tf2meanshape = Tr_Data{k}.tf2meanshape{s};
    meanshape2tf = Tr_Data{k}.meanshape2tf{s};
    %{
    feats_cell   = cell(size(params.meanshape, 1), params.max_numtrees);

    for l = 1:size(params.meanshape, 1)
        % extract feature for each landmark from the random forest
        for t = 1:params.max_numtrees
            % extract feature for each ladnmark from a random tree
            bincode = traversetree(randf{l}{t}, img_gray, bbox, shape(l, :), tf2meanshape, params, stage);
            feats_cell{l, t} = bincode;
        end
    end

    binfeature = zeros(1, dims_binfeat);
    for l = 1:size(params.meanshape, 1)
        % extract feature for each landmark from the random forest
        offset = sum(ind_bincodes(1:l-1, end));
        for t = 1:params.max_numtrees
            ind_s = ind_bincodes(l, t) + 1 + offset;
            ind_e = ind_bincodes(l, t+1) + offset;
            % extract feature for each ladnmark from a random tree
            binfeature(ind_s:ind_e) = feats_cell{l, t};
        end
    end
    %}
    % fast implementation

    binfeature_lmarks = cell(1, size(params.meanshape, 1));
    num_leafnodes = zeros(1, size(params.meanshape, 1));
    for l = 1:size(params.meanshape, 1)
        % concatenate all nodes of all random trees
        [binfeature_lmarks{l}, num_leafnodes(l)]= lbf_fast(randf(l, :), feats{l}, isleaf{l}, threshs{l}, cnodes{l}, img_gray, bbox, shape(l, :), meanshape2tf, params, stage);
    end

    cumnum_leafnodes = [0 cumsum(num_leafnodes)];

    binfeature_alllmarks = zeros(1, cumnum_leafnodes(end));
    for l = 1:size(params.meanshape, 1)
        binfeature_alllmarks(cumnum_leafnodes(l)+1:cumnum_leafnodes(l+1)) = binfeature_lmarks{l};
    end

    %{
    % faster implementation (thinking...)
    binfeature_lmarks = cell(1, size(params.meanshape, 1));
    num_leafnodes = zeros(1, size(params.meanshape, 1));
    for l = 1:size(params.meanshape, 1)
        % concatenate all nodes of all random trees
        [binfeature_lmarks{l}, num_leafnodes(l)]= lbf_faster(randf{l}, img_gray, bbox, shape(l, :), tf2meanshape, params, stage);
    end

    cumnum_leafnodes = [0 cumsum(num_leafnodes)];

    binfeature_alllmarks = zeros(1, cumnum_leafnodes(end));
    for l = 1:size(params.meanshape, 1)
        binfeature_alllmarks(cumnum_leafnodes(l)+1:cumnum_leafnodes(l+1)) = binfeature_lmarks{l};
    end
    %}

    binfeatures(i, :) = binfeature_alllmarks;
end

end

function bincode = traversetree(tree, img_gray, bbox, shape, tf2meanshape, params, stage)

currnode = tree.rtnodes{1};

bincode = zeros(1, tree.num_leafnodes);

width = size(img_gray, 2);
height = size(img_gray, 1);

while(1)

    anglepairs = currnode.feat(1:2);
    radiuspairs = currnode.feat(3:4);

    % calculate the relative location under the coordinate of meanshape
    pixel_a_x_imgcoord = cos(anglepairs(:, 1)).*radiuspairs(:, 1)*params.max_raio_radius(stage)*bbox(3);
    pixel_a_y_imgcoord = sin(anglepairs(:, 1)).*radiuspairs(:, 1)*params.max_raio_radius(stage)*bbox(4);

    pixel_b_x_imgcoord = cos(anglepairs(:, 2)).*radiuspairs(:, 2)*params.max_raio_radius(stage)*bbox(3);
    pixel_b_y_imgcoord = sin(anglepairs(:, 2)).*radiuspairs(:, 2)*params.max_raio_radius(stage)*bbox(4);

    % transform the pixels from image coordinate (meanshape) to coordinate of current shape
    [pixel_a_x_lmcoord, pixel_a_y_lmcoord] = tforminv(tf2meanshape, pixel_a_x_imgcoord, pixel_a_y_imgcoord);

    [pixel_b_x_lmcoord, pixel_b_y_lmcoord] = tforminv(tf2meanshape, pixel_b_x_imgcoord, pixel_b_y_imgcoord);

    pixel_a_x = ceil(pixel_a_x_lmcoord + shape(1));
    pixel_a_y = ceil(pixel_a_y_lmcoord + shape(2));

    pixel_b_x = ceil(pixel_b_x_lmcoord + shape(1));
    pixel_b_y = ceil(pixel_b_y_lmcoord + shape(2));

    pixel_a_x = max(1, min(pixel_a_x, width));
    pixel_a_y = max(1, min(pixel_a_y, height));

    pixel_b_x = max(1, min(pixel_b_x, width));
    pixel_b_y = max(1, min(pixel_b_y, height));

    f = int16(img_gray(pixel_a_y + (pixel_a_x-1)*height)) - int16(img_gray(pixel_b_y + (pixel_b_x-1)*height));

    if f < (currnode.thresh)
        id_chilenode = currnode.cnodes(1);
        currnode = tree.rtnodes{id_chilenode};
    else
        id_chilenode = currnode.cnodes(2);
        currnode = tree.rtnodes{id_chilenode};
    end

    if isempty(currnode.cnodes)
        % if the current node is a leaf node, then stop traversin
        bincode(tree.id_leafnodes == id_chilenode) = 1;
        break;
    end
end

end

function [binfeature, num_leafnodes] = lbf_fast(rf, feats, isleaf, threshs, cnodes, img_gray, bbox, shape, meanshape2tf, params, stage)

%{
num_rfnodes = 0;
for t = 1:params.max_numtrees
    num_rfnodes = num_rfnodes + rf{t}.num_nodes;
end

% fast implementation
feats = zeros(num_rfnodes, 4);
isleaf = zeros(num_rfnodes, 1);
threshs = zeros(num_rfnodes, 1);
cnodes = zeros(num_rfnodes, 2);

id_rfnode = 1;

for t = 1:params.max_numtrees
    feats(id_rfnode:(id_rfnode + rf{t}.num_nodes - 1), :)   = rf{t}.feat(1:rf{t}.num_nodes, :);
    isleaf(id_rfnode:(id_rfnode + rf{t}.num_nodes - 1), :)  = rf{t}.isleafnode(1:rf{t}.num_nodes, :);
    threshs(id_rfnode:(id_rfnode + rf{t}.num_nodes - 1), :) = rf{t}.thresh(1:rf{t}.num_nodes, :);
    cnodes(id_rfnode:(id_rfnode + rf{t}.num_nodes - 1), :)  = rf{t}.cnodes(1:rf{t}.num_nodes, :);

    id_rfnode = id_rfnode + rf{t}.num_nodes;
end
%}

width = size(img_gray, 2);
height = size(img_gray, 1);

anglepairs = feats(:, 1:2);
radiuspairs = feats(:, 3:4);

% calculate the relative location under the coordinate of meanshape
pixel_a_x_imgcoord = cos(anglepairs(:, 1)).*radiuspairs(:, 1)*params.max_raio_radius(stage)*bbox(3);
pixel_a_y_imgcoord = sin(anglepairs(:, 1)).*radiuspairs(:, 1)*params.max_raio_radius(stage)*bbox(4);

pixel_b_x_imgcoord = cos(anglepairs(:, 2)).*radiuspairs(:, 2)*params.max_raio_radius(stage)*bbox(3);
pixel_b_y_imgcoord = sin(anglepairs(:, 2)).*radiuspairs(:, 2)*params.max_raio_radius(stage)*bbox(4);

% no transformaiton
%{
pixel_a_x_lmcoord = pixel_a_x_imgcoord;
pixel_a_y_lmcoord = pixel_a_y_imgcoord;

pixel_b_x_lmcoord = pixel_b_x_imgcoord;
pixel_b_y_lmcoord = pixel_b_y_imgcoord;
%}

% transform the pixels from image coordinate (meanshape) to coordinate of current shape

[pixel_a_x_lmcoord, pixel_a_y_lmcoord] = transformPointsForward(meanshape2tf, pixel_a_x_imgcoord', pixel_a_y_imgcoord');
pixel_a_x_lmcoord = pixel_a_x_lmcoord';
pixel_a_y_lmcoord = pixel_a_y_lmcoord';

[pixel_b_x_lmcoord, pixel_b_y_lmcoord] = transformPointsForward(meanshape2tf, pixel_b_x_imgcoord', pixel_b_y_imgcoord');
pixel_b_x_lmcoord = pixel_b_x_lmcoord';
pixel_b_y_lmcoord = pixel_b_y_lmcoord';

pixel_a_x = ceil(pixel_a_x_lmcoord + shape(1));
pixel_a_y = ceil(pixel_a_y_lmcoord + shape(2));

pixel_b_x = ceil(pixel_b_x_lmcoord + shape(1));
pixel_b_y = ceil(pixel_b_y_lmcoord + shape(2));

pixel_a_x = max(1, min(pixel_a_x, width));
pixel_a_y = max(1, min(pixel_a_y, height));

pixel_b_x = max(1, min(pixel_b_x, width));
pixel_b_y = max(1, min(pixel_b_y, height));

pdfeats = double(img_gray(pixel_a_y + (pixel_a_x-1)*height)) - double(img_gray(pixel_b_y + (pixel_b_x-1)*height));
    % ./ double(img_gray(pixel_a_y + (pixel_a_x-1)*height)) + double(img_gray(pixel_b_y + (pixel_b_x-1)*height));

cind = (pdfeats >= threshs) + 1;

% obtain the indice of child nodes for all nodes, if the current node is
% leaf node, then the indice of its child node is 0
ind_cnodes = cnodes; % diag(cnodes(:, cind));

binfeature = zeros(1, sum(isleaf));

cumnum_nodes = 0;
cumnum_leafnodes = 0;
for t = 1:params.max_numtrees
    num_nodes = (rf{t}.num_nodes);
    id_cnode = 1;
    while(1)
        if isleaf(id_cnode + cumnum_nodes)
            binfeature(cumnum_leafnodes + find(rf{t}.id_leafnodes == id_cnode)) = 1;
            cumnum_nodes = cumnum_nodes + num_nodes;
            cumnum_leafnodes = cumnum_leafnodes + rf{t}.num_leafnodes;
            break;
        end
        id_cnode = ind_cnodes(cumnum_nodes + id_cnode, cind(cumnum_nodes + id_cnode));
    end
end

num_leafnodes = sum(isleaf);

end

author: jwyang

repo: face-alignment


text:
function comparisonTable = compareInchis(model, inchis, met)
% Compare inchi strings. Each inchi string is given a score based on its similarity
% to the chemical formula and charge of the metabolite in the model. Factors such
% as stereochemistry, if it is a standard inchi, and its similarity to the other
% inchis are also considered.
%
% USAGE:
%
%    comparisonTable = compareInchis(model, inchis, met)
%
% INPUT:
%    model:     COBRA model
%                   *.mets - 1xm list of metabolite IDs
%                   *.metFormulas - 1xm list of metabolite Formulas
%                   *.metCharges - 1xm list of metabolite charges
%    inchis:    List of InChIs to compare
%    met:       Metabolite ID in the model
%
% OUTPUTS:
%    comparisonTable: Table containing the information used for each InChI
%

if isrow(inchis)
    inchis = inchis';
end
mets = regexprep(model.mets, '(\[\w\])', '');

% Prepare the table with inchi data
nRows = length(inchis);
varTypes = {'double', 'logical', 'string', 'string', 'logical', 'double', 'logical',...
    'logical', 'logical', 'double', 'double', 'logical', 'double', 'logical', 'string'};
varNames = {'scores', 'rGroup', 'InChI', 'metFormula', 'formulaOkBool', 'netCharge',...
    'chargeOkBool', 'stereochemicalSubLayers', 'standard', 'sourceSimilarity', ...
    'mainLayerSimilarity', 'isotopicLayer', 'layers', 'inchiWithMoreLayers', 'mainLayer'};
comparisonTable = table('Size', [nRows length(varTypes)], 'VariableTypes', varTypes,...
    'VariableNames', varNames);

emptyBool = cellfun(@isempty, inchis);
comparisonTable.InChI = inchis;

% Get InChI data
for i = 1:length(inchis)
    if ~emptyBool(i)
        inchiLayersDetail = getInchiData(inchis{i});
        comparisonTable.layers(i) = inchiLayersDetail.layers;
        comparisonTable.standard(i) = inchiLayersDetail.standard;
        comparisonTable.metFormula(i) = inchiLayersDetail.metFormula;
        comparisonTable.mainLayer(i) = inchiLayersDetail.mainLayer;
        comparisonTable.netCharge(i) = inchiLayersDetail.netCharge;
        comparisonTable.stereochemicalSubLayers(i) = inchiLayersDetail.stereochemicalSubLayers;
        comparisonTable.isotopicLayer(i) = inchiLayersDetail.isotopicLayer;
    end
end

% Charge comparison
% Compares the net charge in the InChI with the charge in the corresponding
% metabolite in the model.
if isfield(model, 'metCharges')
    modelMetCharge = unique(model.metCharges(ismember(mets, met)));
    comparisonTable.chargeOkBool = ~emptyBool & comparisonTable.netCharge == modelMetCharge;
end
if isfield(model, 'metCharge')
    modelMetCharge = unique(model.metCharge(ismember(mets, met)));
    comparisonTable.chargeOkBool = ~emptyBool & comparisonTable.netCharge == modelMetCharge;
end

% Formula comparison
% Compares the formula with the formula in the corresponding metabolite in the
% model. Hydrogens should be ignored because they are considered in the charge
% or isotopic layers.
modelMetFormula = unique(model.metFormulas(ismember(mets, met)));
% Remove the R groups from the formula
rGroup = ["X", "Y", "*", "FULLR"];
if contains(modelMetFormula, rGroup)
    if isempty(regexprep(modelMetFormula, rGroup, ''))
        modelMetFormula = 'H';
    else
        modelMetFormula = editChemicalFormula(modelMetFormula);
    end
end
% Check if the formulas match
comparisonTable.formulaOkBool = ismember(regexprep(comparisonTable.metFormula, ...
    'H\d*', ''), regexprep(modelMetFormula, 'H\d*', ''));

% Source similarity
% For each InChI is calculated the number of equal InChIs divided between
% the number of InChIs
inchis(emptyBool) = {'noData'};
for i = 1:length(inchis)
    if ~emptyBool(i)
        comparisonTable.sourceSimilarity(i) = sum(ismember(inchis, inchis(i))) / sum(~emptyBool);
    end
end

% Main layer similarity
% For each InChI is calculated the the number of equal main layer similarity
% divided between the number of InChIs
for i = 1:length(inchis)
    if ~emptyBool(i)
        comparisonTable.mainLayerSimilarity(i) = sum(ismember(comparisonTable.mainLayer, ...
            comparisonTable.mainLayer(i))) / sum(~emptyBool);
    end
end

% InChI with more layers (stereochemistry and charge)
comparisonTable.inchiWithMoreLayers = comparisonTable.layers == max(comparisonTable.layers);

% Final score
% Priority will be given to molecules that have the same chemical formula as
% the corresponding metabolite in the model, excluding hydrogen and non-chemical
% atoms.
comparisonTable.scores = comparisonTable.chargeOkBool + (comparisonTable.formulaOkBool * 10)...
    + comparisonTable.stereochemicalSubLayers + comparisonTable.standard + ...
    comparisonTable.sourceSimilarity + comparisonTable.mainLayerSimilarity + ...
    comparisonTable.inchiWithMoreLayers;

end

author: opencobra

repo: cobratoolbox


text:
function G =  getGroupOverlapColor27(row, col)

    N = row*col;

    g = sparse(zeros(N,1));
    g = diag(g);

    %% build overlapping group

    % top let corner group

    g(1, 1 ) = 1;
    g(1, 2 ) = 1;
    g(1, 1+row ) = 1;
    g(1, 2+row ) = 1;

     % bottom let corner group

    g(row, row-1 ) = 1;
    g(row, row ) = 1;
    g(row, row-1+row ) = 1;
    g(row, row+row ) = 1;

    % top right corner group

    g((col-1)*row+1, (col-1)*row+1 ) = 1;
    g((col-1)*row+1, (col-2)*row+1 ) = 1;
    g((col-1)*row+1, (col-1)*row+2 ) = 1;
    g((col-1)*row+1, (col-2)*row+2 ) = 1;

    % bottom right corner group

    g((col-1)*row+row , (col-1)*row+row-1 ) = 1;
    g((col-1)*row+row , (col-2)*row+row-1 ) = 1;
    g((col-1)*row+row , (col-1)*row+row ) = 1;
    g((col-1)*row+row , (col-2)*row+row ) = 1;

    % boundary group

    for i=2:col-1
        % top row rgoup

        j = 1;

        g( (i-1)*row+j  , (i-2)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j+1 ) = 1;

        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;

        g( (i-1)*row+j  , i*row+j   ) = 1;
        g( (i-1)*row+j  , i*row+j+1 ) = 1;

        % bottom row group

        j = row;

        g( (i-1)*row+j  , (i-2)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j   ) = 1;

        g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;

        g( (i-1)*row+j  , i*row+j-1 ) = 1;
        g( (i-1)*row+j  , i*row+j   ) = 1;

    end

    for j=2:row-1
          % left column rgoup

        i=1;

        g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;
        g( (i-1)*row+j  , i*row+j-1 ) = 1;
        g( (i-1)*row+j  , i*row+j   ) = 1;
        g( (i-1)*row+j  , i*row+j+1 ) = 1;

        % right  column group

        i=col;

        g( (i-1)*row+j  , (i-2)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-2)*row+j+1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
        g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;

    end

    for i=2:col-1
        for j=2:row-1

            g( (i-1)*row+j  , (i-2)*row+j-1 ) = 1;
            g( (i-1)*row+j  , (i-2)*row+j   ) = 1;
            g( (i-1)*row+j  , (i-2)*row+j+1 ) = 1;
            g( (i-1)*row+j  , (i-1)*row+j-1 ) = 1;
            g( (i-1)*row+j  , (i-1)*row+j   ) = 1;
            g( (i-1)*row+j  , (i-1)*row+j+1 ) = 1;
            g( (i-1)*row+j  , i*row+j-1 ) = 1;
            g( (i-1)*row+j  , i*row+j   ) = 1;
            g( (i-1)*row+j  , i*row+j+1 ) = 1;

        end
    end

     G =[g; g; g];

end

author: andrewssobral

repo: lrslibrary


text:
function Res = dfm(X,Spec,threshold)
%DFM()    Runs the dynamic factor model
%
%  Syntax:
%    Res = DFM(X,Par)
%
%  Description:
%   DFM() inputs the organized and transformed data X and parameter structure Par.
%   Then, the function outputs dynamic factor model structure Res and data
%   summary statistics (mean and standard deviation).
%
%  Input arguments:
%    X: Kalman-smoothed data where missing values are replaced by their expectation
%    Par: A structure containing the following parameters:
%      Par.blocks: Block loadings.
%      Par.nQ: Number of quarterly series
%      Par.p: Number of lags in transition matrix
%      Par.r: Number of common factors for each block
%
% Output Arguments:
%
%   Res - structure of model results with the following fields
%       . X_sm | Kalman-smoothed data where missing values are replaced by their expectation
%       . Z | Smoothed states. Rows give time, and columns are organized according to Res.C.
%       . C | Observation matrix. The rows correspond
%          to each series, and the columns are organized as shown below:
%         - 1-20: These columns give the factor loadings. For example, 1-5
%              give loadings for the first block and are organized in
%              reverse-chronological order (f^G_t, f^G_t-1, f^G_t-2, f^G_t-3,
%              f^G_t-4). Columns 6-10, 11-15, and 16-20 give loadings for
%              the second, third, and fourth blocks respectively.
%       .R: Covariance for observation matrix residuals
%       .A: Transition matrix. This is a square matrix that follows the
%      same organization scheme as Res.C's columns. Identity matrices are
%      used to account for matching terms on the left and righthand side.
%      For example, we place an I4 matrix to account for matching
%      (f_t-1; f_t-2; f_t-3; f_t-4) terms.
%       .Q: Covariance for transition equation residuals.
%       .Mx: Series mean
%       .Wx: Series standard deviation
%       .Z_0: Initial value of state
%       .V_0: Initial value of covariance matrix
%       .r: Number of common factors for each block
%       .p: Number of lags in transition equation
%
% References:
%
%   Marta Banbura, Domenico Giannone and Lucrezia Reichlin
%   Nowcasting (2010)
%   Michael P. Clements and David F. Hendry, editors,
%   Oxford Handbook on Economic Forecasting.

%% Store model parameters ------------------------------------------------

% DFM input specifications: See documentation for details
Par.blocks = Spec.Blocks;                  % Block loading structure
Par.nQ = sum(strcmp('q',Spec.Frequency));  % Number of quarterly series
Par.p = 1;                                 % Number of lags in autoregressive of factor (same for all factors)
Par.r = ones(1,size(Spec.Blocks,2));       % Number of common factors for each block
%Par.r(1) =2;
% Display blocks
try
    fprintf('\n\n\n');
    disp('Table 3: Block Loading Structure')
    disp(array2table(Spec.Blocks,...
         'RowNames', strrep(Spec.SeriesName,' ','_'),...
         'VariableNames',Spec.BlockNames));
    fprintf('\n\n\n');
catch
end

fprintf('Estimating the dynamic factor model (DFM) ... \n\n');

[T,N] = size(X);
r = Par.r;
p = Par.p;
nQ = Par.nQ;
blocks = Par.blocks;

i_idio = logical([ones(N-nQ,1);zeros(nQ,1)]);

%R*Lambda = q; Contraints on the loadings of the quartrly variables

R_mat = [  2 -1  0  0  0;...
           3  0 -1  0  0;...
           2  0  0 -1  0;...
           1  0  0  0 -1];

q = zeros(4,1);

if(nargin < 3)
    threshold = 1e-5;  % EM loop threshold (default value)
end

max_iter = 5000;  % EM loop maximum number of iterations

%% Prepare data -----------------------------------------------------------

Mx = mean(X,'omitnan');
Wx = std(X,'omitnan');
xNaN = (X-repmat(Mx,T,1))./repmat(Wx,T,1);  % Standardize series

%% Initial Conditions -----------------------------------------------------

optNaN.method = 2; % Remove leading and closing zeros
optNaN.k = 3;      % Setting for filter(): See remNaN_spline

[A, C, Q, R, Z_0, V_0] = InitCond(xNaN,r,p,blocks,optNaN,R_mat,q,nQ,i_idio);

% Initialize EM loop values
previous_loglik = -inf;
num_iter = 0;
LL = -inf;
converged = 0;

% y for the estimation is WITH missing data
y = xNaN';

%% EM LOOP ----------------------------------------------------------------

%The model can be written as
%y = C*Z + e;
%Z = A*Z(-1) + v
%where y is NxT, Z is (pr)xT, etc

% Remove the leading and ending nans
optNaN.method = 3;
y_est = remNaNs_spline(xNaN,optNaN)';

while (num_iter < max_iter) & ~converged % Loop until converges or max iter.

    [C_new, R_new, A_new, Q_new, Z_0, V_0, loglik] = ...  % Applying EM algorithm
        EMstep(y_est, A, C, Q, R, Z_0, V_0, r,p,R_mat,q,nQ,i_idio,blocks);

    C = C_new;
    R = R_new;
    A = A_new;
    Q = Q_new;

    if num_iter > 2  % Checking convergence
        [converged, decrease(num_iter + 1)] = ...
            em_converged(loglik, previous_loglik, threshold, 1);
    end

    if (mod(num_iter,10) == 0) && (num_iter > 0)  % Print updates to command window
        disp(['Now running the ',num2str(num_iter),...
              'th iteration of max ', num2str(max_iter)]);
        disp(['  Loglik','   (% Change)'])
        disp([num2str(loglik),'   (', sprintf('%6.2f',100*((loglik-previous_loglik)/previous_loglik)) '%)'])
    end

    LL = [LL loglik];
    previous_loglik = loglik;
    num_iter =  num_iter + 1;

end

if(num_iter < max_iter)
    disp(['Successful: Convergence at ', num2str(num_iter), ' iterations'])
else
   disp('Stopped because maximum iterations reached')
end

% Final run of the Kalman filter
Zsmooth = runKF(y, A, C, Q, R, Z_0, V_0)';

x_sm = Zsmooth(2:end,:) * C';  % Get smoothed X

%%  Loading the structure with the results --------------------------------
Res.x_sm = x_sm;
Res.X_sm = repmat(Wx,T,1) .* x_sm + repmat(Mx,T,1);  % Unstandardized, smoothed
Res.Z = Zsmooth(2:end,:);
Res.C = C;
Res.R = R;
Res.A = A;
Res.Q = Q;
Res.Mx = Mx;
Res.Wx = Wx;
Res.Z_0 = Z_0;
Res.V_0 = V_0;
Res.r = r;
Res.p = p;

%% Display output
% Table with names and factor loadings

nQ       = Par.nQ;                      % Number of quarterly series
nM       = size(Spec.SeriesID,1) - nQ;  % Number monthly series
nLags    = max(Par.p, 5);               % 5 comes from monthly-quarterly aggregation
nFactors = sum(Par.r);

fprintf('\n\n\n');

try
disp('Table 4: Factor Loadings for Monthly Series');
disp(array2table(Res.C(1:nM, 1:5:nFactors*5),...  % Only select lag(0) terms
     'RowNames', strrep(Spec.SeriesName(1:nM),' ','_'), ...
     'VariableNames', Spec.BlockNames));
fprintf('\n\n\n');
disp('Table 5: Quarterly Loadings Sample (Global Factor)')
disp(array2table(Res.C(end-nQ+1:end, 1:5), ...  % Select only quarterly series
     'RowNames', strrep(Spec.SeriesName(end-nQ+1:end),' ','_'), ...
     'VariableNames', {'f1_lag0', 'f1_lag1', 'f1_lag2', 'f1_lag3', 'f1_lag4'}));
fprintf('\n\n\n');
catch
end

% Table with AR model on factors (factors with AR parameter and variance of residuals)

A_terms = diag(Res.A);  % Transition equation terms
Q_terms = diag(Res.Q);  % Covariance matrix terms

try
disp('Table 6: Autoregressive Coefficients on Factors')
disp(table(A_terms(1:5:nFactors*5), ...  % Only select lag(0) terms
           Q_terms(1:5:nFactors*5), ...
           'VariableNames', {'AR_Coefficient', 'Variance_Residual'}, ...
           'RowNames',      strrep(Spec.BlockNames,' ','_')));
fprintf('\n\n\n');
catch
end

% Table with AR model idiosyncratic errors (factors with AR parameter and variance of residuals)
try
disp('Table 7: Autoregressive Coefficients on Idiosyncratic Component')
disp(table(A_terms([nFactors*5+1:nFactors*5+nM nFactors*5+nM+1:5:end]),...  % 21:50 give monthly, 51:5:61 give quarterly
           Q_terms([nFactors*5+1:nFactors*5+nM nFactors*5+nM+1:5:end]), ...
           'VariableNames', {'AR_Coefficient', 'Variance_Residual'}, ...
           'RowNames', strrep(Spec.SeriesName,' ','_')));
catch
end

end

%% PROCEDURES -------------------------------------------------------------
% Note: Kalman filter (runKF()) is in the 'functions' folder

function  [C_new, R_new, A_new, Q_new, Z_0, V_0, loglik] = EMstep(y, A, C, Q, R, Z_0, V_0, r,p,R_mat,q,nQ,i_idio,blocks)
%EMstep    Applies EM algorithm for parameter reestimation
%
%  Syntax:
%    [C_new, R_new, A_new, Q_new, Z_0, V_0, loglik]
%    = EMstep(y, A, C, Q, R, Z_0, V_0, r, p, R_mat, q, nQ, i_idio, blocks)
%
%  Description:
%    EMstep reestimates parameters based on the Estimation Maximization (EM)
%    algorithm. This is a two-step procedure:
%    (1) E-step: the expectation of the log-likelihood is calculated using
%        previous parameter estimates.
%    (2) M-step: Parameters are re-estimated through the maximisation of
%        the log-likelihood (maximize result from (1)).
%
%    See "Maximum likelihood estimation of factor models on data sets with
%    arbitrary pattern of missing data" for details about parameter
%    derivation (Banbura & Modugno, 2010). This procedure is in much the
%    same spirit.
%
%  Input:
%    y:      Series data
%    A:      Transition matrix
%    C:      Observation matrix
%    Q:      Covariance for transition equation residuals
%    R:      Covariance for observation matrix residuals
%    Z_0:    Initial values of factors
%    V_0:    Initial value of factor covariance matrix
%    r:      Number of common factors for each block (e.g. vector [1 1 1 1])
%    p:      Number of lags in transition equation
%    R_mat:  Estimation structure for quarterly variables (i.e. "tent")
%    q:      Constraints on loadings
%    nQ:     Number of quarterly series
%    i_idio: Indices for monthly variables
%    blocks: Block structure for each series (i.e. for a series, the structure
%            [1 0 0 1] indicates loadings on the first and fourth factors)
%
%  Output:
%    C_new: Updated observation matrix
%    R_new: Updated covariance matrix for residuals of observation matrix
%    A_new: Updated transition matrix
%    Q_new: Updated covariance matrix for residuals for transition matrix
%    Z_0:   Initial value of state
%    V_0:   Initial value of covariance matrix
%    loglik: Log likelihood
%
% References:
%   "Maximum likelihood estimation of factor models on data sets with
%   arbitrary pattern of missing data" by Banbura & Modugno (2010).
%   Abbreviated as BM2010
%
%

%% Initialize preliminary values

% Store series/model values
[n, T] = size(y);
nM = n - nQ;  % Number of monthly series
pC = size(R_mat,2);
ppC = max(p,pC);
num_blocks = size(blocks,2);  % Number of blocks

%% ESTIMATION STEP: Compute the (expected) sufficient statistics for a single
%Kalman filter sequence

% Running the Kalman filter and smoother with current parameters
% Note that log-liklihood is NOT re-estimated after the runKF step: This
% effectively gives the previous iteration's log-likelihood
% For more information on output, see runKF
[Zsmooth, Vsmooth, VVsmooth, loglik] = runKF(y, A, C, Q, R, Z_0, V_0);

%% MAXIMIZATION STEP (TRANSITION EQUATION)
% See (Banbura & Modugno, 2010) for details.

% Initialize output
A_new = A;
Q_new = Q;
V_0_new = V_0;

%%% 2A. UPDATE FACTOR PARAMETERS INDIVIDUALLY ----------------------------

for i = 1:num_blocks  % Loop for each block: factors are uncorrelated

    % SETUP INDEXING
    r_i = r(i);  % r_i = 1 if block is loaded
    rp = r_i*p;
    rp1 = sum(r(1:i-1))*ppC;
    b_subset = rp1+1:rp1+rp;  % Subset blocks: Helps for subsetting Zsmooth, Vsmooth
    t_start = rp1+1;          % Transition matrix factor idx start
    t_end = rp1+r_i*ppC;      % Transition matrix factor idx end

    % ESTIMATE FACTOR PORTION OF Q, A
    % Note: EZZ, EZZ_BB, EZZ_FB are parts of equations 6 and 8 in BM 2010

    % E[f_t*f_t' | Omega_T]
    EZZ = Zsmooth(b_subset, 2:end) * Zsmooth(b_subset, 2:end)'...
        +sum(Vsmooth(b_subset, b_subset, 2:end) ,3);

    % E[f_{t-1}*f_{t-1}' | Omega_T]
    EZZ_BB = Zsmooth(b_subset, 1:end-1)*Zsmooth(b_subset, 1:end-1)'...
            +sum(Vsmooth(b_subset, b_subset, 1:end-1), 3);

    % E[f_t*f_{t-1}' | Omega_T]
    EZZ_FB = Zsmooth(b_subset, 2:end)*Zsmooth(b_subset, 1:end-1)'...
        +sum(VVsmooth(b_subset, b_subset, :), 3);

    % Select transition matrix/covariance matrix for block i
    A_i = A(t_start:t_end, t_start:t_end);
    Q_i = Q(t_start:t_end, t_start:t_end);

    % Equation 6: Estimate VAR(p) for factor
    A_i(1:r_i,1:rp) = EZZ_FB(1:r_i,1:rp) * inv(EZZ_BB(1:rp,1:rp));

    % Equation 8: Covariance matrix of residuals of VAR
    Q_i(1:r_i,1:r_i) = (EZZ(1:r_i,1:r_i) - A_i(1:r_i,1:rp)* EZZ_FB(1:r_i,1:rp)') / T;

    % Place updated results in output matrix
    A_new(t_start:t_end, t_start:t_end) = A_i;
    Q_new(t_start:t_end, t_start:t_end) = Q_i;
    V_0_new(t_start:t_end, t_start:t_end) =...
        Vsmooth(t_start:t_end, t_start:t_end,1);
end

%%% 2B. UPDATING PARAMETERS FOR IDIOSYNCRATIC COMPONENT ------------------

rp1 = sum(r)*ppC;           % Col size of factor portion
niM = sum(i_idio(1:nM));    % Number of monthly values
t_start = rp1+1;            % Start of idiosyncratic component index
i_subset = t_start:rp1+niM; % Gives indices for monthly idiosyncratic component values

% Below 3 estimate the idiosyncratic component (for eqns 6, 8 BM 2010)

% E[f_t*f_t' | \Omega_T]
EZZ = diag(diag(Zsmooth(t_start:end, 2:end) * Zsmooth(t_start:end, 2:end)'))...
    + diag(diag(sum(Vsmooth(t_start:end, t_start:end, 2:end), 3)));

% E[f_{t-1}*f_{t-1}' | \Omega_T]
EZZ_BB = diag(diag(Zsmooth(t_start:end, 1:end-1)* Zsmooth(t_start:end, 1:end-1)'))...
       + diag(diag(sum(Vsmooth(t_start:end, t_start:end, 1:end-1), 3)));

% E[f_t*f_{t-1}' | \Omega_T]
EZZ_FB = diag(diag(Zsmooth(t_start:end, 2:end)*Zsmooth(t_start:end, 1:end-1)'))...
       + diag(diag(sum(VVsmooth(t_start:end, t_start:end, :), 3)));

A_i = EZZ_FB * diag(1./diag((EZZ_BB)));  % Equation 6
Q_i = (EZZ - A_i*EZZ_FB') / T;           % Equation 8

% Place updated results in output matrix
A_new(i_subset, i_subset) = A_i(1:niM,1:niM);
Q_new(i_subset, i_subset) = Q_i(1:niM,1:niM);
V_0_new(i_subset, i_subset) = diag(diag(Vsmooth(i_subset, i_subset, 1)));

%% 3 MAXIMIZATION STEP (observation equation)

%%% INITIALIZATION AND SETUP ----------------------------------------------
Z_0 = Zsmooth(:,1); %zeros(size(Zsmooth,1),1); %

% Set missing data series values to 0
nanY = isnan(y);
y(nanY) = 0;

% LOADINGS
C_new = C;

% Blocks
bl = unique(blocks,'rows');  % Gives unique loadings
n_bl = size(bl,1);           % Number of unique loadings

% Initialize indices: These later help with subsetting
bl_idxM = [];  % Indicator for monthly factor loadings
bl_idxQ = [];  % Indicator for quarterly factor loadings
R_con = [];    % Block diagonal matrix giving monthly-quarterly aggreg scheme
q_con = [];

% Loop through each block
for i = 1:num_blocks
    bl_idxQ = [bl_idxQ repmat(bl(:,i),1,r(i)*ppC)];
    bl_idxM = [bl_idxM repmat(bl(:,i),1,r(i)) zeros(n_bl,r(i)*(ppC-1))];
    R_con = blkdiag(R_con, kron(R_mat,eye(r(i))));
    q_con = [q_con;zeros(r(i)*size(R_mat,1),1)];
end

% Indicator for monthly/quarterly blocks in observation matrix
bl_idxM = logical(bl_idxM);
bl_idxQ = logical(bl_idxQ);

i_idio_M = i_idio(1:nM);            % Gives 1 for monthly series
n_idio_M = length(find(i_idio_M));  % Number of monthly series
c_i_idio = cumsum(i_idio);          % Cumulative number of monthly series

for i = 1:n_bl  % Loop through unique loadings (e.g. [1 0 0 0], [1 1 0 0])

    bl_i = bl(i,:);
    rs = sum(r(logical(bl_i)));                    % Total num of blocks loaded
    idx_i = find(ismember(blocks, bl_i, 'rows'));  % Indices for bl_i
    idx_iM = idx_i(idx_i<nM+1);                    % Only monthly
    n_i = length(idx_iM);                          % Number of monthly series

    % Initialize sums in equation 13 of BGR 2010
    denom = zeros(n_i*rs,n_i*rs);
    nom = zeros(n_i,rs);

    % Stores monthly indicies. These are done for input robustness
    i_idio_i = i_idio_M(idx_iM);
    i_idio_ii = c_i_idio(idx_iM);
    i_idio_ii = i_idio_ii(i_idio_i);

    %%% UPDATE MONTHLY VARIABLES: Loop through each period ----------------
    for t = 1:T
        Wt = diag(~nanY(idx_iM, t));  % Gives selection matrix (1 for nonmissing values)

        denom = denom +...  % E[f_t*t_t' | Omega_T]
                kron(Zsmooth(bl_idxM(i, :), t+1) * Zsmooth(bl_idxM(i, :), t+1)' + ...
                Vsmooth(bl_idxM(i, :), bl_idxM(i, :), t+1), Wt);

        nom = nom + ...  E[y_t*f_t' | \Omega_T]
              y(idx_iM, t) * Zsmooth(bl_idxM(i, :), t+1)' - ...
              Wt(:, i_idio_i) * (Zsmooth(rp1 + i_idio_ii, t+1) * ...
              Zsmooth(bl_idxM(i, :), t+1)' + ...
              Vsmooth(rp1 + i_idio_ii, bl_idxM(i, :), t+1));
    end

    vec_C = inv(denom)*nom(:);  % Eqn 13 BGR 2010

    % Place updated monthly results in output matrix
    C_new(idx_iM,bl_idxM(i,:)) = reshape(vec_C, n_i, rs);

   %%% UPDATE QUARTERLY VARIABLES -----------------------------------------

   idx_iQ = idx_i(idx_i > nM);  % Index for quarterly series
   rps = rs * ppC;

   % Monthly-quarterly aggregation scheme
   R_con_i = R_con(:,bl_idxQ(i,:));
   q_con_i = q_con;

   no_c = ~(any(R_con_i,2));
   R_con_i(no_c,:) = [];
   q_con_i(no_c,:) = [];

   % Loop through quarterly series in loading. This parallels monthly code
   for j = idx_iQ'
       % Initialization
       denom = zeros(rps,rps);
       nom = zeros(1,rps);

       idx_jQ = j-nM;  % Ordinal position of quarterly variable
       % Loc of factor structure corresponding to quarterly var residuals
       i_idio_jQ = (rp1 + n_idio_M + 5*(idx_jQ-1)+1:rp1+ n_idio_M + 5*idx_jQ);

       % Place quarterly values in output matrix
       V_0_new(i_idio_jQ, i_idio_jQ) = Vsmooth(i_idio_jQ, i_idio_jQ,1);
       A_new(i_idio_jQ(1), i_idio_jQ(1)) = A_i(i_idio_jQ(1)-rp1, i_idio_jQ(1)-rp1);
       Q_new(i_idio_jQ(1), i_idio_jQ(1)) = Q_i(i_idio_jQ(1)-rp1, i_idio_jQ(1)-rp1);

       for t=1:T
           Wt = diag(~nanY(j,t));  % Selection matrix for quarterly values

           % Intermediate steps in BGR equation 13
           denom = denom + ...
                   kron(Zsmooth(bl_idxQ(i,:), t+1) * Zsmooth(bl_idxQ(i,:), t+1)'...
                 + Vsmooth(bl_idxQ(i,:), bl_idxQ(i,:), t+1), Wt);
           nom = nom + y(j, t)*Zsmooth(bl_idxQ(i,:), t+1)';
           nom = nom -...
                Wt * ([1 2 3 2 1] * Zsmooth(i_idio_jQ,t+1) * ...
                Zsmooth(bl_idxQ(i,:),t+1)'+...
                [1 2 3 2 1]*Vsmooth(i_idio_jQ,bl_idxQ(i,:),t+1));
       end

        C_i = inv(denom) * nom';
        C_i_constr = C_i - ...  % BGR equation 13
                     inv(denom) * R_con_i'*inv(R_con_i*inv(denom)*R_con_i') * (R_con_i*C_i-q_con_i);

        % Place updated values in output structure
        C_new(j,bl_idxQ(i,:)) = C_i_constr;
   end
end

%%% 3B. UPDATE COVARIANCE OF RESIDUALS FOR OBSERVATION EQUATION -----------
% Initialize covariance of residuals of observation equation
R_new = zeros(n,n);
for t=1:T
    Wt = diag(~nanY(:,t));  % Selection matrix
    R_new = R_new + (y(:,t) - ...  % BGR equation 15
            Wt * C_new * Zsmooth(:, t+1)) * (y(:,t) - Wt*C_new*Zsmooth(:,t+1))'...
           + Wt*C_new*Vsmooth(:,:,t+1)*C_new'*Wt + (eye(n)-Wt)*R*(eye(n)-Wt);
end

R_new = R_new/T;
RR = diag(R_new); %RR(RR<1e-2) = 1e-2;
RR(i_idio_M) = 1e-04;  % Ensure non-zero measurement error. See Doz, Giannone, Reichlin (2012) for reference.
RR(nM+1:end) = 1e-04;
R_new = diag(RR);

end

%--------------------------------------------------------------------------

function [converged, decrease] = em_converged(loglik, previous_loglik, threshold, check_decreased)
%em_converged    checks whether EM algorithm has converged
%
%  Syntax:
%    [converged, decrease] = em_converged(loglik, previous_loglik, threshold, check_increased)
%
%  Description:
%    em_converged() checks whether EM has converged. Convergence occurs if
%    the slope of the log-likelihood function falls below 'threshold'(i.e.
%    f(t) - f(t-1)| / avg < threshold) where avg = (|f(t)| + |f(t-1)|)/2
%    and f(t) is log lik at iteration t. 'threshold' defaults to 1e-4.
%
%    This stopping criterion is from Numerical Recipes in C (pg. 423).
%    With MAP estimation (using priors), the likelihood can decrease
%    even if the mode of the posterior increases.
%
%  Input arguments:
%    loglik: Log-likelihood from current EM iteration
%    previous_loglik: Log-likelihood from previous EM iteration
%    threshold: Convergence threshhold. The default is 1e-4.
%    check_decreased: Returns text output if log-likelihood decreases.
%
%  Output:
%    converged (numeric): Returns 1 if convergence criteria satisfied, and 0 otherwise.
%    decrease (numeric): Returns 1 if loglikelihood decreased.

%% Instantiate variables

% Threshhold arguments: Checks default behavior
if nargin < 3, threshold = 1e-4; end
if nargin < 4, check_decreased = 1; end

% Initialize output
converged = 0;
decrease = 0;

%% Check if log-likelihood decreases (optional)

if check_decreased
    if loglik - previous_loglik < -1e-3 % allow for a little imprecision
        fprintf(1, '******likelihood decreased from %6.4f to %6.4f!\n', previous_loglik, loglik);
        decrease = 1;
    end
end

%% Check convergence criteria

delta_loglik = abs(loglik - previous_loglik);  % Difference in loglik
avg_loglik = (abs(loglik) + abs(previous_loglik) + eps)/2;

if (delta_loglik / avg_loglik) < threshold,
    converged = 1;  % Check convergence
end

end

%--------------------------------------------------------------------------

%InitCond()      Calculates initial conditions for parameter estimation
%
%  Description:
%    Given standardized data and model information, InitCond() creates
%    initial parameter estimates. These are intial inputs in the EM
%    algorithm, which re-estimates these parameters using Kalman filtering
%    techniques.
%
%Inputs:
%  - x:      Standardized data
%  - r:      Number of common factors for each block
%  - p:      Number of lags in transition equation
%  - blocks: Gives series loadings
%  - optNaN: Option for missing values in spline. See remNaNs_spline() for details.
%  - Rcon:   Incorporates estimation for quarterly series (i.e. "tent structure")
%  - q:      Constraints on loadings for quarterly variables
%  - NQ:     Number of quarterly variables
%  - i_idio: Logical. Gives index for monthly variables (1) and quarterly (0)
%
%Output:
%  - A:   Transition matrix
%  - C:   Observation matrix
%  - Q:   Covariance for transition equation residuals
%  - R:   Covariance for observation equation residuals
%  - Z_0: Initial value of state
%  - V_0: Initial value of covariance matrix

function [ A, C, Q, R, Z_0, V_0] = InitCond(x,r,p,blocks,optNaN,Rcon,q,nQ,i_idio)

pC = size(Rcon,2);  % Gives 'tent' structure size (quarterly to monthly)
ppC = max(p,pC);
n_b = size(blocks,2);  % Number of blocks

OPTS.disp=0;  % Turns off diagnostic information for eigenvalue computation
[xBal,indNaN] = remNaNs_spline(x,optNaN);  % Spline without NaNs

[T,N] = size(xBal);  % Time T series number N
nM = N-nQ;           % Number of monthly series

xNaN = xBal;
xNaN(indNaN) = nan;  % Set missing values equal to NaNs
res = xBal;          % Spline output equal to res: Later this is used for residuals
resNaN = xNaN;       % Later used for residuals

% Initialize model coefficient output
C = [];
A = [];
Q = [];
V_0 = [];

% Set the first observations as NaNs: For quarterly-monthly aggreg. scheme
indNaN(1:pC-1,:) = true;

for i = 1:n_b  % Loop for each block

    r_i = r(i);  % r_i = 1 when block is loaded

    %% Observation equation -----------------------------------------------

    C_i = zeros(N,r_i*ppC);     % Initialize state variable matrix helper
    idx_i = find(blocks(:,i));  % Returns series index loading block i
    idx_iM = idx_i(idx_i<nM+1); % Monthly series indicies for loaded blocks
    idx_iQ = idx_i(idx_i>nM);   % Quarterly series indicies for loaded blocks

    % Returns eigenvector v w/largest eigenvalue d
    [v, d] = eigs(cov(res(:,idx_iM)), r_i, 'lm');

    % Flip sign for cleaner output. Gives equivalent results without this section
    if(sum(v) < 0)
        v = -v;
    end

    % For monthly series with loaded blocks (rows), replace with eigenvector
    % This gives the loading
    C_i(idx_iM,1:r_i) = v;
    f = res(:,idx_iM)*v;  % Data projection for eigenvector direction
    F = [];

    % Lag matrix using loading. This is later used for quarterly series
    for kk = 0:max(p+1,pC)-1
        F = [F f(pC-kk:end-kk,:)];
    end

    Rcon_i = kron(Rcon,eye(r_i));  % Quarterly-monthly aggregation scheme
    q_i = kron(q,zeros(r_i,1));

    % Produces projected data with lag structure (so pC-1 fewer entries)
    ff = F(:, 1:r_i*pC);

    for j = idx_iQ'      % Loop for quarterly variables

        % For series j, values are dropped to accommodate lag structure
        xx_j = resNaN(pC:end,j);

        if sum(~isnan(xx_j)) < size(ff,2)+2
            xx_j = res(pC:end,j);  % Replaces xx_j with spline if too many NaNs

        end

        ff_j = ff(~isnan(xx_j),:);
        xx_j = xx_j(~isnan(xx_j));

        iff_j = inv(ff_j'*ff_j);
        Cc = iff_j*ff_j'*xx_j;  % Least squares

        % Spline data monthly to quarterly conversion
        Cc = Cc - iff_j*Rcon_i'*inv(Rcon_i*iff_j*Rcon_i')*(Rcon_i*Cc-q_i);

        C_i(j,1:pC*r_i)=Cc';  % Place in output matrix
    end

    ff = [zeros(pC-1,pC*r_i);ff];  % Zeros in first pC-1 entries (replace dropped from lag)

    % Residual calculations
    res = res - ff*C_i';
    resNaN = res;
    resNaN(indNaN) = nan;

    C = [C C_i];  % Combine past loadings together

    %% Transition equation ------------------------------------------------

    z = F(:,1:r_i);            % Projected data (no lag)
    Z = F(:,r_i+1:r_i*(p+1));  % Data with lag 1

    A_i = zeros(r_i*ppC,r_i*ppC)';  % Initialize transition matrix

    A_temp = inv(Z'*Z)*Z'*z;  % OLS: gives coefficient value AR(p) process
    A_i(1:r_i,1:r_i*p) = A_temp';
    A_i(r_i+1:end,1:r_i*(ppC-1)) = eye(r_i*(ppC-1));

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Q_i = zeros(ppC*r_i,ppC*r_i);
    e = z  - Z*A_temp;         % VAR residuals
    Q_i(1:r_i,1:r_i) = cov(e); % VAR covariance matrix

    initV_i = reshape(inv(eye((r_i*ppC)^2)-kron(A_i,A_i))*Q_i(:),r_i*ppC,r_i*ppC);

    % Gives top left block for the transition matrix
    A = blkdiag(A,A_i);
    Q = blkdiag(Q,Q_i);
    V_0 = blkdiag(V_0,initV_i);
end

eyeN = eye(N);  % Used inside observation matrix
eyeN(:,~i_idio) = [];

C=[C eyeN];
C = [C [zeros(nM,5*nQ); kron(eye(nQ),[1 2 3 2 1])]];  % Monthly-quarterly agreggation scheme
R = diag(var(resNaN,'omitnan'));  % Initialize covariance matrix for transition matrix

ii_idio = find(i_idio);    % Indicies for monthly variables
n_idio = length(ii_idio);  % Number of monthly variables
BM = zeros(n_idio);        % Initialize monthly transition matrix values
SM = zeros(n_idio);        % Initialize monthly residual covariance matrix values

for i = 1:n_idio;  % Loop for monthly variables
    % Set observation equation residual covariance matrix diagonal
    R(ii_idio(i),ii_idio(i)) = 1e-04;

    % Subsetting series residuals for series i
    res_i = resNaN(:,ii_idio(i));

    % Returns number of leading/ending zeros
    leadZero = max( find( (1:T)' == cumsum(isnan(res_i)) ) );
    endZero  = max( find( (1:T)' == cumsum(isnan(res_i(end:-1:1))) ) );

    % Truncate leading and ending zeros
    res_i = res(:,ii_idio(i));
    res_i(end-endZero + 1:end) = [];
    res_i(1:leadZero) = [];

    % Linear regression: AR 1 process for monthly series residuals
    BM(i,i) = inv(res_i(1:end-1)'*res_i(1:end-1))*res_i(1:end-1)'*res_i(2:end,:);
    SM(i,i) = cov(res_i(2:end)-res_i(1:end-1)*BM(i,i));  % Residual covariance matrix

end

Rdiag = diag(R);
sig_e = Rdiag(nM+1:N)/19;
Rdiag(nM+1:N) = 1e-04;
R = diag(Rdiag);  % Covariance for obs matrix residuals

% For BQ, SQ
rho0 = 0.1;
temp = zeros(5);
temp(1,1) = 1;

% Blocks for covariance matrices
SQ = kron(diag((1-rho0^2)*sig_e),temp);
BQ = kron(eye(nQ),[[rho0 zeros(1,4)];[eye(4),zeros(4,1)]]);
initViQ = reshape(inv(eye((5*nQ)^2)-kron(BQ,BQ))*SQ(:),5*nQ,5*nQ);
initViM = diag(1./diag(eye(size(BM,1))-BM.^2)).*SM;

% Output
A = blkdiag(A, BM, BQ);                % Observation matrix
Q = blkdiag(Q, SM, SQ);                % Residual covariance matrix (transition)
Z_0 = zeros(size(A,1),1);              % States
V_0 = blkdiag(V_0, initViM, initViQ);  % Covariance of states

end

function [zsmooth, Vsmooth, VVsmooth, loglik] = runKF(Y, A, C, Q, R, Z_0, V_0);
%runKF()    Applies Kalman filter and fixed-interval smoother
%
%  Syntax:
%    [zsmooth, Vsmooth, VVsmooth, loglik] = runKF(Y, A, C, Q, R, Z_0, V_0)
%
%  Description:
%    runKF() applies a Kalman filter and fixed-interval smoother. The
%    script uses the following model:
%           Y_t = C_t Z_t + e_t for e_t ~ N(0, R)
%           Z_t = A Z_{t-1} + mu_t for mu_t ~ N(0, Q)

%  Throughout this file:
%    'm' denotes the number of elements in the state vector Z_t.
%    'k' denotes the number of elements (observed variables) in Y_t.
%    'nobs' denotes the number of time periods for which data are observed.
%
%  Input parameters:
%    Y: k-by-nobs matrix of input data
%    A: m-by-m transition matrix
%    C: k-by-m observation matrix
%    Q: m-by-m covariance matrix for transition equation residuals (mu_t)
%    R: k-by-k covariance for observation matrix residuals (e_t)
%    Z_0: 1-by-m vector, initial value of state
%    V_0: m-by-m matrix, initial value of state covariance matrix
%
%  Output parameters:
%    zsmooth: k-by-(nobs+1) matrix, smoothed factor estimates
%             (i.e. zsmooth(:,t+1) = Z_t|T)
%    Vsmooth: k-by-k-by-(nobs+1) array, smoothed factor covariance matrices
%             (i.e. Vsmooth(:,:,t+1) = Cov(Z_t|T))
%    VVsmooth: k-by-k-by-nobs array, lag 1 factor covariance matrices
%              (i.e. Cov(Z_t,Z_t-1|T))
%    loglik: scalar, log-likelihood
%
%  References:
%  - QuantEcon's "A First Look at the Kalman Filter"
%  - Adapted from replication files for:
%    "Nowcasting", 2010, (by Marta Banbura, Domenico Giannone and Lucrezia
%    Reichlin), in Michael P. Clements and David F. Hendry, editors, Oxford
%    Handbook on Economic Forecasting.
%
% The software can be freely used in applications.
% Users are kindly requested to add acknowledgements to published work and
% to cite the above reference in any resulting publications

S = SKF(Y, A, C, Q, R, Z_0, V_0);  % Kalman filter
S = FIS(A, S);                     % Fixed interval smoother

% Organize output
zsmooth = S.ZmT;
Vsmooth = S.VmT;
VVsmooth = S.VmT_1;
loglik = S.loglik;

end

%______________________________________________________________________
function S = SKF(Y, A, C, Q, R, Z_0, V_0)
%SKF    Applies Kalman filter
%
%  Syntax:
%    S = SKF(Y, A, C, Q, R, Z_0, V_0)
%
%  Description:
%    SKF() applies the Kalman filter

%  Input parameters:
%    Y: k-by-nobs matrix of input data
%    A: m-by-m transition matrix
%    C: k-by-m observation matrix
%    Q: m-by-m covariance matrix for transition equation residuals (mu_t)
%    R: k-by-k covariance for observation matrix residuals (e_t)
%    Z_0: 1-by-m vector, initial value of state
%    V_0: m-by-m matrix, initial value of state covariance matrix
%
%  Output parameters:
%    S.Zm: m-by-nobs matrix, prior/predicted factor state vector
%          (S.Zm(:,t) = Z_t|t-1)
%    S.ZmU: m-by-(nobs+1) matrix, posterior/updated state vector
%           (S.Zm(t+1) = Z_t|t)
%    S.Vm: m-by-m-by-nobs array, prior/predicted covariance of factor
%          state vector (S.Vm(:,:,t) = V_t|t-1)
%    S.VmU: m-by-m-by-(nobs+1) array, posterior/updated covariance of
%           factor state vector (S.VmU(:,:,t+1) = V_t|t)
%    S.loglik: scalar, value of likelihood function
%    S.k_t: k-by-m Kalman gain

%% INITIALIZE OUTPUT VALUES ---------------------------------------------
  % Output structure & dimensions of state space matrix
  [~, m] = size(C);

  % Outputs time for data matrix. "number of observations"
  nobs  = size(Y,2);

  % Instantiate output
  S.Zm  = nan(m, nobs);       % Z_t | t-1 (prior)
  S.Vm  = nan(m, m, nobs);    % V_t | t-1 (prior)
  S.ZmU = nan(m, nobs+1);     % Z_t | t (posterior/updated)
  S.VmU = nan(m, m, nobs+1);  % V_t | t (posterior/updated)
  S.loglik = 0;

%% SET INITIAL VALUES ----------------------------------------------------
  Zu = Z_0;  % Z_0|0 (In below loop, Zu gives Z_t | t)
  Vu = V_0;  % V_0|0 (In below loop, Vu guvse V_t | t)

  % Store initial values
  S.ZmU(:,1)    = Zu;
  S.VmU(:,:,1)  = Vu;

%% KALMAN FILTER PROCEDURE ----------------------------------------------
  for t = 1:nobs
      %%% CALCULATING PRIOR DISTIBUTION----------------------------------

      % Use transition eqn to create prior estimate for factor
      % i.e. Z = Z_t|t-1
      Z   = A * Zu;

      % Prior covariance matrix of Z (i.e. V = V_t|t-1)
      %   Var(Z) = Var(A*Z + u_t) = Var(A*Z) + Var(\epsilon) =
      %   A*Vu*A' + Q
      V   = A * Vu* A' + Q;
      V   =  0.5 * (V+V');  % Trick to make symmetric

      %%% CALCULATING POSTERIOR DISTRIBUTION ----------------------------

      % Removes missing series: These are removed from Y, C, and R
      [Y_t, C_t, R_t, ~] = MissData(Y(:,t), C, R);

      % Check if y_t contains no data. If so, replace Zu and Vu with prior.
      if isempty(Y_t)
          Zu = Z;
          Vu = V;
      else
          % Steps for variance and population regression coefficients:
          % Var(c_t*Z_t + e_t) = c_t Var(A) c_t' + Var(u) = c_t*V *c_t' + R
          VC  = V * C_t';
          iF  = inv(C_t * VC + R_t);

          % Matrix of population regression coefficients (QuantEcon eqn #4)
          VCF = VC*iF;

          % Gives difference between actual and predicted observation
          % matrix values
          innov  = Y_t - C_t*Z;

          % Update estimate of factor values (posterior)
          Zu  = Z  + VCF * innov;

          % Update covariance matrix (posterior) for time t
          Vu  = V  - VCF * VC';
          Vu   =  0.5 * (Vu+Vu'); % Approximation trick to make symmetric

          % Update log likelihood
          S.loglik = S.loglik + 0.5*(log(det(iF))  - innov'*iF*innov);
      end

      %%% STORE OUTPUT----------------------------------------------------

      % Store covariance and observation values for t-1 (priors)
      S.Zm(:,t)   = Z;
      S.Vm(:,:,t) = V;

      % Store covariance and state values for t (posteriors)
      % i.e. Zu = Z_t|t   & Vu = V_t|t
      S.ZmU(:,t+1)    = Zu;
      S.VmU(:,:,t+1)  = Vu;
  end

  % Store Kalman gain k_t
  if isempty(Y_t)
      S.k_t = zeros(m,m);
  else
      S.k_t = VCF * C_t;
  end

end

%______________________________________________________________________
function S = FIS(A, S)
%FIS()    Applies fixed-interval smoother
%
%  Syntax:
%    S = FIS(A, S)
%
%  Description:
%    SKF() applies a fixed-interval smoother, and is used in conjunction
%    with SKF(). See  page 154 of 'Forecasting, structural time series models
%    and the Kalman filter' for more details (Harvey, 1990).
%
%  Input parameters:
%    A: m-by-m transition matrix
%    S: structure returned by SKF()
%
%  Output parameters:
%    S: FIS() adds the following smoothed estimates to the S structure:
%    - S.ZmT: m-by-(nobs+1) matrix, smoothed states
%             (S.ZmT(:,t+1) = Z_t|T)
%    - S.VmT: m-by-m-by-(nobs+1) array, smoothed factor covariance
%             matrices (S.VmT(:,:,t+1) = V_t|T = Cov(Z_t|T))
%    - S.VmT_1: m-by-m-by-nobs array, smoothed lag 1 factor covariance
%               matrices (S.VmT_1(:,:,t) = Cov(Z_t Z_t-1|T))
%
%  Model:
%   Y_t = C_t Z_t + e_t for e_t ~ N(0, R)
%   Z_t = A Z_{t-1} + mu_t for mu_t ~ N(0, Q)

%% ORGANIZE INPUT ---------------------------------------------------------

% Initialize output matrices
  [m, nobs] = size(S.Zm);
  S.ZmT = zeros(m,nobs+1);
  S.VmT = zeros(m,m,nobs+1);

  % Fill the final period of ZmT, VmT with SKF() posterior values
  S.ZmT(:,nobs+1)   = squeeze(S.ZmU(:, nobs+1));
  S.VmT(:,:,nobs+1) = squeeze(S.VmU(:,:, nobs+1));

  % Initialize VmT_1 lag 1 covariance matrix for final period
  S.VmT_1(:,:,nobs) = (eye(m)-S.k_t) *A*squeeze(S.VmU(:,:,nobs));

  % Used for recursion process. See companion file for details
  J_2 = squeeze(S.VmU(:,:,nobs)) * A' * pinv(squeeze(S.Vm(:,:,nobs)));

  %% RUN SMOOTHING ALGORITHM ----------------------------------------------

  % Loop through time reverse-chronologically (starting at final period nobs)
    for t = nobs:-1:1

        % Store posterior and prior factor covariance values
        VmU = squeeze(S.VmU(:,:,t));
        Vm1 = squeeze(S.Vm(:,:,t));

        % Store previous period smoothed factor covariance and lag-1 covariance
        V_T = squeeze(S.VmT(:,:,t+1));
        V_T1 = squeeze(S.VmT_1(:,:,t));

        J_1 = J_2;

        % Update smoothed factor estimate
        S.ZmT(:,t) = S.ZmU(:,t) + J_1 * (S.ZmT(:,t+1) - A * S.ZmU(:,t)) ;

        % Update smoothed factor covariance matrix
        S.VmT(:,:,t) = VmU + J_1 * (V_T - Vm1) * J_1';

        if t>1
            % Update weight
            J_2 = squeeze(S.VmU(:, :, t-1)) * A' * pinv(squeeze(S.Vm(:,:,t-1)));

            % Update lag 1 factor covariance matrix
            S.VmT_1(:,:,t-1) = VmU * J_2'+J_1 * (V_T1 - A * VmU) * J_2';
        end
    end

end

function [y,C,R,L]  = MissData(y,C,R)
% Syntax:
% Description:
%   Eliminates the rows in y & matrices C, R that correspond to missing
%   data (NaN) in y
%
% Input:
%   y: Vector of observations at time t
%   C: Observation matrix
%   R: Covariance for observation matrix residuals
%
% Output:
%   y: Vector of observations at time t (reduced)
%   C: Observation matrix (reduced)
%   R: Covariance for observation matrix residuals
%   L: Used to restore standard dimensions(n x #) where # is the nr of
%      available data in y

  % Returns 1 for nonmissing series
  ix = ~isnan(y);

  % Index for columns with nonmissing variables
  e  = eye(size(y,1));
  L  = e(:,ix);

  % Removes missing series
  y  = y(ix);

  % Removes missing series from observation matrix
  C  =  C(ix,:);

  % Removes missing series from transition matrix
  R  =  R(ix,ix);

end

author: FRBNY-TimeSeriesAnalysis

repo: Nowcasting


%% Etract IDE features
clear;clc;
addpath('../caffe/matlab/');
addpath(genpath('utils/'));
% load model and creat network
caffe.set_device(0);
caffe.set_mode_gpu();
netname = 'ResNet_50'; % network: CaffeNet  or ResNet_50

% set your path to the prototxt and model
model =  ['../models/market/' netname '/' netname '_test.prototxt'];
weights = ['../output/market_train/IDE_' netname '.caffemodel'];
net = caffe.Net(model, weights, 'test');

if strcmp(netname, 'CaffeNet')
    im_size = 227;
    feat_dim = 4096;
else
    im_size = 224;
    feat_dim = 2048;
end

% mean data
mean_data = importdata('../caffe/matlab/+caffe/imagenet/ilsvrc_2012_mean.mat');
image_mean = mean_data;
off = floor((size(image_mean,1) - im_size)/2)+1;
image_mean = image_mean(off:off+im_size-1, off:off+im_size-1, :);

ef_path = {'dataset/bounding_box_train/', 'dataset/bounding_box_test/', 'dataset/query/'};
ef_name = {'train', 'test', 'query'};

if ~exist('feat')
    mkdir('feat')
end

% extract features
for i = 1:3
    img_path = ef_path{i};
    img_file = dir([img_path '*.jpg']);
    feat = single(zeros(feat_dim, length(img_file)));

    for n = 1:length(img_file)
        if mod(n, 1000) ==0
            fprintf('%s: %d/%d\n',ef_name{i}, n, length(img_file))
        end
        img_name = [img_path  img_file(n).name];
        im = imread(img_name);
        im = prepare_img( im, image_mean, im_size);
        feat_img = net.forward({im});
        feat(:, n) = single(feat_img{1}(:));
    end

    save(['feat/IDE_'  netname  '_' ef_name{i} '.mat'], 'feat');
    feat = [];
end

caffe.reset_all();


text:
function mean = binomial_mean ( a, b )

%*****************************************************************************80
%
%% BINOMIAL_MEAN returns the mean of the Binomial PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A, the number of trials.
%    1 <= A.
%
%    Input, real B, the probability of success on one trial.
%    0.0 <= B <= 1.0.
%
%    Output, real MEAN, the expected value of the number of
%    successes in A trials.
%
  mean = a * b;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pyramid_rule ( legendre_order, jacobi_order, filename )

%*****************************************************************************80
%
%% PYRAMID_RULE generates a quadrature rule for a pyramid.
%
%  Discussion:
%
%    This program computes a quadrature rule for a pyramid
%    and writes it to a file.
%
%    The user specifies:
%    * the LEGENDRE_ORDER (number of points in the X and Y dimensions)
%    * the JACOBI_ORDER (number of points in the Z dimension)
%    * FILENAME< the root name of the output files.
%
%    The integration region is:
%
%      - ( 1 - Z ) <= X <= 1 - Z
%      - ( 1 - Z ) <= Y <= 1 - Z
%                0 <= Z <= 1.
%
%    When Z is zero, the integration region is a square lying in the (X,Y)
%    plane, centered at (0,0,0) with "radius" 1.  As Z increases to 1, the
%    radius of the square diminishes, and when Z reaches 1, the square has
%    contracted to the single point (0,0,1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 July 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PYRAMID_RULE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute a quadrature rule for approximating\n' );
  fprintf ( 1, '  the integral of a function over a pyramid.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The user specifies:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEGENDRE_ORDER, the order of the Legendre rule for X and Y.\n' );
  fprintf ( 1, '  JACOBI_ORDER, the order of the Jacobi rule for Z,\n' );
  fprintf ( 1, '  FILENAME, the prefix for the three output files:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    filename_w.txt - the weight file\n' );
  fprintf ( 1, '    filename_x.txt - the abscissa file.\n' );
  fprintf ( 1, '    filename_r.txt - the region file.\n' );
%
%  Get the Legendre order.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    legendre_order = input ( '  Enter the Legendre rule order:' );
  elseif ( ischar ( legendre_order ) )
    legendre_order = str2num ( legendre_order );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The requested Legendre order of the rule is %d.\n', ...
    legendre_order );
%
%  Get the Jacobi order.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    jacobi_order = input ( '  Enter the Jacobi rule order:' );
  elseif ( ischar ( jacobi_order ) )
    jacobi_order = str2num ( jacobi_order );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The requested Jacobi order of the rule is %d.\n', ...
    jacobi_order );
%
%  Get the output option or quadrature file root name:
%
  if ( nargin < 3 )
    fprintf ( 1, '\n' );
    filename = input ( '  Enter the "root name" of the quadrature files).' );
  end

  pyramid_handle ( legendre_order, jacobi_order, filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PYRAMID_RULE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ x, w ] = jacobi_compute ( order, alpha, beta )

%*****************************************************************************80
%
%% JACOBI_COMPUTE computes the abscissa and weights for Jacobi quadrature.
%
%  Discussion:
%
%    The integration interval is [ -1, 1 ].
%
%    The weight function is w(x) = (1-X)^ALPHA * (1+X)^BETA.
%
%    The integral to approximate is:
%
%      Integral ( -1 <= X <= 1 ) (1-X)^ALPHA * (1+X)^BETA * F(X) dX
%
%    The quadrature formula is:
%
%      Sum ( 1 <= I <= ORDER ) W(I) * F ( X(I) )
%
%    Thanks to Xu Xiang of Fudan University for pointing out that
%    an earlier implementation of this routine was incorrect!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2008
%
%  Author:
%
%    Original FORTRAN77 version by Arthur Stroud, Don Secrest.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, integer ORDER, the order of the quadrature rule to be computed.
%
%    Input, real ALPHA, BETA, the exponents of (1-X) and
%    (1+X) in the quadrature rule.  For simple Legendre quadrature,
%    set ALPHA = BETA = 0.0.  -1.0 < ALPHA and -1.0 < BETA are required.
%
%    Output, real X(ORDER), the abscissas.
%
%    Output, real W(ORDER), the weights.
%
  if ( order < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'JACOBI_COMPUTE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of ORDER = %d\n', order );
    error ( 'JACOBI_COMPUTE - Fatal error!' );
  end
%
%  Check ALPHA and BETA.
%
  if ( alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'JACOBI_COMPUTE - Fatal error!\n' );
    fprintf ( 1, '  -1.0 < ALPHA is required.\n' );
    error ( 'JACOBI_COMPUTE - Fatal error!' );
  end

  if ( beta <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'JACOBI_COMPUTE - Fatal error!\n' );
    fprintf ( 1, '  -1.0 < BETA is required.\n' );
    error ( 'JACOBI_COMPUTE - Fatal error!' );
  end
%
%  Set the recursion coefficients.
%
  b = zeros ( order, 1 );
  c = zeros ( order, 1 );
  w = zeros ( order, 1 );
  x = zeros ( order, 1 );

  for i = 1 : order

    if ( alpha + beta == 0.0 || beta - alpha == 0.0 )

      b(i) = 0.0;

    else

      b(i) = ( alpha + beta ) * ( beta - alpha ) / ...
            ( ( alpha + beta + 2 * i ) ...
            * ( alpha + beta + 2 * i - 2 ) );

    end

    if ( i == 1 )

      c(i) = 0.0;

    else

      c(i) = 4.0 * ( i - 1 ) * ( alpha + i - 1 ) * ( beta + i - 1 ) ...
        * ( alpha + beta + i - 1 ) / ( ( alpha + beta + 2 * i - 1 ) ...
        * ( alpha + beta + 2 * i - 2 )^2 * ( alpha + beta + 2 * i - 3 ) );

    end

  end

  delta = r8_gamma ( alpha        + 1.0 ) ...
        * r8_gamma (         beta + 1.0 ) ...
        / r8_gamma ( alpha + beta + 2.0 );

  cc = delta * 2.0^( alpha + beta + 1.0 ) * prod ( c(2:order) );

  for i = 1 : order

    if ( i == 1 )

      an = alpha / order;
      bn = beta / order;

      r1 = ( 1.0 + alpha ) * ( 2.78 / ( 4.0 + order * order ) ...
        + 0.768 * an / order );

      r2 = 1.0 + 1.48 * an + 0.96 * bn + 0.452 * an^2 + 0.83 * an * bn;

      x0 = ( r2 - r1 ) / r2;

    elseif ( i == 2 )

      r1 = ( 4.1 + alpha ) / ...
        ( ( 1.0 + alpha ) * ( 1.0 + 0.156 * alpha ) );

      r2 = 1.0 + 0.06 * ( order - 8.0 ) * ( 1.0 + 0.12 * alpha ) / order;

      r3 = 1.0 + 0.012 * beta * ...
        ( 1.0 + 0.25 * abs ( alpha ) ) / order;

      x0 = x0 - r1 * r2 * r3 * ( 1.0 - x0 );

    elseif ( i == 3 )

      r1 = ( 1.67 + 0.28 * alpha ) / ( 1.0 + 0.37 * alpha );

      r2 = 1.0 + 0.22 * ( order - 8.0 ) / order;

      r3 = 1.0 + 8.0 * beta / ( ( 6.28 + beta ) * order * order );

      x0 = x0 - r1 * r2 * r3 * ( x(1) - x0 );

    elseif ( i < order - 1 )

      x0 = 3.0 * x(i-1) - 3.0 * x(i-2) + x(i-3);

    elseif ( i == order - 1 )

      r1 = ( 1.0 + 0.235 * beta ) / ( 0.766 + 0.119 * beta );

      r2 = 1.0 / ( 1.0 + 0.639 * ( order - 4.0 ) ...
        / ( 1.0 + 0.71 * ( order - 4.0 ) ) );

      r3 = 1.0 / ( 1.0 + 20.0 * alpha / ( ( 7.5 + alpha ) * order * order ) );

      x0 = x0 + r1 * r2 * r3 * ( x0 - x(i-2) );

    elseif ( i == order )

      r1 = ( 1.0 + 0.37 * beta ) / ( 1.67 + 0.28 * beta );

      r2 = 1.0 / ( 1.0 + 0.22 * ( order - 8.0 ) / order );

      r3 = 1.0 / ( 1.0 + 8.0 * alpha / ( ( 6.28 + alpha ) * order * order ) );

      x0 = x0 + r1 * r2 * r3 * ( x0 - x(i-2) );

    end

    [ x0, dp2, p1 ] = jacobi_root ( x0, order, alpha, beta, b, c );

    x(i) = x0;
    w(i) = cc / ( dp2 * p1 );

  end
%
%  Reverse the order of the values.
%
  x(1:order) = x(order:-1:1);
  w(1:order) = w(order:-1:1);

  return
end
function [ p2, dp2, p1 ] = jacobi_recur ( x, order, alpha, beta, b, c )

%*****************************************************************************80
%
%% JACOBI_RECUR finds the value and derivative of a Jacobi polynomial.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2005
%
%  Author:
%
%    Original FORTRAN77 version by Arthur Stroud, Don Secrest.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, real X, the point at which polynomials are evaluated.
%
%    Input, integer ORDER, the order of the polynomial to be computed.
%
%    Input, real ALPHA, BETA, the exponents of (1+X) and
%    (1-X) in the quadrature rule.
%
%    Input, real B(ORDER), C(ORDER), the recursion
%    coefficients.
%
%    Output, real P2, the value of J(ORDER)(X).
%
%    Output, real DP2, the value of J'(ORDER)(X).
%
%    Output, real P1, the value of J(ORDER-1)(X).
%
  p1 = 1.0;
  dp1 = 0.0;

  p2 = x + ( alpha - beta ) / ( alpha + beta + 2.0 );
  dp2 = 1.0;

  for i = 2 : order

    p0 = p1;
    dp0 = dp1;

    p1 = p2;
    dp1 = dp2;

    p2 = ( x - b(i) ) * p1 - c(i) * p0;
    dp2 = ( x - b(i) ) * dp1 + p1 - c(i) * dp0;

  end

  return
end
function [ x, dp2, p1 ] = jacobi_root ( x, order, alpha, beta, b, c )

%*****************************************************************************80
%
%% JACOBI_ROOT improves an approximate root of a Jacobi polynomial.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 October 2005
%
%  Author:
%
%    Original FORTRAN77 version by Arthur Stroud, Don Secrest.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, real X, the approximate root.
%
%    Input, integer ORDER, the order of the polynomial to be computed.
%
%    Input, real ALPHA, BETA, the exponents of (1+X) and
%    (1-X) in the quadrature rule.
%
%    Input, real B(ORDER), C(ORDER), the recursion coefficients.
%
%    Output, real X, the improved approximate root.
%
%    Output, real DP2, the value of J'(ORDER)(X).
%
%    Output, real P1, the value of J(ORDER-1)(X).
%
  maxstep = 10;

  eps = r8_epsilon ( );

  for i = 1 : maxstep

    [ p2, dp2, p1 ] = jacobi_recur ( x, order, alpha, beta, b, c );

    d = p2 / dp2;
    x = x - d;

    if ( abs ( d ) <= eps * ( abs ( x ) + 1.0 ) )
      return
    end

  end

  return
end
function [ x, w ] = legendre_compute ( order )

%*****************************************************************************80
%
%% LEGENDRE_COMPUTE computes a Legendre quadrature rule.
%
%  Discussion:
%
%    The integration interval is [ -1, 1 ].
%
%    The weight function is w(x) = 1.0.
%
%    The integral to approximate:
%
%      Integral ( -1 <= X <= 1 ) F(X) dX
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= ORDER ) W(I) * F ( X(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by Philip Davis, Philip Rabinowitz.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%    ORDER must be greater than 0.
%
%    Output, real X(ORDER), the abscissas of the rule.
%
%    Output, real W(ORDER), the weights of the rule.
%    The weights are positive, symmetric, and should sum to 2.
%
  if ( order < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LEGENDRE_COMPUTE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of ORDER = %d\n', order );
    error ( 'LEGENDRE_COMPUTE - Fatal error!' );
  end

  w = zeros ( order, 1 );
  x = zeros ( order, 1 );

  e1 = order * ( order + 1 );

  m = floor ( ( order + 1 ) / 2 );

  for i = 1 : floor ( ( order + 1 ) / 2 )

    mp1mi = m + 1 - i;

    t = ( 4 * i - 1 ) * pi / ( 4 * order + 2 );

    x0 = cos(t) * ( 1.0 - ( 1.0 - 1.0 / ( order ) ) / ( 8 * order * order ) );

    pkm1 = 1.0;
    pk = x0;

    for k = 2 : order
      pkp1 = 2.0 * x0 * pk - pkm1 - ( x0 * pk - pkm1 ) / k;
      pkm1 = pk;
      pk = pkp1;
    end

    d1 = order * ( pkm1 - x0 * pk );

    dpn = d1 / ( 1.0 - x0 * x0 );

    d2pn = ( 2.0 * x0 * dpn - e1 * pk ) / ( 1.0 - x0 * x0 );

    d3pn = ( 4.0 * x0 * d2pn + ( 2.0 - e1 ) * dpn ) / ( 1.0 - x0 * x0 );

    d4pn = ( 6.0 * x0 * d3pn + ( 6.0 - e1 ) * d2pn ) / ( 1.0 - x0 * x0 );

    u = pk / dpn;
    v = d2pn / dpn;
%
%  Initial approximation H:
%
    h = - u * ( 1.0 + 0.5 * u * ( v + u * ( v * v - d3pn / ( 3.0 * dpn ) ) ) );
%
%  Refine H using one step of Newton's method:
%
    p = pk + h * ( dpn + 0.5 * h * ( d2pn + h / 3.0 ...
      * ( d3pn + 0.25 * h * d4pn ) ) );

    dp = dpn + h * ( d2pn + 0.5 * h * ( d3pn + h * d4pn / 3.0 ) );

    h = h - p / dp;

    xtemp = x0 + h;

    x(mp1mi) = xtemp;

    fx = d1 - h * e1 * ( pk + 0.5 * h * ( dpn + h / 3.0 ...
      * ( d2pn + 0.25 * h * ( d3pn + 0.2 * h * d4pn ) ) ) );

    w(mp1mi) = 2.0 * ( 1.0 - xtemp * xtemp ) / fx / fx;

  end

  if ( mod ( order, 2 ) == 1 )
    x(1) = 0.0;
  end
%
%  Shift the data up.
%
  nmove = floor ( ( order + 1 ) / 2 );
  ncopy = order - nmove;

  for i = 1 : nmove
    iback = order + 1 - i;
    x(iback) = x(iback-ncopy);
    w(iback) = w(iback-ncopy);
  end
%
%  Reflect values for the negative abscissas.
%
  for i = 1 : order - nmove
    x(i) = - x(order+1-i);
    w(i) = w(order+1-i);
  end

  return
end
function pyramid_handle ( legendre_order, jacobi_order, filename )

%*****************************************************************************80
%
%% PYRAMID_HANDLE computes the requested pyramid rule and outputs it.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 July 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer LEGENDRE_ORDER, JACOBI_ORDER, the orders
%    of the component Legendre and Jacobi rules.
%
%    Input, string FILENAME, the rootname for the files,
%    write files 'file_w.txt' and 'file_x.txt', and 'file_r.txt', weights,
%    abscissas, and region.
%
  dim_num = 3;
%
%  Compute the factor rules.
%
  [ legendre_x, legendre_w ] = legendre_compute ( legendre_order );

  jacobi_alpha = 2.0;
  jacobi_beta = 0.0;

  [ jacobi_x, jacobi_w ] = jacobi_compute ( jacobi_order, jacobi_alpha, ...
    jacobi_beta );
%
%  Compute the pyramid rule.
%
  pyramid_order = legendre_order * legendre_order * jacobi_order;

  volume = 4.0 / 3.0;

  pyramid_w = zeros ( pyramid_order, 1 );
  pyramid_x = zeros ( dim_num, pyramid_order );

  l = 0;
  for k = 1 : jacobi_order
    xk = ( jacobi_x(k) + 1.0 ) / 2.0;
    wk = jacobi_w(k) / 2.0;
    for j = 1 : legendre_order
      xj = legendre_x(j);
      wj = legendre_w(j);
      for i = 1 : legendre_order
        xi = legendre_x(i);
        wi = legendre_w(i);
        l = l + 1;
        pyramid_w(l) = wi * wj * wk / 4.0 / volume;
        pyramid_x(1,l) = xi * ( 1.0 - xk );
        pyramid_x(2,l) = xj * ( 1.0 - xk );
        pyramid_x(3,l) =              xk;
      end
    end
  end

  pyramid_r(1:dim_num,1) = [ -1.0, -1.0, 0.0 ]';
  pyramid_r(1:dim_num,2) = [ +1.0, -1.0, 0.0 ]';
  pyramid_r(1:dim_num,3) = [ -1.0, +1.0, 0.0 ]';
  pyramid_r(1:dim_num,4) = [ +1.0, +1.0, 0.0 ]';
  pyramid_r(1:dim_num,5) = [  0.0,  0.0, 1.0 ]';
%
%  Write the rule to files.
%
  filename_w = strcat ( filename, '_w.txt' );
  filename_x = strcat ( filename, '_x.txt' );
  filename_r = strcat ( filename, '_r.txt' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Creating quadrature files.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  "Root" file name is   "%s".\n', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Weight file will be   "%s".\n', filename_w );
  fprintf ( 1, '  Abscissa file will be "%s".\n', filename_x );
  fprintf ( 1, '  Region file will be   "%s".\n', filename_r );

  r8mat_write ( filename_w, 1,       pyramid_order, pyramid_w' );
  r8mat_write ( filename_x, dim_num, pyramid_order, pyramid_x );
  r8mat_write ( filename_r, dim_num, 5,             pyramid_r );

  return
end
function value = r8_epsilon ( )

%*****************************************************************************80
%
%% R8_EPSILON returns the R8 roundoff unit.
%
%  Discussion:
%
%    The roundoff unit is a number R which is a power of 2 with the
%    property that, to the precision of the computer's arithmetic,
%      1 < 1 + R
%    but
%      1 = ( 1 + R / 2 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 August 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VALUE, the roundoff unit.
%
  value = eps;

  return
end
function value = r8_gamma ( x )

%*****************************************************************************80
%
%% R8_GAMMA evaluates Gamma(X) for a real argument.
%
%  Discussion:
%
%    This routine calculates the gamma function for a real argument X.
%
%    Computation is based on an algorithm outlined in reference 1.
%    The program uses rational functions that approximate the gamma
%    function to at least 20 significant decimal digits.  Coefficients
%    for the approximation over the interval (1,2) are unpublished.
%    Those for the approximation for 12 <= X are from reference 2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by William Cody, Laura Stoltz.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    William Cody,
%    An Overview of Software Development for Special Functions,
%    in Numerical Analysis Dundee, 1975,
%    edited by GA Watson,
%    Lecture Notes in Mathematics 506,
%    Springer, 1976.
%
%    John Hart, Ward Cheney, Charles Lawson, Hans Maehly,
%    Charles Mesztenyi, John Rice, Henry Thatcher,
%    Christoph Witzgall,
%    Computer Approximations,
%    Wiley, 1968,
%    LC: QA297.C64.
%
%  Parameters:
%
%    Input, real X, the argument of the function.
%
%    Output, real VALUE, the value of the function.
%

%
%  Coefficients for minimax approximation over (12, INF).
%
  c = [ ...
   -1.910444077728E-03, ...
    8.4171387781295E-04, ...
   -5.952379913043012E-04, ...
    7.93650793500350248E-04, ...
   -2.777777777777681622553E-03, ...
    8.333333333333333331554247E-02, ...
    5.7083835261E-03 ];
%
%  Mathematical constants
%
  one = 1.0;
  half = 0.5;
  twelve = 12.0;
  two = 2.0;
  zero = 0.0;
  sqrtpi = 0.9189385332046727417803297;
%
%  Machine dependent parameters
%
  xbig = 171.624E+00;
  xminin = 2.23E-308;
  eps = 2.22E-16;
  xinf = 1.79E+308;
%
%  Numerator and denominator coefficients for rational minimax
%  approximation over (1,2).
%
  p = [ ...
   -1.71618513886549492533811E+00, ...
    2.47656508055759199108314E+01, ...
   -3.79804256470945635097577E+02, ...
    6.29331155312818442661052E+02, ...
    8.66966202790413211295064E+02, ...
   -3.14512729688483675254357E+04, ...
   -3.61444134186911729807069E+04, ...
    6.64561438202405440627855E+04 ];

  q = [ ...
   -3.08402300119738975254353E+01, ...
    3.15350626979604161529144E+02, ...
   -1.01515636749021914166146E+03, ...
   -3.10777167157231109440444E+03, ...
    2.25381184209801510330112E+04, ...
    4.75584627752788110767815E+03, ...
   -1.34659959864969306392456E+05, ...
   -1.15132259675553483497211E+05 ];

  parity = 0;
  fact = one;
  n = 0;
  y = x;
%
%  Argument is negative.
%
  if ( y <= zero )

    y = - x;
    y1 = floor ( y );
    res = y - y1;

    if ( res ~= zero )

      if ( y1 ~= floor ( y1 * half ) * two )
        parity = 1;
      end

      fact = - pi / sin ( pi * res );
      y = y + one;

    else

      res = xinf;
      value = res;
      return

    end

  end
%
%  Argument is positive.
%
  if ( y < eps )
%
%  Argument < EPS.
%
    if ( xminin <= y )
      res = one / y;
    else
      res = xinf;
      value = res;
      return
    end

  elseif ( y < twelve )

    y1 = y;
%
%  0.0 < argument < 1.0.
%
    if ( y < one )

      z = y;
      y = y + one;
%
%  1.0 < argument < 12.0.
%  Reduce argument if necessary.
%
    else

      n = floor ( y ) - 1;
      y = y - n;
      z = y - one;

    end
%
%  Evaluate approximation for 1.0 < argument < 2.0.
%
    xnum = zero;
    xden = one;
    for i = 1 : 8
      xnum = ( xnum + p(i) ) * z;
      xden = xden * z + q(i);
    end

    res = xnum / xden + one;
%
%  Adjust result for case  0.0 < argument < 1.0.
%
    if ( y1 < y )

      res = res / y1;
%
%  Adjust result for case 2.0 < argument < 12.0.
%
    elseif ( y < y1 )

      for i = 1 : n
        res = res * y;
        y = y + one;
      end

    end

  else
%
%  Evaluate for 12.0 <= argument.
%
    if ( y <= xbig )

      ysq = y * y;
      sum = c(7);
      for i = 1 : 6
        sum = sum / ysq + c(i);
      end
      sum = sum / y - y + sqrtpi;
      sum = sum + ( y - half ) * log ( y );
      res = exp ( sum );

    else

      res = xinf;
      value = res;
      return

    end

  end
%
%  Final adjustments and return.
%
  if ( parity )
    res = - res;
  end

  if ( fact ~= one )
    res = fact / res;
  end

  value = res;

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function varargout = mode(varargin)
% Allows inputs to be char
ic = ischar(varargin{1});
if ic
    varargin{1} = double(varargin{1});
end
[varargout{1:nargout}] = builtin('mode', varargin{:});
if ic
    varargout{1} = char(varargout{1});
    if nargout>2
        varargout{3} = cellfun(@char, varargout{3}, 'UniformOutput', false);
    end
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
% im2pointsLocally    computes image statistics from several images
%                     and finds the projections of laser points
%
% The version for local computation on each machine
%
% The computation core is taken from  im2points.m
% computes average image and image of standard deviations
% requires read_configuration.m
% The name of the experiment has to be specified
% it determines the location of files etc ...
%
% the scripts serves as a template for a multiprocessing
% it assumes a vector of camera IDs CamIds to be known
% indexes in the CamIds are supposed to be correct
%

donefile = '.done';

addpath /home/svoboda/Work/BlueCCal/BlueCFindingPoints

% config = localconfig('BlueCHoengg')

% Read configuration from whatever is specified on command-line (via --config=FILENAME)
config = read_configuration();

STEP4STAT = 5; % step for computing average and std images, if 1 then all images taken

im.dir = config.paths.data;
im.ext = config.files.imgext;

% get the information about machine
machinename = getenv('HOST');
CamsIds = str2num(machinename(find(machinename>47 & machinename<58)));

NoCams = size(CamsIds,2);
CamsIds

% load image names
for i=1:NoCams,
  seq(i).camId = CamsIds(i);
  if seq(i).camId > -1
	seq(i).data = dir([sprintf(im.dir,seq(i).camId),sprintf(config.files.imnames,seq(i).camId),im.ext]);
	[sprintf(im.dir,seq(i).camId),sprintf(config.files.imnames,seq(i).camId),im.ext]
  else
	seq(i).data = dir([im.dir,sprintf(config.files.imnames),im.ext]);
  end
  seq(i).size = size(seq(i).data,1);
  if seq(i).size<4
	i, seq
	error('Not enough images found. Wrong image path or name pattern?');
  end
end

% Expected number of 3D points is equal to the number of frames.
% In fact, some frames might be without any calibration point

% Because of non-consistent stopping of capturing, the sequences might
% have different number of images, select the minimal value as the right one
NoPoints = min([seq.size]);

% compute the average images that will be used for finding LEDs
% if not already computed

pointsIdx = [1:STEP4STAT:NoPoints];

t = cputime;
for i=1:NoCams,
  if ~exist(sprintf(config.files.avIM,seq(i).camId)),
	disp(sprintf('The average image of the camera %d is being computed',seq(i).camId));
	avIM = zeros(size(imread([sprintf(im.dir,seq(i).camId),seq(i).data(1).name])));
	for j=pointsIdx,
	  IM = imread([sprintf(im.dir,seq(i).camId),seq(i).data(j).name]);
	  avIM = avIM + double(IM);
	end
	avIM = uint8(round(avIM./size(pointsIdx,2)));
	imwrite(avIM,sprintf(config.files.avIM,seq(i).camId));
  else	disp('Average files already exist');
  end
end
disp(sprintf('Elapsed time for average images: %d [sec]',cputime-t))
% compute the standard deviations images that will be used for finding LEDs
% if not already computed
t = cputime;
for i=1:NoCams,
  if ~exist(sprintf(config.files.stdIM,seq(i).camId)),
	avIM = double(imread(sprintf(config.files.avIM,seq(i).camId)));
	disp(sprintf('The image of standard deviations of the camera %d is being computed',seq(i).camId));
	stdIM = zeros(size(imread([sprintf(im.dir,seq(i).camId),seq(i).data(1).name])));
	for j=pointsIdx,
	  IM = imread([sprintf(im.dir,seq(i).camId),seq(i).data(j).name]);
	  stdIM = stdIM + (double(IM)-avIM).^2;
	end
	stdIM = uint8(round(sqrt(stdIM./(size(pointsIdx,2)-1))));
	imwrite(stdIM,sprintf(config.files.stdIM,seq(i).camId));
  else
	disp('Image of standard deviations already exist')
  end
end
disp(sprintf('Elapsed time for computation of images [sec]: %d',cputime-t))

% find points in the images
Ws    = [];
IdWs  = [];
Res	  = [];

IdMat = ones(NoCams,NoPoints);
% IdMat is very important for Martinec&Pajdla filling [ECCV2002]
% it is a NoCams x NoPoints matrix,
% IdMat(i,j) = 0 -> no j-th point in i-th
% IdMat(i,j) = 1 -> point successfully detected

for i=1:NoCams,
  Points = [];
  avIM  = imread(sprintf(config.files.avIM,seq(i).camId));
  stdIM	= imread(sprintf(config.files.stdIM,seq(i).camId));
  for j=1:NoPoints,
	[pos,err] = getpoint([sprintf(im.dir,seq(i).camId),seq(i).data(j).name], 0, config.imgs, avIM, stdIM);
	if err
	  IdMat(i,j) = 0;
	  Points = [Points, [NaN; NaN; NaN]];
	else
	  Points = [Points, [pos; 1]];
	end
  end
  Ws = [Ws; Points];
  Res= [Res; size(avIM,2), size(avIM,1)];
end

idx = '.';
for i=CamsIds,
  idx = sprintf('%s%02d',idx,i);
end

save([config.files.points,idx], 'Ws','-ASCII')
% save(config.files.IdPoints,'IdWs','-ASCII')
save([config.files.Res,idx], 'Res', '-ASCII')
save([config.files.IdMat,idx], 'IdMat', '-ASCII')

% write auxiliary file that is done
done=1;
save(donefile,'done','-ascii');

% exit the Matlab
% this script is to be used in the batch mode
% hence exit at the end is necessary
exit;

author: strawlab

repo: MultiCamSelfCal


text:
function [M] = spm_nwpost (M,w)
% Get posterior distribution over m,Lambda
% FORMAT [M] = spm_nwpost (M,w)
%
% M     M.prior - params of Normal-Wishart prior
% w     Multivariate data samples
%
% M     M.post - params of Normal-Wishart posterior
%
% Bernardo and Smith, Bayesian Theory, 2000 (p.441)
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_nwpost.m 6548 2015-09-11 12:39:47Z will $

mw=mean(w,2);
Sw=cov(w',1);

N=size(w,2);

prior=M.prior;
P=prior.P;
a0=prior.a;
B0=prior.B;
beta0=prior.beta;
m0=prior.m;

post.beta=beta0+N;
post.m=(beta0*m0+N*mw)/post.beta;

post.a=a0+N/2;
post.B=B0+0.5*N*Sw+0.5*(beta0*N/post.beta)*(mw-m0)*(mw-m0)';
post.a=a0+N/2;

% Quantities for predictive density (over new samples)
post.mu_w=post.m;
w_s=(post.beta/(post.beta+1))*(post.a-0.5*(P-1));
post.Lambda_w=w_s*inv(post.B);
post.v_w=2*post.a-P+1;

% Wrap up
post.P=P;
post.N=N;
M.post=post;

author: spm


text:
function neuron = fix_mat_files(root)

    %% Get path of all .mat files in root directory
    % get filenames to process
    files = dir([root, '/**/*-*_*_*_*.mat']);
    not = dir([root, '/**/*-*-*_*_*_*.mat']);
    [~,ia] = setdiff({files.name}, {not.name});
    files = files(ia);

    %% Open up each .mat file and save as a cleaned file (and format v6)
    for ff = 1:length(files)
        if ~exist([files(ff).folder, '/out.mat'], 'file')
            load([files(ff).folder, '/', files(ff).name]);
            neuron = struct(neuron);
            save([files(ff).folder, '/out.mat'],'-v6','-struct','neuron');
            clear neuron;
        end
    end

end

author: zhoupc


text:
%%%%% Representacion en esfera (radio variable, posicion (0,0,0),escribe por el "polo sur"):

function [xyz, lin]=en_esfera(P2)
    global radio;
    global precision;

    %precision=0.001;
    P=[P2(1) -P2(2) -radio];
    V=[0,0,radio]-P;
    for i=0:1:10000;
        recta(i+1,:)=P+V.*(i*precision/10);
        if abs((recta(i+1,1))^2+recta(i+1,2)^2+(recta(i+1,3))^2-radio^2) < precision
            xyz=recta(i+1,:);
            %display('exito')
            lin=recta;
            break
        end
    end
%     figure
%     hold on
%     [x,y,z] = sphere;
%     surf(x,y,z)
%    plot3(recta(:,1),recta(:,2),recta(:,3),'*')
    %hold on
    %plot3(xyz(1),xyz(2),xyz(3),'*')
end

repo: ARTE


text:
classdef NomadStopIfFeasibleEnum < matlab.mixin.SetGet
    %NomadStopIfFeasibleEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        StopIfFeasible('Stop If Feasible', 1)
        DoNotStopIfFeasible('Do Not Stop On Feasible Solution', 0)
    end

    properties
        optVal(1,1) double
        name char = '';
    end

    methods
        function obj = NomadStopIfFeasibleEnum(name, optVal)
            obj.name = name;
            obj.optVal = optVal;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('NomadStopIfFeasibleEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('NomadStopIfFeasibleEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('NomadStopIfFeasibleEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot