text:
function save_model_path = rfcn_train(conf, imdb_train, roidb_train, varargin)
% save_model_path = rfcn_train(conf, imdb_train, roidb_train, varargin)
% --------------------------------------------------------
% R-FCN implementation
% Modified from MATLAB Faster R-CNN (https://github.com/shaoqingren/faster_rcnn)
% Copyright (c) 2016, Jifeng Dai
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

%% inputs
    ip = inputParser;
    ip.addRequired('conf',                              @isstruct);
    ip.addRequired('imdb_train',                        @iscell);
    ip.addRequired('roidb_train',                       @iscell);
    ip.addParamValue('do_val',          false,          @isscalar);
    ip.addParamValue('imdb_val',        struct(),       @isstruct);
    ip.addParamValue('roidb_val',       struct(),       @isstruct);
    ip.addParamValue('val_iters',       500,            @isscalar);
    ip.addParamValue('val_interval',    5000,           @isscalar);
    ip.addParamValue('snapshot_interval',...
                                        10000,          @isscalar);
    ip.addParamValue('solver_def_file', fullfile(pwd, 'models', 'rfcn_prototxts', 'ResNet-50L_res3a', 'solver_80k120k_lr1_3.prototxt'), ...
                                                        @isstr);
    ip.addParamValue('cache_name',      'ResNet-50L_res3a', ...
                                                        @isstr);
    ip.addParamValue('caffe_version',   'Unkonwn',      @isstr);

    ip.addParamValue('offline_roidb',          false,          @isscalar);
    ip.addParamValue('output_dir',                       @isstr);
    ip.addParameter('resume_iter',          0,             @isscalar);
    ip.addParameter('visualize_interval',          0,             @isscalar);
    ip.addParameter('net_file',    '' );

    ip.parse(conf, imdb_train, roidb_train, varargin{:});
    opts = ip.Results;

%% try to find trained model
    imdbs_name = cell2mat(cellfun(@(x) x.name, imdb_train, 'UniformOutput', false));
    cache_dir = fullfile(opts.output_dir, 'output', 'rfcn_cachedir', opts.cache_name, imdbs_name);
    save_model_path = fullfile(cache_dir, 'final');
    if exist(save_model_path, 'file')
        return;
    end

%% init
    % set random seed
    prev_rng = seed_rand(conf.rng_seed);
    caffe.set_random_seed(conf.rng_seed);

    % init caffe solver
    mkdir_if_missing(cache_dir);
    caffe_log_file_base = fullfile(cache_dir, 'caffe_log');
    caffe.init_log(caffe_log_file_base);
    caffe_solver = caffe.Solver(opts.solver_def_file);

    if (opts.resume_iter == 0)
        if iscell(opts.net_file)
          arrayfun(@(i) caffe_solver.net.copy_from(opts.net_file{i}), 1:numel(opts.net_file));
        else
          caffe_solver.net.copy_from(opts.net_file);
        end
    else
        % loading solverstate, resume
        caffe_solver.net.copy_from(fullfile(cache_dir, ...
            sprintf('iter_%d', opts.resume_iter )));
    end
    if isfield(conf, 'nFrames') && strcmp(conf.input_modality, 'rgbflow')
      weights_conv1 = caffe_solver.net.layer_vec(caffe_solver.net.name2layer_index('conv1')).params(1).get_data();
      conv1_flow  = caffe_solver.net.layer_vec(caffe_solver.net.name2layer_index('conv1_flow')).params(1).get_data();
      diff =   10 - size(weights_conv1,3) ;
      conv1_flow = padarray(weights_conv1, [0 0 diff 0], 'symmetric', 'post');
      caffe_solver.net.layer_vec(caffe_solver.net.name2layer_index('conv1_flow')).params(1).set_data(conv1_flow);
    end

    % init log
    timestamp = datestr(datevec(now()), 'yyyymmdd_HHMMSS');
    mkdir_if_missing(fullfile(cache_dir, 'log'));
    log_file = fullfile(cache_dir, 'log', ['train_', timestamp, '.txt']);
    diary(log_file);

    % set gpu/cpu
    if conf.use_gpu
        caffe.set_mode_gpu();
    else
        caffe.set_mode_cpu();
    end

    disp('conf:');
    disp(conf);
    disp('opts:');
    disp(opts);

%% making tran/val data
    if opts.offline_roidb
      load('bbox_moments.mat')
      bbox_means(:,1) = bbox_means(:,1) * 0; % set horizontal mean to zero for online target flipping
    else
      fprintf('Preparing training data...');
      [image_roidb_train, bbox_means, bbox_stds] = rfcn_prepare_image_roidb(conf, opts.imdb_train, opts.roidb_train);
      fprintf('Done.\n');
    end
      if opts.do_val
      % linear sample val vids
         inds_val =  1:length(opts.imdb_val.image_ids);
         inds_val = inds_val(:,floor(linspace(1,size(inds_val,2), opts.val_iters*conf.ims_per_batch)));

         inds_val = reshape(inds_val , conf.ims_per_batch, []);
         inds_val = num2cell(inds_val, 1);
         shuffled_inds_val = inds_val;
        if isfield(conf,'sample_vid') && conf.sample_vid
            [shuffled_inds_val,shuffled_frames_val, ~] = video_generate_random_minibatch([], [], opts.imdb_val, conf);
            shuffled_frames_val = cellfun(@(x) x(1:conf.nFramesPerVid), shuffled_frames_val, 'UniformOutput' , false); % sample first frames
        end
      end
%% training
    shuffled_inds = []; shuffled_vids  = []; shuffled_frames = [];
    train_results = [];  shuffled_inds2 = [];
    val_results = [];
    iter_ = caffe_solver.iter();
    max_iter = caffe_solver.max_iter();

    p = new_parpool(2);
    parfor i=1:2
        seed_rand(conf.rng_seed);
    end
    if isfield(conf,'sample_vid') && conf.sample_vid
      [shuffled_vids,shuffled_frames, sub_db_inds] = video_generate_random_minibatch(shuffled_vids, shuffled_frames, imdb_train{1}, conf);
    else
      [shuffled_inds, sub_db_inds] = generate_random_minibatch(shuffled_inds, imdb_train{1}, conf.ims_per_batch);
    end
    if numel(imdb_train) > 1
      [shuffled_inds2, ~] = generate_random_minibatch(shuffled_inds2, imdb_train{2}, conf.ims_per_batch);
    end
    [image_roidb_train, ~, ~] = rfcn_prepare_image_roidb_offline(conf, opts.imdb_train{1}, roidb_train{1}.rois(sub_db_inds) , opts.output_dir, sub_db_inds, bbox_means, bbox_stds);
    parHandle = parfeval(p, @rfcn_get_minibatch, 1, conf, image_roidb_train);
    tic

    while (iter_ < max_iter)
      caffe_solver.net.set_phase('train');
      % generate minibatch training data
      [~, net_inputs] = fetchNext(parHandle);
      %  [net_inputs] = feval(@rfcn_get_minibatch, conf, image_roidb_train);

      if iscell(net_inputs) && min([size(net_inputs{1},1) size(net_inputs{1},2)]) > max(conf.scales)
       fprintf('Problem with getbatch...pausing\n');
       pause;
      end
      % generate minibatch training data
      if numel(imdb_train) == 1 || ~mod(iter_,2) % sample from db1
      if isfield(conf,'sample_vid') && conf.sample_vid
        [shuffled_vids,shuffled_frames, sub_db_inds] = video_generate_random_minibatch(shuffled_vids, shuffled_frames, imdb_train{1}, conf);
      else
        [shuffled_inds, sub_db_inds] = generate_random_minibatch(shuffled_inds, image_roidb_train, conf.ims_per_batch);
      end
        [image_roidb_train, ~, ~] = rfcn_prepare_image_roidb_offline(conf, opts.imdb_train{1}, roidb_train{1}.rois(sub_db_inds) , opts.output_dir, sub_db_inds, bbox_means, bbox_stds);
      else
        [shuffled_inds2, sub_db_inds2] = generate_random_minibatch(shuffled_inds2, image_roidb_train, conf.ims_per_batch);
        [image_roidb_train, ~, ~] = rfcn_prepare_image_roidb_offline(conf, opts.imdb_train{2}, roidb_train{2}.rois(sub_db_inds2) , opts.output_dir, sub_db_inds2, bbox_means, bbox_stds);
      end

      parHandle = parfeval(p, @rfcn_get_minibatch, 1, conf, image_roidb_train);

      if  ~iscell(net_inputs) %|| isempty(image_roidb_train)
        continue;
      end

      if opts.visualize_interval >0 && ~mod(iter_, opts.visualize_interval)
        vis_inputs;
      end
      caffe_solver.net.reshape_as_input(net_inputs);

      % one iter SGD update
      caffe_solver.net.set_input_data(net_inputs);
      caffe_solver.step(1);
      rst = caffe_solver.net.get_output();
      train_results = parse_rst(train_results, rst);

      % do valdiation per val_interval iterations
      if ~mod(iter_, opts.val_interval)  && iter_ >0
        if opts.do_val
            caffe_solver.net.set_phase('test');
            % use all regoins for testing
            bt_thresh_lo = conf.bg_thresh_lo; conf.bg_thresh_lo = 0;
            batch_size = conf.batch_size; conf.batch_size = -1;
            use_flipped = conf.use_flipped; conf.use_flipped = false;
            for i = 1:length(shuffled_inds_val)
                sub_db_inds = shuffled_inds_val{i};
                 if isfield(conf,'sample_vid') && conf.sample_vid
                   [shuffled_inds_val, shuffled_frames_val, sub_db_inds] = video_generate_random_minibatch(shuffled_inds_val, shuffled_frames_val, opts.imdb_val, conf);
                 end

                [image_roidb_val, ~, ~] = rfcn_prepare_image_roidb_offline(conf, opts.imdb_val, opts.roidb_val.rois(sub_db_inds) , opts.output_dir, sub_db_inds, bbox_means, bbox_stds);
                if isempty(image_roidb_val)
                  continue;
                end

                ind = arrayfun(@(x) sum(x.overlap,2) == 1, image_roidb_val, 'UniformOutput', false);
                for b = 1:numel(sub_db_inds)
                  for f = {'overlap','class','bbox_targets',}
                    image_roidb_val(b).(char(f)) = image_roidb_val(b).(char(f))(ind{b},:) ;
                  end
                end

                net_inputs = rfcn_get_minibatch(conf, image_roidb_val);
                if ~iscell(net_inputs)
                  continue;
                end
                caffe_solver.net.reshape_as_input(net_inputs);

                caffe_solver.net.forward(net_inputs);
                rst = caffe_solver.net.get_output();
                val_results = parse_rst(val_results, rst);

            end
            % use all regoins for testing
            conf.bg_thresh_lo = bt_thresh_lo;
            conf.batch_size = batch_size;
            conf.use_flipped = use_flipped ;
        end

        show_state(iter_, train_results, val_results);
        toc;tic;
        train_results = [];
        val_results = [];
        diary; diary; % flush diary
      end

      % snapshot
      if ~mod(iter_, opts.snapshot_interval) && iter_ >0
        snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, sprintf('iter_%d', iter_));
      end

      iter_ = caffe_solver.iter();

    end

    % final snapshot
    snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, sprintf('iter_%d', iter_));
    save_model_path = snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, 'final');

    diary off;
    caffe.reset_all();
    rng(prev_rng);
end

function [shuffled_inds, sub_inds] = generate_random_minibatch(shuffled_inds, image_roidb_train, ims_per_batch)

    % shuffle training data per batch
    if isempty(shuffled_inds)
        % make sure each minibatch, only has horizontal images or vertical
        % images, to save gpu memory

        if isfield(image_roidb_train, 'sizes')
          hori_image_inds = image_roidb_train.sizes(:,2) > image_roidb_train.sizes(:,1);
        else
          hori_image_inds = arrayfun(@(x) x.im_size(2) >= x.im_size(1), image_roidb_train, 'UniformOutput', true);
        end
        vert_image_inds = ~hori_image_inds;
        hori_image_inds = find(hori_image_inds);
        vert_image_inds = find(vert_image_inds);

        % random perm
        lim = floor(length(hori_image_inds) / ims_per_batch) * ims_per_batch;
        hori_image_inds = hori_image_inds(randperm(length(hori_image_inds), lim));
        lim = floor(length(vert_image_inds) / ims_per_batch) * ims_per_batch;
        vert_image_inds = vert_image_inds(randperm(length(vert_image_inds), lim));

        % combine sample for each ims_per_batch
        hori_image_inds = reshape(hori_image_inds, ims_per_batch, []);
        vert_image_inds = reshape(vert_image_inds, ims_per_batch, []);

        shuffled_inds = [hori_image_inds, vert_image_inds];
        shuffled_inds = shuffled_inds(:, randperm(size(shuffled_inds, 2)));
        shuffled_inds = repmat(shuffled_inds, 1, 10);
        shuffled_inds = num2cell(shuffled_inds, 1);
    end

    if nargout > 1
        % generate minibatch training data
        sub_inds = shuffled_inds{1};
        assert(length(sub_inds) == ims_per_batch);
        shuffled_inds(1) = [];
    end
end

function model_path = snapshot(caffe_solver, bbox_means, bbox_stds, cache_dir, file_name)
    bbox_pred_layer_name = 'rfcn_bbox';
    idx = find(~cellfun(@isempty, strfind(caffe_solver.net.layer_names, bbox_pred_layer_name)));
    bbox_pred_layer_name = caffe_solver.net.layer_names{idx};
    weights = caffe_solver.net.params(bbox_pred_layer_name, 1).get_data();
    biase = caffe_solver.net.params(bbox_pred_layer_name, 2).get_data();
    weights_back = weights;
    biase_back = biase;

    rep_time = size(weights, 4)/length(bbox_means(:));

    bbox_stds_flatten = bbox_stds';
    bbox_stds_flatten = bbox_stds_flatten(:);
    bbox_stds_flatten = repmat(bbox_stds_flatten, [1,rep_time])';
    bbox_stds_flatten = bbox_stds_flatten(:);
    bbox_stds_flatten = permute(bbox_stds_flatten, [4,3,2,1]);

    bbox_means_flatten = bbox_means';
    bbox_means_flatten = bbox_means_flatten(:);
    bbox_means_flatten = repmat(bbox_means_flatten, [1,rep_time])';
    bbox_means_flatten = bbox_means_flatten(:);
    bbox_means_flatten = permute(bbox_means_flatten, [4,3,2,1]);

    % merge bbox_means, bbox_stds into the model
    weights = bsxfun(@times, weights, bbox_stds_flatten); % weights = weights * stds;
    biase = biase .* bbox_stds_flatten(:) + bbox_means_flatten(:); % bias = bias * stds + means;

    caffe_solver.net.set_params_data(bbox_pred_layer_name, 1, weights);
    caffe_solver.net.set_params_data(bbox_pred_layer_name, 2, biase);

    model_path = fullfile(cache_dir, file_name);
    caffe_solver.net.save(model_path);
    fprintf('Saved as %s\n', model_path);

    % restore net to original state
    caffe_solver.net.set_params_data(bbox_pred_layer_name, 1, weights_back);
    caffe_solver.net.set_params_data(bbox_pred_layer_name, 2, biase_back);
end

function show_state(iter, train_results, val_results)
    fprintf('\n------------------------- Iteration %d -------------------------\n', iter);
    fprintf('Training : accuracy %.3g, loss (cls %.3g, reg %.3g)\n', ...
        mean(train_results.accuarcy.data), ...
        mean(train_results.loss_cls.data), ...
        mean(train_results.loss_bbox.data));
    if exist('val_results', 'var') && ~isempty(val_results)
        fprintf('Testing  : accuracy %.3g, loss (cls %.3g, reg %.3g)\n', ...
            mean(val_results.accuarcy.data), ...
            mean(val_results.loss_cls.data), ...
            mean(val_results.loss_bbox.data));
    end
end

author: feichtenhofer

repo: Detect-Track


text:
function [fMc, fBvalue, fBStd, fAvalue, fSigmaLow, fSigmaHi] = calc_McduebBst(magnitudes, varargin)
    % calc_McduebBst  Calculate Mc using the function b-value vs. cut-off-magnitude: Bootstrap approach
    % [fMc, fBvalue, fBStd, fAvalue, fSigmaLow, fSigmaHi] = calc_McduebBst(mCatalog, fBinning, nWindowSize, nMinNumberEvents, nSample)
    %-------------------------------------------------------------------------------------------------------
    % Calculate Mc using the function b-value vs. cut-off-magnitude: Bootstrap approach
    % Decision criterion for b and Mc: b_i-std_Bst(b_i) <= b_ave <= b_i+std_Bst(b_i)

    % Relevant reference: Cao A., Gao, S.S., Temporal variation of seismic b-values
    % beneath northeastern Japan island arc, GRL, 29, 9, 2002
    %
    % Incoming variables:
    % magnitudes         : EQ catalog magnitudes
    % fBinning         : Bin size
    % nWindowSize      : Window size
    % nMinNumberEvents : Minimum number of events
    % nSample          : Number of bootstrap samples
    %
    % Outgoing variables:
    % fMc              : Magnitude of completeness
    % fBvalue          : b-value
    % fBStd            : 2nd moment of b-value-distribution (comparable to standard deviation)
    % fAvalue          :a-value
    % fSigmaLow        : 16-percentile of b-value distribution
    % fSigmaHi         : 84-percentile of b-value distribution
    % Author: J. Woessner
    % updated: 04.06.03

    % Check input
    p=inputParser;
    p.addOptional('fBinning',0.1);
    p.addOptional('nWindowSize',5);
    p.addOptional('nMinNumberEvents',50);
    p.addOptional('nSample',100);
    p.parse(varargin{:});

    PCTL_RANGE = [16, 84];

    fBinning            = p.Results.fBinning;
    nWindowSize         = p.Results.nWindowSize;
    nMinNumberEvents    = p.Results.nMinNumberEvents;
    nSample             = p.Results.nSample;

    % Set fix values
    fMinMag = min(magnitudes);
    fMaxMag = max(magnitudes);

    % Create bootstrap samples using bootstrap matlab toolbox
    deepdim = 1; % for deepdim = 1:size(magnitudes,2)
    mMag_bstsamp = bootrsp(magnitudes,nSample);
    % Calculate b-with magnitude
    magBins = fMinMag:fBinning:fMaxMag;
    nBins = numel(magBins);

    % mBvalue = nan(nBins, 7); % contains [meanb_value, meanb_value_std, meana_value vSigma(1) vSigma(2), fStdBst]
    mean_b_value     = nan(nBins,1);
    mean_b_value_std = nan(nBins,1);
    mean_a_value     = nan(nBins,1);
    mean_mag         = nan(nBins,1);
    sigmas           = nan(nSample,2);
    std_bst          = nan(nBins,1);

    for i = 1 : nBins
        fMag = magBins(i);
        mBvalue_bst = nan(nSample, 4); % resets each pass, contains [b_value, b_value_std, a_value, magvalue]

        magMask = mMag_bstsamp >= fMag - 0.05;

        idxs        = (i-1).* nSample + (1:nSample);
        b_value     = nan(nSample,1);
        b_value_std = nan(nSample,1);
        a_value     = nan(nSample,1);

        mBvalue_bst(idxs,4) = fMag;
        for nSamp=1:nSample
            bvalIsNan = magMask(:, nSamp);
            mCat = mMag_bstsamp(bvalIsNan, nSamp);

            % Check for minimum number of events

            if numel(mCat) >= nMinNumberEvents
                [ b_value(nSamp), b_value_std(nSamp), a_value(nSamp)] =  calc_bmemag(mCat, fBinning);
            end
        end
        mBvalue_bst(idxs, 1:3) = [b_value, b_value_std, a_value];

        % Check for Nan and create output for [16 84]-percentile
        bvalIsNan   = isnan(b_value);
        bval_nonans = b_value(~bvalIsNan);

        if ~isempty(bval_nonans)
            vSigma = prctile(bval_nonans, PCTL_RANGE);
            if length(bval_nonans)==1
                vSigma = vSigma';
            end
            fStdBst = std(bval_nonans, 1, 'omitnan'); % Calculate 2nd moment

            sigmas(i,:)         = vSigma;
            std_bst(i)          = fStdBst;
        end

        mean_b_std_a_mag    = mean(mBvalue_bst,'omitnan');
        mean_b_value(i)     = mean_b_std_a_mag(1);
        mean_b_value_std(i) = mean_b_std_a_mag(2);
        mean_a_value(i)     = mean_b_std_a_mag(3);
        mean_mag(i)         = mean_b_std_a_mag(4);

        % mBvalue(i,:) = [mean(mBvalue_bst,'omitnan'), vSigma(1), vSigma(2), fStdBst]; % contains [meanb_value, meanb_value_std, meana_value, meanmag, vSigma(1), vSigma(2), fStdBst]

    end % END of FOR fMag

    % Use bootstrap percentiles to decide for Mc
    totalSteps = nBins - nWindowSize + 1;

    b_avg       = movmean(mean_b_value, nWindowSize, 'Endpoints','discard');
    keeprows    = sigmas(1:totalSteps, 1) <= b_avg & b_avg <= sigmas(1:totalSteps, 2);

    % b_avg       = b_avg(keeprows); % apply to b_avg before changing it's size

    keeprows(nBins) = false; % OK, because len(keeprows) (aka. totalsteps) is always less than nBins
    fMc         = mean_mag(keeprows);
    fBvalue     = mean_b_value(keeprows);
    fAvalue     = mean_a_value(keeprows);
    fBStd       = mean_b_value_std(keeprows);
    fSigmaLow   = sigmas(keeprows, 1);
    fSigmaHi    = sigmas(keeprows, 2);

end

author: CelsoReyes


text:
function [ hBlobs_valid, hBlobs_centroid ] = regionShapeValidation( hBlobs, H, W, min_size, aspectratio, arearatio, iscuboid )
%REGIONSHAPEVALIDATION Summary of this function goes here
%   check whether to find a cuboid or rectangle fitting based on
%   segmentation shape
numSeg = length(hBlobs);
invalidMask = false(H,W);
if iscuboid
    invalidMask(:,1) = true;   invalidMask(:,W) = true;
    invalidMask(1,:) = true;   invalidMask(H,:) = true;
    invalidMask(round(H/2), 1) = true;
    invalidMask(round(H/2), round(W/4)) = true;
    invalidMask(round(H/2), round(W/2)) = true;
    invalidMask(round(H/2), round(3*W/4)) = true;
    invalidMask(round(H/2), round(W)) = true;
    invalidMask = bwdist(invalidMask)<10;
end
SE = strel('disk', 5, 4);
SS = strel('disk', 2, 4);

hBlobs_valid = false(numSeg, 1);
hBlobs_centroid = zeros(numSeg,2);

parfor bid = 1:length(hBlobs)
%     fprintf('%d/%d\n', bid, length(hBlobs));

    if hBlobs{bid}.size<min_size
        continue;
    end
    seg = false(H,W);
    rect = hBlobs{bid}.rect;
    seg(rect(1):rect(3),rect(2):rect(4)) = hBlobs{bid}.mask;

    binaryMap = imdilate( seg, SE);
    binaryMap = imerode( binaryMap, SE);
    binaryMap = imerode( binaryMap, SS);
    binaryMap = imdilate( binaryMap, SS);
    binaryMap = imfill(binaryMap, 'hole');

    CC = bwconncomp(binaryMap);
    valid = true(CC.NumObjects,1);
    for j = 1:CC.NumObjects
        if length(CC.PixelIdxList{j})<min_size
            valid(j) = false;
        end
        if any(invalidMask(CC.PixelIdxList{j}))
            valid(j) = false;
        end
    end
    CC.NumObjects = sum(valid);
    CC.PixelIdxList = CC.PixelIdxList(valid);

    if CC.NumObjects ~= 1
        continue;
    end

    R = regionprops( CC, 'BoundingBox', 'Area', 'ConvexArea', 'Centroid');
    boundingbox = R(1).BoundingBox;
    asrt = boundingbox(4)/boundingbox(3);
    if asrt>aspectratio || asrt<1/aspectratio
        continue;
    end
    if  R(1).Area/(R(1).ConvexArea+0.0001)<arearatio
        continue;
    end
    hBlobs_valid(bid) = true;
    hBlobs_centroid(bid,:) = R(1).Centroid;
end

end

author: zouchuhang

repo: LayoutNet


text:
function [km2] = in22km2(in2)
% Convert area from square inches to square kilometers.
% Chad A. Greene 2012
km2 = in2*6.4516E-10;

author: Sable

repo: mcbench-benchmarks


text:
function obj = addTimeVariable(obj, bounds)
    % Adds time as the NLP decision variables to the problem
    %
    % Parameters:
    %  bounds: a structed data stores the boundary information of the
    %  NLP variables @type struct

    t_var = struct();
    t_var.Name = 'T';
    t_var.Dimension = 2;
    t_var.lb = zeros(2,1);
    t_var.ub = ones(2,1);
    t_var.x0 = [0;1];
    if isfield(bounds, 't0')
        if isfield(bounds.t0,'lb')
            t_var.lb(1) = bounds.t0.lb;
        end
        if isfield(bounds.t0,'ub')
            t_var.ub(1) = bounds.t0.ub;
        end
        if isfield(bounds.t0,'x0')
            t_var.x0(1) = bounds.t0.x0;
        end
    end

    if isfield(bounds, 'tf')
        if isfield(bounds.tf,'lb')
            t_var.lb(2) = bounds.tf.lb;
        end
        if isfield(bounds.tf,'ub')
            t_var.ub(2) = bounds.tf.ub;
        end
        if isfield(bounds.tf,'x0')
            t_var.x0(2) = bounds.tf.x0;
        end
    end

    % determines the nodes at which the variables to be defined.
    if obj.Options.DistributeTimeVariable
        % time variables are defined at all nodes if distributes the weightes
        obj = addVariable(obj, 'T', 'all', t_var);

        % add an equality constraint between the time variable at
        % neighboring nodes to make sure they are same
        Ti  = SymVariable('ti',[2,1]);
        Tn  = SymVariable('tn',[2,1]);
        t_cont = SymFunction('tCont',Ti-Tn,{Ti,Tn});

        % create an array of constraints structure
        t_cstr(obj.NumNode-1) = struct();
        [t_cstr.Name] = deal(t_cont.Name);
        [t_cstr.Dimension] = deal(2);
        [t_cstr.lb] = deal(0);
        [t_cstr.ub] = deal(0);
        [t_cstr.Type] = deal('Linear');
        [t_cstr.SymFun] = deal(t_cont);
        for i=1:obj.NumNode-1
            t_cstr(i).DepVariables = [obj.OptVarTable.T(i);obj.OptVarTable.T(i+1)];
        end

        % add to the NLP constraints table
        obj = addConstraint(obj,'tCont','except-last',t_cstr);
    else
        % otherwise only define at the first node
        obj = addVariable(obj, 'T', 'first', t_var);
    end

    if isfield(bounds,'duration')
        % only impose at the first node
        T  = SymVariable('t',[2,1]);
        timeDuration = SymFunction('timeDuration',flatten(T(2)-T(1)),{T});

        if isfield(bounds.duration','lb')
            lb = bounds.duration.lb;
        else
            lb = 0;
        end
        if isfield(bounds.duration','ub')
            ub = bounds.duration.ub;
        else
            ub = inf;
        end

        addNodeConstraint(obj, timeDuration, 'T', 'first', lb, ub, 'Linear');

    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function varargout = plotsparsemarkers(varargin)
%
%  [h =] plotsparsemarkers(hp, hl, markerStyle[, numberOfMarkers = 6 [,staggered = true]])
%
%  Adds only a limited number of markers to a set of lines given by
%  the handle array hp. If hl is not [] it should be the handle of
%  the legend corresponding to the lines.
%
%  markerStyle is a cell array of styles, e.g. {'o', 'x'} describing
%  the markers for the lines in hp, the length of the cell array has
%  to be the same as hp (or shorter in which case not all lines will get
%  markers).
%
%  numberOfMarkers is an optional parameter, default = 6, that sets
%  the number of markers on each line, this number must be equal or less
%  than the number of points plotted for each line.
%
%  staggered is a boolean that tries to stagger the markers so that the
%  are not at the same x value for all lines. This defaults to true and
%  is intended to avoid markers overwriting each other if the x-values
%  are identical for the lines.
%

%%
% Check input parameters
if nargin <3
    disp('[sparsemarker]: wrong number of arguments')
    help sparsemarker
    return
end

hp = varargin{1};
hl = varargin{2};
markStyle = varargin{3};

if nargin >3
    numberOfMarkers = varargin{4};
else
    numberOfMarkers = 6;
end

if nargin >4
    staggered = varargin{5};
else
    staggered = true;
end

%%
% Define marker lines
nrLines =    length(markStyle);

%Create stagger distance
if staggered
    x  = 1:nrLines;
    stag = -1.^(mod(x,2)+1).*x./nrLines./2;
else
    stag = zeros(1,nrLines);
end

for ii = 1:nrLines
    % Determine number of points in line ii
    X = get(hp(ii),'XDATA');
    Y = get(hp(ii),'YDATA');
    interval = (X(end)-X(1))./(numberOfMarkers+2);
    targetx  = stag(ii)*interval + ...  %Stagger offset
                  linspace(X(1)+interval,X(end)-interval,...
                           numberOfMarkers);
    for jj = 1:numberOfMarkers
       [~,ind] = min(abs(targetx(jj)-X));
       mx(ii,jj) = X(ind);
       my(ii,jj) = Y(ind);
    end
end

whos mx my

%%
% plot marker lines
hpp = get(hp(ii),'Parent');
set(hpp,'LineStyleOrder',markStyle); %This is a bad solution
                                     %that works
hold on
h = plot(mx',my');
hold off
%Change markers to same color as the line
for ii = 1:nrLines
    set(h(ii),'MarkerEdgeColor',get(hp(ii),'Color'));
    set(h(ii),'Marker',markStyle{ii});
    set(h(ii),'MarkerSize',12);  %These are sizes I like
    set(h(ii),'LineWidth',2);    %you can change them using
end                              %the handle that is returned

%%
% update legend
if ~isempty(hl)
    hlc = get(hl,'Children');
    for ii = 1:nrLines
        jj = (nrLines-ii)*3+1;
        set(hlc(jj),'Marker',markStyle{ii});
        set(hlc(jj),'Marker',markStyle{ii});
        set(hlc(jj),'MarkerEdgeColor',get(hp(ii),'Color'));
        set(hlc(jj),'MarkerSize',12);
    end
end

%%
% Return handle to markers
if nargout > 0
    varargout{1} = h;
end


text:
function fx = p03_fun ( n, x )

%*****************************************************************************80
%
%% P03_FUN evaluates the integrand for problem 3.
%
%  Discussion:
%
%    D&R gives "exact" value as 13.628...
%    Mathematica returns        13.440045415012575106...
%    D&R gives Laguerre(16) as   0.44996932...
%
%  Integral:
%
%    exp ( -2 ) Integral ( 2 <= x < +oo ) 1 / ( x^1.01 ) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 December 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis, Philip Rabinowitz,
%    Methods of Numerical Integration,
%    Second Edition,
%    Dover, 2007,
%    ISBN: 0486453391,
%    LC: QA299.3.D28.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the function values.
%
  fx(1:n) = exp ( - 2.0 ) * 1.0 ./ x(1:n).^1.01;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [y,deriv] = dotprod_of_functions(w,f,g)
% This is an MV2DF (see MV2DF_API_DEFINITION.readme) which
% represents the new function,
%
%    g(w) = f(w)'g(w)
%
% where f(w) and g(w) are column vectors of the same size.
%
% Here f,g are function handles to MV2DF's.

if nargin==0
    test_this();
    return;
end

function A = extractA(w)
    A = w(1:length(w)/2);
    A = A(:).';
end

function B = extractB(w)
    B = w(1+length(w)/2:end);
    B = B(:)    ;
end

if isempty(w)

    s = stack(w,f,g);
    y = mm_special(s,@(w)extractA(w),@(w)extractB(w));
    return;
end

if isa(w,'function_handle')
    f = dotprod_of_functions([],f,g);
    y = compose_mv(f,w,[]);
    return;
end

f = dotprod_of_functions([],f,g);
if nargout==1
    y = f(w);
else
    [y,deriv] = f(w);
end

end

function test_this()

m = 5;
w = randn(2*m,1);
f = subvec([],2*m,1,m);
g = subvec([],2*m,m+1,m);

h = dotprod_of_functions([],f,g);
test_MV2DF(h,w);
end

author: nesl


text:
% Author:
% - Mehrtash Harandi (mehrtash.harandi at gmail dot com)
%
% This file is provided without any warranty of
% fitness for any purpose. You can redistribute
% this file and/or modify it under the terms of
% the GNU General Public License (GPL) as published
% by the Free Software Foundation, either version 3
% of the License or (at your option) any later version.

function outS = Stein_Divergence(X1,X2)
l1 = size(X1,3);
l2 = size(X2,3);

outS = zeros(l1,l2);
for tmpC1 = 1:l1
    for tmpC2 = 1:l2
        X = X1(:,:,tmpC1);
        Y = X2(:,:,tmpC2);

        %outS(tmpC1,tmpC2) = log(det(0.5*(X+Y))) -   0.5*log(det(X*Y)); % Comment outed by HK

        S = log(det(0.5*(X+Y))) -   0.5*log(det(X*Y));

        real_flag = isreal(S);
        if real_flag
            outS(tmpC1,tmpC2) = S;
        else
            outS(tmpC1,tmpC2) = real(S);
        end

        if  (outS(tmpC1,tmpC2) < 1e-10)
            outS(tmpC1,tmpC2) = 0.0;
        end
    end
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
verts = Target.vertices;
normals = Target.normals;
faceMap = [find(verts(:,2) < 0)'; 1:numel(find(verts(:,2) < 0))]';
newVerts = verts(verts(:,2) < 0, :);
newNormals = normals(verts(:,2) < 0, :);

flatFaces = reshape(Target.faces, numel(Target.faces), 1);
flatFacesNew = zeros(size(flatFaces));
% Flattens by cols
for v=1:numel(flatFaces)
    oldVert = flatFaces(v);
    newVert = faceMap(faceMap(:,1) == oldVert, 2);
    if isempty(newVert)
        flatFacesNew(v) = NaN;
    else
        flatFacesNew(v) = newVert;
    end
    disp(v);
end
newFaces = reshape(flatFacesNew, numel(flatFacesNew)/3, 3);
newFaces = newFaces(~any(isnan(newFaces),2),:);

TargetMissing.vertices = newVerts;
TargetMissing.faces = newFaces;
TargetMissing.normals = newNormals;

save('data/faceTargetMissing.mat', 'TargetMissing')

author: charlienash

repo: nricp

save_path: github-repos/MATLAB/charlienash-nricp


text:
function linpack_s_test27 ( )

%*****************************************************************************80
%
%% TEST27 tests SSPFA and SSPSL.
%
%  Discussion:
%
%    SSPFA and SSPSL are for packed symmetric indefinite matrices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 100;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST27\n' );
  fprintf ( 1, '  For a symmetric indefinite packed matrix,\n' );
  fprintf ( 1, '  SSPFA factors the matrix,\n' );
  fprintf ( 1, '  SSPSL solves a factored linear system.\n' );
  fprintf ( 1, '  The matrix size is N = %d\n', n );
%
%  Assign values to the matrix A and the right hand side B.
%
  b(1:n-1) = 0.0;
  b(n)= n + 1;

  k = 0;
  for j = 1 : n
    for i = 1 : j
      k = k + 1;
      if ( i == j )
        a(k) = 2.0;
      elseif ( j == i+1 )
        a(k) = -1.0;
      else
        a(k) = 0.0;
      end
    end
  end
%
%  Factor the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Factor the matrix.\n' );

  [ a, ipvt, info ] = sspfa ( a, n );

  if ( info ~= 0 )
    fprintf ( 1, '  Error!  SSPFA returns INFO = %d\n', info );
    return
  end
%
%  Solve the linear system.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the linear system.\n' );

  b = sspsl ( a, n, ipvt, b );
%
%  Print the result.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The first and last five entries of the solution:\n' );
  fprintf ( 1, '  (Should be 1,2,3,4,5,...,n-1,n):\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    if ( i <= 5 | n-5 < i )
      fprintf ( 1, '  %6d  %14f\n', i, b(i) );
    end
    if ( i == 5 )
      fprintf ( 1, '  ......  ..............\n' );
    end
  end

  return
end


text:
function job = mergeInclusions(job, varargin)
% merge inclusions into grains
%
% Syntax
%   job.mergeInclusions
%
% Input
%  job - @parentGrainReconstructor
%
% Output
%  job.grains - reconstructed parent grains
%  job.mergeId -
%
% Options
%  maxSize - maximum pixel size of an inclusion to be merged
%

[job.grains, mergeId] = merge(job.grains, 'inclusions', job.isChild | job.isParent, ...
  'host', job.isParent, varargin{:});
job.mergeId = mergeId(job.mergeId); %#ok<*PROPLC>

end

author: mtex-toolbox

repo: mtex

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


text:
function graphs = makeConnectivityCircle(hmm,k,labels,...
    centergraphs,scalegraphs,partialcorr,threshold,ColorMap,figure_number)
% Plot HMM connectomes in connectivity circle format
%   (not to be used on an HMM-TDE model or an HMM-MAR model directly;
%   for this use makeSpectralConnectivityCircle)
%
% hmm: hmm struct as comes out of hmmmar
% k: which state or states to plot, e.g. 1:4. If left empty, then all of them
% labels : names of the regions
% centergraphs: whether to center the graphs according to the across-map average
% scalegraphs: whether to scale the graphs so that each voxel has variance
%       equal 1 across maps
% partialcorr: whether to use a partial correlation matrix or a correlation
%   matrix (default: 0)
% threshold: proportion threshold above which graph connections are
%       displayed (between 0 and 1, the higher the fewer displayed connections)
% ColorMap: a (nodes by 3) matrix of [r g b] triples, with colours
%           to be used for the nodes
% figure_number: number of figure to display
%
% OUTPUT:
% graph: (regions by regions by state) array with the estimated connectivity maps
%
% Notes:
% It needs the circularGraph toolbox from Matlab in the path:
%   https://github.com/paul-kassebaum-mathworks/circularGraph
%
% Diego Vidaurre (2020)

if nargin < 4 || isempty(centergraphs), centergraphs = 0; end
if nargin < 5 || isempty(scalegraphs), scalegraphs = 0; end
if nargin < 6 || isempty(partialcorr), partialcorr = 0; end
if nargin < 7 || isempty(threshold), threshold = 0.95; end
if nargin < 8 || isempty(ColorMap), ColorMap = []; end
if nargin < 9 || isempty(figure_number), figure_number = randi(100); end

do_HMM_pca = strcmpi(hmm.train.covtype,'pca');
if ~do_HMM_pca && ~strcmp(hmm.train.covtype,'full')
    error('Cannot great a brain graph because the states do not contain any functional connectivity')
end

K = length(hmm.state);
if ~isempty(k), index_k = k; else, index_k = 1:K; end
if do_HMM_pca
    ndim = size(hmm.state(1).W.Mu_W,1);
else
    if isfield(hmm.train,'A'), ndim = size(hmm.train.A,1);
    else, ndim = size(hmm.state(1).Omega.Gam_rate,1);
    end
end

if nargin < 3 || isempty(labels)
    labels = cell(ndim,1);
    for j = 1:ndim
        if ndim < 200
            labels{j} = ['Parcel ' num2str(j)];
        else
            labels{j} = num2str(j);
        end
    end
end

graphs = zeros(ndim,ndim,K);

for k = 1:K
    if partialcorr
        [~,~,~,C] = getFuncConn(hmm,k,1);
    else
        [~,C,] = getFuncConn(hmm,k,1);
    end
    C(eye(ndim)==1) = 0;
    graphs(:,:,k) = C;
end

if centergraphs
    graphs = graphs - repmat(mean(graphs,3),[1 1 K]);
end
if scalegraphs
    graphs = graphs ./ repmat(std(graphs,[],3),[1 1 K]);
end

for ik = 1:length(index_k)
    k = index_k(ik);
    C = graphs(:,:,k);
    c = C(triu(true(ndim),1)==1); c = sort(c); c = c(end-1:-1:1);
    th = c(round(length(c)*(1-threshold)));
    C(C<th) = 0;
    try
        figure(k+figure_number)
        if isempty(ColorMap)
            circularGraph(C,'Label',labels);
        else
            circularGraph(C,'Label',labels,'Colormap',ColorMap);
        end
    catch
        error('Please get the Matlab circularGraph toolbox first')
    end
end

graphs = graphs(:,:,index_k);

end
%
%
%
% function nets_netweb_2(netF,netP,sumpics,outputdir)
% % borrowed from Steve's fslnets
%
% % replicate functionality from nets_hierarchy
% grot=prctile(abs(netF(:)),99); netmatL=netF/grot; netmatH=netP/grot;
% usenet=netmatL;  usenet(usenet<0)=0;
% N=size(netmatL,1);  grot=prctile(abs(usenet(:)),99); usenet=max(min(usenet/grot,1),-1)/2;
% DD = 1:N;
% for J = 1:N, for I = 1:J-1,   y((I-1)*(N-I/2)+J-I) = 0.5 - usenet(I,J);  end; end;
% linkages=linkage(y,'ward');
% set(0,'DefaultFigureVisible','off');
% figure;[~,~,hier]=dendrogram(linkages,0,'colorthreshold',0.75);
% close;set(0,'DefaultFigureVisible','on');
% clusters=cluster(linkages,'maxclust',10)';
%
% mkdir(outputdir)
% netjs = [fileparts(which('makeConnectivityCircle')) '/netjs'];
% copyfile(netjs,outputdir) % copy javascript stuff into place
% NP=sprintf('%s/data/dataset1',outputdir);
% save(sprintf('%s/Znet1.txt',NP),'netF','-ascii');
% save(sprintf('%s/Znet2.txt',NP),'netP','-ascii');
% save(sprintf('%s/hier.txt',NP),'hier','-ascii');
% save(sprintf('%s/linkages.txt',NP),'linkages','-ascii');
% save(sprintf('%s/clusters.txt',NP),'clusters','-ascii');
% mkdir(sprintf('%s/melodic_IC_sum.sum',NP));
% for i=1:length(DD)
%   system(sprintf('/bin/cp %s.sum/%.4d.png %s/melodic_IC_sum.sum/%.4d.png',sumpics,DD(i)-1,NP,i-1));
% end
%
% end
%

author: OHBA-analysis

repo: HMM-MAR


function price = PROJ_GMDB_ComboExpos_Linear( P, Pbar, S_0, W, call, r, params_levy, params_mort, T)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Pricing Function for Gauranteed Minimum Death Benefit Options using PROJ method
%        This version assumes a comination of exponentials model of mortality (see first reference below)
% Models Supported: Levy Processes, including jump diffusions and Black-Scholes model
% Returns: price of contract
%
% Author: Zhimin Zhang  (Original Code)
%         Justin Lars Kirkby (Convert into common framework)
%
% References:  (1) Valuing Equity-Linked Death Benefits in General Exponential
%               Levy Models, J. Comput. and Appl. Math. 2019 (Z. Zhang, Y. Yong, W. Yu)
%              (2) Efficient Option Pricing By Frame Duality with The Fast
%               Fourier Transform, SIAM J. Financial Math., 2015 (J.L. Kirkby)
%
% ----------------------
% Contract/Model Params
% ----------------------
% S_0 = initial stock price (e.g. 100)
% W   = strike  (e.g. 100)
% r   = interest rate (e.g. 0.05)
% T   = time remaining until maturity (in years, e.g. T=10)
%       NOTE: set T=-1 to price a perpetual contract (no expiry)
% call  = 1 for call (else put)
% params_levy = parameters of Levy Model
% params_mort = mortality params
%
% ----------------------
% Numerical (PROJ) Params
% ----------------------
% P  = resolution parameter (increase P to use more basis elements)
% Pbar = gridwidth parameter (increase Pbar to increase the truncated density support)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set Contract / Numerical Params
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kk=log(W/S_0);

a=2^P;
hat_a=2^Pbar;

N=hat_a*a;
Delta=1/a;Delta_s=2*pi/hat_a;

% Deterime if there is a Time to expiry in contract, else it's perpetual
if nargin < 8
    T = -1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parse Death Model Params
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lambda = params_mort.lambda;  % expo rates
A = params_mort.A; % coefficients in expo combo

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Levy Parse Model Params, set model inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model = params_levy.model;

if model == 1 %BSM (Black Scholes Merton)
    sigma = params_levy.sigmaBSM;

    Psi=@(x)-1/2*sigma^2*x.^2;
    mu=r-Psi(-1i);

    phi=@(s)mu*1i*s-1/2*sigma^2*s.^2;

    if T > 0
        diffFgc=@(c)sum(A.*lambda*(mu+sigma^2*c)./(r+lambda-phi(-1i*c)).^2.*(1-exp(-(r+lambda-phi(-1i*c))*T)-T*exp(-(r+lambda-phi(-1i*c))*T).*(r+lambda-phi(-1i*c))));
    else
        diffFgc=@(c)sum(A.*lambda*(mu+sigma^2*c)./(r+lambda-phi(-1i*c)).^2);
    end

elseif model == 3 %NIG
    alpha = params_levy.alpha;
    beta = params_levy.beta;
    NIG_delta = params_levy.delta;
    sigma = params_levy.sigma;

    gamma=sqrt(alpha^2-beta^2);

    Psi=@(x)-1/2*sigma^2*x.^2-NIG_delta*(sqrt(alpha^2-(beta+1i*x).^2)-gamma);
    mu=r-Psi(-1i);

    phi=@(s)mu*1i*s-1/2*sigma^2*s.^2-NIG_delta*(sqrt(alpha^2-(beta+1i*s).^2)-gamma);

    if T > 0
        diffFgc=@(c)sum(A.*lambda*(mu+NIG_delta*sqrt(alpha^2-(beta+c).^2).*(beta+c))./(r+lambda-phi(-1i*c)).^2.*(1-exp(-(r+lambda-phi(-1i*c))*T)-T*exp(-(r+lambda-phi(-1i*c))*T).*(r+lambda-phi(-1i*c))));
    else
        diffFgc=@(c)sum(A.*lambda*(mu+NIG_delta*sqrt(alpha^2-(beta+c).^2).*(beta+c))./(r+lambda-phi(-1i*c)).^2);
    end

elseif model == 4 %MJD (Merton Jump Diffusion)
    sigma = params_levy.sigma;
    lambda_J = params_levy.lam;
    mu_J = params_levy.muj;
    sigma_J = params_levy.sigmaj;

    Psi=@(x)-sigma^2/2*x.^2+lambda_J*(exp(1i*x*mu_J-sigma_J^2/2*x.^2)-1);
    mu=r-Psi(-1i);

    phi=@(s)1i*s*mu-sigma^2/2*s.^2+lambda_J*(exp(1i*s*mu_J-sigma_J^2/2*s.^2)-1);

    if T > 0
        diffFgc=@(c)sum(A.*lambda*(mu+sigma^2*c+lambda_J*exp(mu_J*c+1/2*sigma_J^2*c).*(mu_J+sigma_J^2*c))./(r+lambda-phi(-1i*c)).^2.*(1-exp(-(r+lambda-phi(-1i*c))*T)-T*exp(-(r+lambda-phi(-1i*c))*T).*(r+lambda-phi(-1i*c))));
    else
        diffFgc=@(c)sum(A.*lambda*(mu+sigma^2*c+lambda_J*exp(mu_J*c+1/2*sigma_J^2*c).*(mu_J+sigma_J^2*c))./(r+lambda-phi(-1i*c)).^2);
    end

elseif model == 5 %Kou Double Expo

    sigma = params_levy.sigma;
    lam_pois = params_levy.lam;
    p = params_levy.p_up;    omega=p*lam_pois;  nv=(1-p)*lam_pois;
    v = params_levy.eta1;
    w = params_levy.eta2;

    D=1/2*sigma^2;

    %risk-neutral condition
    Psi=@(x)-D*x.^2-nv.*1i.*x./(v+1i.*x)+omega.*1i.*x./(w-1i.*x);
    mu=r-Psi(-1i);

    phi=@(s)mu.*1i.*s-D.*s.^2-nv.*1i.*s./(v+1i.*s)+omega.*1i.*s./(w-1i.*s);

    if T > 0
        diffFgc=@(c)sum(A.*lambda*(mu+sigma^2*c-nv*v./(v+c).^2+omega*w./(w-c).^2)./(r+lambda-phi(-1i*c)).^2.*(1-exp(-(r+lambda-phi(-1i*c))*T)-T*exp(-(r+lambda-phi(-1i*c))*T).*(r+lambda-phi(-1i*c))));
    else
        diffFgc=@(c)sum(A.*lambda*(mu+sigma^2*c-nv*v./(v+c).^2+omega*w./(w-c).^2)./(r+lambda-phi(-1i*c)).^2);
    end

elseif model == 8 % Variance Gamma
    sigma = params_levy.sigmaGBM; % geometric brownian motion add-on
    VG_mu = params_levy.theta;
    VG_sigma = params_levy.sigma;
    nv = params_levy.nu;

    Psi=@(x)-1/2*sigma^2*x.^2-1/nv*log(1-1i*nv*VG_mu*x+nv*VG_sigma^2/2*x.^2);
    mu=r-Psi(-1i);

    phi=@(s)mu*1i*s-1/2*sigma^2*s.^2-1/nv*log(1-1i*nv*VG_mu*s+nv*VG_sigma^2/2*s.^2);

    if T > 0
        diffFgc=@(c)sum(A.*lambda*(mu+sigma^2*c-(VG_mu*1i*c-VG_sigma^2*c)./(1-nv*VG_mu*c-1/2*nv*VG_sigma^2*c.^2))./(r+lambda-phi(-1i*c)).^2.*(1-exp(-(r+lambda-phi(-1i*c))*T)-T*exp(-(r+lambda-phi(-1i*c))*T).*(r+lambda-phi(-1i*c))));
    else
        diffFgc=@(c)sum(A.*lambda*(mu+sigma^2*c-(VG_mu*1i*c-VG_sigma^2*c)./(1-nv*VG_mu*c-1/2*nv*VG_sigma^2*c.^2))./(r+lambda-phi(-1i*c)).^2);
    end
end

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

delta1=@(j,k)(j==k);

S=(0:(N-1))*Delta_s;
v=@(j)1-(delta1(j,1)+delta1(j,N))/2;
Fphi=@(s)12.*sin(s./2).^2./(s.^2.*(2+cos(s)));
Fphi1=Fphi(S/a); Fphi1(1)=1;

%coefs
if T > 0
    Ak=sum(repmat(A.*lambda,length(S),1)./(r+repmat(lambda,length(S),1)-repmat(phi(S).',1,2)).*(1-exp(-(r+repmat(lambda,length(S),1)-repmat(phi(S).',1,2))*T)),2);
    Ck=sum(repmat(A.*lambda,length(S),1)./(r+repmat(lambda,length(S),1)-repmat(phi(S-1i).',1,2)).*(1-exp(-(r+repmat(lambda,length(S),1)-repmat(phi(S-1i).',1,2))*T)),2);

    Fgc=@(c)sum(A.*lambda./(r+lambda-phi(-1i*c)).*(1-exp(-(r+lambda-phi(-1i*c))*T)));
else
    Ak=sum(repmat(A.*lambda,length(S),1)./(r+repmat(lambda,length(S),1)-repmat(phi(S).',1,2)),2);
    Ck=sum(repmat(A.*lambda,length(S),1)./(r+repmat(lambda,length(S),1)-repmat(phi(S-1i).',1,2)),2);

    Fgc=@(c)sum((A.*lambda)./(r+lambda-phi(-1i*c)));
end

% Compute Projections
x10=real(double(diffFgc(0))/Fgc(0))-Delta*(N/2-1);
X0=x10+(0:(N-1))*Delta;
beta0=a^(-1/2)/pi*real(fft((Ak.').*Fphi1.*v(1:N).*Delta_s.*exp(-1i*x10.*S)));

x11=real(double(diffFgc(1))/Fgc(1))-Delta*(N/2-1);
X1=x11+(0:(N-1))*Delta;
beta1=a^(-1/2)/pi*real(fft((Ck.').*Fphi1.*v(1:N).*Delta_s.*exp(-1i*x11.*S)));

% Compute Final Payoff and Value
nn=find((kk<X0)==1);
n=nn(1);
PHi0=zeros(1,length(S));

nn1=find((kk<X1)==1);
n1=nn1(1);
PHi1=zeros(1,length(S));

if call == 1
    PHi0(n)=a^(1/2)*(X0(n)-kk)-a^(3/2)/2*(kk-X0(n))^2+1/2*a^(-1/2);
    PHi0(n-1)=a^(1/2)*(X0(n-1)+1/a-kk)-a^(3/2)/2*(1/(a^2)-(kk-X0(n-1))^2);
    PHi0(n+1:length(S))=1/sqrt(a);

    PHi1(n1)=a^(1/2)*(X1(n1)-kk)-a^(3/2)/2*(kk-X1(n1))^2+1/2*a^(-1/2);
    PHi1(n1-1)=a^(1/2)*(X1(n1-1)+1/a-kk)-a^(3/2)/2*(1/(a^2)-(kk-X1(n1-1))^2);
    PHi1(n1+1:length(S))=1/sqrt(a);

    price=S_0*beta1*PHi1.'-W*beta0*PHi0.';
else
    PHi0(n)=a^(1/2)*(kk-X0(n)+1/a)-a^(3/2)/2*(1/(a^2)-(kk-X0(n))^2);
    PHi0(n-1)=a^(1/2)*(kk-X0(n-1))-a^(3/2)/2*(kk-X0(n-1))^2+1/2*a^(-1/2);
    PHi0(1:n-2)=1/sqrt(a);

    PHi1(n1)=a^(1/2)*(kk-X1(n1)+1/a)-a^(3/2)/2*(1/(a^2)-(kk-X1(n1))^2);
    PHi1(n1-1)=a^(1/2)*(kk-X1(n1-1))-a^(3/2)/2*(kk-X1(n1-1))^2+1/2*a^(-1/2);
    PHi1(1:n1-2)=1/sqrt(a);

    price=W*beta0*PHi0.'-S_0*beta1*PHi1.';
end

end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%  This script is used to test whether the functions in this  %%%%%%%%%%
%%%%%%%%%%  package can run smoothly.                                  %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%  H.D. Li, lhdcsu@gmail.com                                  %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%+++ Import data;
load DM2;
%+++ Cross validation
A=6;
K=5;
method='autoscaling';
N=500;
Nmcs=50;
CV=plsldacv(Xcal,ycal,A,K,method)
MCCV=plsldamccv(Xcal,ycal,A,method,N)
DCV=plsldardcv(Xcal,ycal,A,K,method,Nmcs)

%+++ Build a PLS-LDA model
nLV=3;
LDA=plslda(Xcal,ycal,nLV);
[ScoresTest]=plsldaproj(LDA,Xcal(1:3,:))

%+++ Scores plot
plotlda(LDA,1,0,[2 3 1]);
figure;
plotlda(LDA,1,1,[1  2 3]);
figure;
plotlda(LDA,0,0,[1 2]);
figure;
plotlda(LDA,1,1,[2 3 ]);
figure;
plotlda(LDA,2,1,[3 1]);

%+++ CARS-PLSLDA for variable selection
CARS=carsplslda(Xcal,ycal,A,K,method,50);
figure;
plotcars(CARS);
%+++ simplified version of CARS-PLSLDA for variable selection
sCARS=scarsplslda(Xcal,ycal,A,K,method,50);
figure;
plotcars(sCARS);

%+++ SPA for vairable selection: based on Model Population Analysis
SPA=spa(Xcal,ycal,A,K,method,N,0.7,15);
figure;
bar(SPA.COSS,'b','edgecolor','w');
figure;
plotspa(SPA,SPA.RankedVariable(1));
p=SPA.p(SPA.RankedVariable(1))
%+++ Test ended

author: viggin

repo: domain-adaptation-toolbox


t = 0:0.1:10;
x = exp(-0.2*t) .* cos(2*t);
y = exp(-0.2*t) .* sin(2*t);
plot3(x,y,t,'LineWidth',2);
title('\bfThree-Dimensional Line Plot');
xlabel('\bfx');
ylabel('\bfy');
zlabel('\bftime');
grid on;


text:
function write_neuralynx_ncs(filename, ncs);

% WRITE_NEURALYNX_NCS writes continuous data to a NCS file
% The input data should be scaled in uV.
%
% Use as
%   write_neuralynx_ncs(filename, ncs)
%
% See also READ_NEURALYNX_NCS

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

if ~isa(ncs.TimeStamp, 'uint64')
  error('timestamps should be uint64');
end

% convert the data from uV into V
ncs.dat = ncs.dat * 1e-6;
% scale the data and convert to 16 bits,
% this has to be done prior to writing the header
ADMaxValue = double(intmax('int16'));
ADMaxVolts = max(abs(ncs.dat(:)));
if ADMaxVolts>0
  ADBitVolts = ADMaxVolts / ADMaxValue;
else
  ADBitVolts = 1;
end
ncs.dat = int16(ncs.dat / ADBitVolts);
% update the header with the calibration values
ncs.hdr.ADBitVolts = ADBitVolts;
ncs.hdr.ADMaxValue = ADMaxValue;

% construct the header
buf  = [];
buf  = [buf sprintf('######## Neuralynx Data File Header\r\n')];
buf  = [buf sprintf('## File Name: %s\r\n', filename)];
buf  = [buf sprintf('## Time Opened: (m/d/y): %s  At Time: %s\r\n', datestr(clock, 'mm/dd/yy'), datestr(clock, 'HH:MM:SS'))];
f = fieldnames(ncs.hdr);
for i=1:length(f)
  v = getfield(ncs.hdr, f{i});
  switch class(v)
    case 'char'
      buf = [buf sprintf('-%s\t%s\r\n', f{i}, v)];
    case 'double'
      buf = [buf sprintf('-%s\t%s\r\n', f{i}, num2str(v))];
    otherwise
      error('unknown class in writing header');
  end
end

% pad the rest of the header with zeros
buf((end+1):16384) = 0;

% open the file and write the header
fid  = fopen(filename, 'wb', 'ieee-le');
fwrite(fid, buf);

% The format of a continuous sampled record is
%   int64 TimeStamp
%   int32 ChanNumber
%   int32 SampFreq
%   int32 NumValidSamp
%   int16 Samp[0] ... int16 Samp[511]
% Note that if NumValidSamp < 512, Samp[n], where n >= NumValidSamp, will
% contain random data.

for i=1:size(ncs.dat,2)
  % write a single continuous data record
  fwrite(fid, ncs.TimeStamp(i)   , 'uint64');
  fwrite(fid, ncs.ChanNumber(i)  , 'int32');
  fwrite(fid, ncs.SampFreq(i)    , 'int32');
  fwrite(fid, ncs.NumValidSamp(i), 'int32');
  fwrite(fid, ncs.dat(:,i)       , 'int16');
end

% close the file
fclose(fid);

author: PatternRecognition

repo: OpenBMI


text:
function F = in_fread_msr(sFile, sfid, SamplesBounds)
% IN_FREAD_MSR:  Read a block of recordings from a ANT ASA .msm file
%
% USAGE:  F = in_fread_brainamp(sFile, sfid, SamplesBounds=[])

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

% FORMAT: linear matrix [nChannels x nTime]
nChannels = sFile.header.numchannels;
nTime     = round((sFile.prop.times(2) - sFile.prop.times(1)) .* sFile.prop.sfreq) + 1;
bytesize  = 4;
% Get start position
offsetTime  = SamplesBounds(1) * nChannels * bytesize;
% Number of time values to read for each channel
nTimeToRead = SamplesBounds(2) - SamplesBounds(1) + 1;
% Number of values to skip after each channel
nSkipTimeEnd = (nTime - SamplesBounds(2) - 1) * nChannels * bytesize;
nSkip        = nSkipTimeEnd + offsetTime;
% Position file at the beginning of the data block
fseek(sfid, double(offsetTime), 'bof');
% Read everything at once
% => WARNING: CALL TO FREAD WITH SKIP=0 DOES NOT WORK PROPERLY
if (nSkip == 0)
    F = fread(sfid, [nChannels, nTime], '*float32');
else
    precision = sprintf('%d*float32=>float32', nTimeToRead * nChannels);
    F = fread(sfid, [nChannels, nTimeToRead], precision, nSkip);
end

% Convert from microVolts to Volts
F = 1e-6 * F;

author: brainstorm-tools


text:
function X = FDDL_updateX(Y, Y_range, D, D_range, X, opts)
    % X = argmin_X 0.5\|Yhat - Dhat X\| +
    %    + 0.5*lambda2(\sum (normF2(Xi - Mi) - normF2(Mi - M)) + normF2(X) + normF2(X0 - M0))
    %    + lambda1*norm1(X)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% ================== block: test module ==========================
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if nargin == 0
        addpath('../utils');
        addpath('../sparse_coding');
        tic
        d       = 30;
        N       = 7;
        k       = 5;
        C       = 3 ;
        Y       = normc(rand(d,N*C));
        D       = normc(rand(d,k*C));
        Y_range = N* (0:C);
        D_range = k* (0:C);
        X       = randn(size(D,2), size(Y,2));

        opts.k0       = k0;
        opts.lambda1  = 0.01;
        opts.lambda2  = 0.002;
        opts.lambda3  = 0.1;
        opts.max_iter = 250;
        opts.show     = true;
        opts.check_grad = true;
        opts          = initOpts(opts); % other attributes
    end

    lambda1  = opts.lambda1;
    lambda2  = opts.lambda2;
    DtD      = D'*D;
    D_0      = buildMhat(DtD, D_range, D_range);
    Dhat     = D_0 + 2*opts.lambda2*eye(size(D_0,1));
    DtY     = D'*Y;
    Y_0     = buildMhat(DtY, D_range, Y_range);
    %% cost w.r.t. X, X0, not includeing norm1 terms
    function cost = calc_f(X)
        cost = 0.5*(normF2(Y - D*X) + ...
                    FDDL_fidelity(Y, Y_range, D, D_range, X)) + ...
               0.5*opts.lambda2* (FDDL_discriminative(X, Y_range));
    end
    %% Total cose
    function cost = calc_F(X)
        cost = calc_f(X) + lambda1*norm1(X);
    end
    %%
    %% Gradient for FISTA
    function g = grad(X)
        g       = Dhat*X - Y_0 + buildM_2Mbar(X, Y_range, lambda2);
    end
    %% check gradient
    if opts.check_grad &&~check_grad(@calc_f, @grad, X)
        fprintf('Check gradient or cost again!\n')
        pause
    end
    %% ========= Main FISTA ==============================
    optsXX0          = opts;
    optsXX0.max_iter = 300;
    L = max(eig(Dhat)) + 6*lambda2;
    X      = fista(@grad, X, L, opts.lambda1, opts, @calc_F);
    %%
    if nargin == 0
        fprintf('done, press any key to see results\n');
        pause;
    end
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function [step, lowerBound, upperBound] = defaultBoundaries(nPoints)
%Generates default lower bound, upper bound, and step for qth derivatives.

    step = 1;
    lowerBound = 1;
    upperBound = nPoints;

end


text:
% FORMAT:
%
% sampval = ms2sample(timems, fs, rounding, offset)
%
% INPUTS:
%
% timesam     - time value in samples
% fs          - sample rate
%
% Optional inputs:
%
% rounding    - 1 means round the result {default}, 0 means do not round
% offset      - offset time in milliseconds (you can use sample2ms recursively here). Default 0
%
% OUTPUT:
%
% sampval     - time in milliseconds
%
% EXAMPLES:
%
% 1) For a time serie recorded from 0 to 5 secs, at fs=500 sps, get the time in ms for the sample # 337.
%
% msval = sample2ms(337, 500)
%
% msval =
%
%    674
%
% 2) For a time serie recorded from -1 to 5 secs, at fs=500 sps, get the absolute time in ms for the sample # 337.
%
% msval = sample2ms(337, 500, 0, sample2ms(500,500))
%
% msval =
%
%    1674
%
%
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2013

function msval = sample2ms(timesam, fs, rounding, offset)
if nargin<1
        help sample2ms
        return
end
if nargin<4
        offset = 0;
end
if nargin<3
        rounding = 1;
end
if nargin<2
        error('Two inputs are requiered at least.')
end
msval = offset + 1000*timesam/fs;
if rounding
        msval = round(msval);
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function DoDisableButton(varargin)

handles=varargin{3};
handles=guidata(handles.SimuPanel_figure);

set(handles.Scan_pushbutton,'Enable','off');
set(handles.Batch_pushbutton,'Enable','off');

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function count = cprintf(style,format,varargin)
% CPRINTF displays styled formatted text in the Command Window
%
% Syntax:
%    count = cprintf(style,format,...)
%
% Description:
%    CPRINTF processes the specified text using the exact same FORMAT
%    arguments accepted by the built-in SPRINTF and FPRINTF functions.
%
%    CPRINTF then displays the text in the Command Window using the
%    specified STYLE argument. The accepted styles are those used for
%    Matlab's syntax highlighting (see: File / Preferences / Colors /
%    M-file Syntax Highlighting Colors), and also user-defined colors.
%
%    The possible pre-defined STYLE names are:
%
%       'Text'                 - default: black
%       'Keywords'             - default: blue
%       'Comments'             - default: green
%       'Strings'              - default: purple
%       'UnterminatedStrings'  - default: dark red
%       'SystemCommands'       - default: orange
%       'Errors'               - default: light red
%       'Hyperlinks'           - default: underlined blue
%
%       'Black','Cyan','Magenta','Blue','Green','Red','Yellow','White'
%
%    STYLE beginning with '-' or '_' will be underlined. For example:
%          '-Blue' is underlined blue, like 'Hyperlinks';
%          '_Comments' is underlined green etc.
%
%    STYLE beginning with '*' will be bold (R2011b+ only). For example:
%          '*Blue' is bold blue;
%          '*Comments' is bold green etc.
%    Note: Matlab does not currently support both bold and underline,
%          only one of them can be used in a single cprintf command. But of
%          course bold and underline can be mixed by using separate commands.
%
%    STYLE also accepts a regular Matlab RGB vector, that can be underlined
%    and bolded: -[0,1,1] means underlined cyan, '*[1,0,0]' is bold red.
%
%    STYLE is case-insensitive and accepts unique partial strings just
%    like handle property names.
%
%    CPRINTF by itself, without any input parameters, displays a demo
%
% Example:
%    cprintf;   % displays the demo
%    cprintf('text',   'regular black text');
%    cprintf('hyper',  'followed %s','by');
%    cprintf('key',    '%d colored', 4);
%    cprintf('-comment','& underlined');
%    cprintf('err',    'elements\n');
%    cprintf('cyan',   'cyan');
%    cprintf('_green', 'underlined green');
%    cprintf(-[1,0,1], 'underlined magenta');
%    cprintf([1,0.5,0],'and multi-\nline orange\n');
%    cprintf('*blue',  'and *bold* (R2011b+ only)\n');
%    cprintf('string');  % same as fprintf('string') and cprintf('text','string')
%
% Bugs and suggestions:
%    Please send to Yair Altman (altmany at gmail dot com)
%
% Warning:
%    This code heavily relies on undocumented and unsupported Matlab
%    functionality. It works on Matlab 7+, but use at your own risk!
%
%    A technical description of the implementation can be found at:
%    http://UndocumentedMatlab.com/blog/cprintf/
%
% Limitations:
%    1. In R2011a and earlier, a single space char is inserted at the
%       beginning of each CPRINTF text segment (this is ok in R2011b+).
%
%    2. In R2011a and earlier, consecutive differently-colored multi-line
%       CPRINTFs sometimes display incorrectly on the bottom line.
%       As far as I could tell this is due to a Matlab bug. Examples:
%         >> cprintf('-str','under\nline'); cprintf('err','red\n'); % hidden 'red', unhidden '_'
%         >> cprintf('str','regu\nlar'); cprintf('err','red\n'); % underline red (not purple) 'lar'
%
%    3. Sometimes, non newline ('\n')-terminated segments display unstyled
%       (black) when the command prompt chevron ('>>') regains focus on the
%       continuation of that line (I can't pinpoint when this happens).
%       To fix this, simply newline-terminate all command-prompt messages.
%
%    4. In R2011b and later, the above errors appear to be fixed. However,
%       the last character of an underlined segment is not underlined for
%       some unknown reason (add an extra space character to make it look better)
%
%    5. In old Matlab versions (e.g., Matlab 7.1 R14), multi-line styles
%       only affect the first line. Single-line styles work as expected.
%       R14 also appends a single space after underlined segments.
%
%    6. Bold style is only supported on R2011b+, and cannot also be underlined.
%
% Change log:
%    2015-06-24: Fixed a few discoloration issues (some other issues still remain)
%    2015-03-20: Fix: if command window isn't defined yet (startup) use standard fprintf as suggested by John Marozas
%    2012-08-09: Graceful degradation support for deployed (compiled) and non-desktop applications; minor bug fixes
%    2012-08-06: Fixes for R2012b; added bold style; accept RGB string (non-numeric) style
%    2011-11-27: Fixes for R2011b
%    2011-08-29: Fix by Danilo (FEX comment) for non-default text colors
%    2011-03-04: Performance improvement
%    2010-06-27: Fix for R2010a/b; fixed edge case reported by Sharron; CPRINTF with no args runs the demo
%    2009-09-28: Fixed edge-case problem reported by Swagat K
%    2009-05-28: corrected nargout behavior suggested by Andreas G?b
%    2009-05-13: First version posted on MathWorks File Exchange
%
% See also:
%    sprintf, fprintf

% License to use and modify this code is granted freely to all interested, as long as the original author is
% referenced and attributed as such. The original author maintains the right to be solely associated with this work.

% Programmed and Copyright by Yair M. Altman: altmany(at)gmail.com
% $Revision: 1.10 $  $Date: 2015/06/24 01:29:18 $

  persistent majorVersion minorVersion
  if isempty(majorVersion)
      %v = version; if str2double(v(1:3)) <= 7.1
      %majorVersion = str2double(regexprep(version,'^(\d+).*','$1'));
      %minorVersion = str2double(regexprep(version,'^\d+\.(\d+).*','$1'));
      %[a,b,c,d,versionIdStrs]=regexp(version,'^(\d+)\.(\d+).*');  %#ok unused
      v = sscanf(version, '%d.', 2);
      majorVersion = v(1); %str2double(versionIdStrs{1}{1});
      minorVersion = v(2); %str2double(versionIdStrs{1}{2});
  end

  % The following is for debug use only:
  %global docElement txt el
  if ~exist('el','var') || isempty(el),  el=handle([]);  end  %#ok mlint short-circuit error ("used before defined")
  if nargin<1, showDemo(majorVersion,minorVersion); return;  end
  if isempty(style),  return;  end
  if all(ishandle(style)) && length(style)~=3
      dumpElement(style);
      return;
  end

  % Process the text string
  if nargin<2, format = style; style='text';  end
  %error(nargchk(2, inf, nargin, 'struct'));
  %str = sprintf(format,varargin{:});

  % In compiled mode
  try useDesktop = usejava('desktop'); catch, useDesktop = false; end
  if isdeployed | ~useDesktop %#ok<OR2> - for Matlab 6 compatibility
      % do not display any formatting - use simple fprintf()
      % See: http://undocumentedmatlab.com/blog/bold-color-text-in-the-command-window/#comment-103035
      % Also see: https://mail.google.com/mail/u/0/?ui=2&shva=1#all/1390a26e7ef4aa4d
      % Also see: https://mail.google.com/mail/u/0/?ui=2&shva=1#all/13a6ed3223333b21
      count1 = fprintf(format,varargin{:});
  else
      % Else (Matlab desktop mode)
      % Get the normalized style name and underlining flag
      [underlineFlag, boldFlag, style, debugFlag] = processStyleInfo(style);

      % Set hyperlinking, if so requested
      if underlineFlag
          format = ['' format ''];

          % Matlab 7.1 R14 (possibly a few newer versions as well?)
          % have a bug in rendering consecutive hyperlinks
          % This is fixed by appending a single non-linked space
          if majorVersion < 7 || (majorVersion==7 && minorVersion <= 1)
              format(end+1) = ' ';
          end
      end

      % Set bold, if requested and supported (R2011b+)
      if boldFlag
          if (majorVersion > 7 || minorVersion >= 13)
              format = ['' format ''];
          else
              boldFlag = 0;
          end
      end

      % Get the current CW position
      cmdWinDoc = com.mathworks.mde.cmdwin.CmdWinDocument.getInstance;
      lastPos = cmdWinDoc.getLength;

      % If not beginning of line
      bolFlag = 0;  %#ok
      %if docElement.getEndOffset - docElement.getStartOffset > 1
          % Display a hyperlink element in order to force element separation
          % (otherwise adjacent elements on the same line will be merged)
          if majorVersion<7 || (majorVersion==7 && minorVersion<13)
              if ~underlineFlag
                  fprintf(' ');  %fprintf(' \b');
              elseif format(end)~=10  % if no newline at end
                  fprintf(' ');  %fprintf(' \b');
              end
          end
          %drawnow;
          bolFlag = 1;
      %end

      % Get a handle to the Command Window component
      mde = com.mathworks.mde.desk.MLDesktop.getInstance;
      cw = mde.getClient('Command Window');

      % Fix: if command window isn't defined yet (startup), use standard fprintf()
      if (isempty(cw))
         count1 = fprintf(format,varargin{:});
         if nargout
             count = count1;
         end
         return;
      end

      xCmdWndView = cw.getComponent(0).getViewport.getComponent(0);

      % Store the CW background color as a special color pref
      % This way, if the CW bg color changes (via File/Preferences),
      % it will also affect existing rendered strs
      com.mathworks.services.Prefs.setColorPref('CW_BG_Color',xCmdWndView.getBackground);

      % Display the text in the Command Window
      % Note: fprintf(2,...) is required in order to add formatting tokens, which
      % ^^^^  can then be updated below (no such tokens when outputting to stdout)
      count1 = fprintf(2,format,varargin{:});

      % Repaint the command window
      % awtinvoke(cmdWinDoc,'remove',lastPos,1);   % TODO: find out how to remove the extra '_'
      drawnow  % this is necessary for the following to work properly (refer to Evgeny Pr in FEX comment 16/1/2011)
      xCmdWndView.repaint;
      %hListeners = cmdWinDoc.getDocumentListeners; for idx=1:numel(hListeners), try hListeners(idx).repaint; catch, end, end

      docElement = cmdWinDoc.getParagraphElement(lastPos+1);
      if majorVersion<7 || (majorVersion==7 && minorVersion<13)
          if bolFlag && ~underlineFlag
              % Set the leading hyperlink space character ('_') to the bg color, effectively hiding it
              % Note: old Matlab versions have a bug in hyperlinks that need to be accounted for...
              %disp(' '); dumpElement(docElement)
              setElementStyle(docElement,'CW_BG_Color',1+underlineFlag,majorVersion,minorVersion); %+getUrlsFix(docElement));
              %disp(' '); dumpElement(docElement)
              el(end+1) = handle(docElement);  % #ok used in debug only
          end

          % Fix a problem with some hidden hyperlinks becoming unhidden...
          fixHyperlink(docElement);
          %dumpElement(docElement);
      end

      % Get the Document Element(s) corresponding to the latest fprintf operation
      while docElement.getStartOffset < cmdWinDoc.getLength
          % Set the element style according to the current style
          if debugFlag, dumpElement(docElement); end
          specialFlag = underlineFlag | boldFlag;
          setElementStyle(docElement,style,specialFlag,majorVersion,minorVersion);
          if debugFlag, dumpElement(docElement); end
          docElement2 = cmdWinDoc.getParagraphElement(docElement.getEndOffset+1);
          if isequal(docElement,docElement2),  break;  end
          docElement = docElement2;
      end
      if debugFlag, dumpElement(docElement); end

      % Force a Command-Window repaint
      % Note: this is important in case the rendered str was not '\n'-terminated
      xCmdWndView.repaint;

      % The following is for debug use only:
      el(end+1) = handle(docElement);  %#ok used in debug only
      %elementStart  = docElement.getStartOffset;
      %elementLength = docElement.getEndOffset - elementStart;
      %txt = cmdWinDoc.getText(elementStart,elementLength);
  end

  if nargout
      count = count1;
  end
  return;  % debug breakpoint

% Process the requested style information
function [underlineFlag,boldFlag,style,debugFlag] = processStyleInfo(style)
  underlineFlag = 0;
  boldFlag = 0;
  debugFlag = 0;

  % First, strip out the underline/bold markers
  if ischar(style)
      % Styles containing '-' or '_' should be underlined (using a no-target hyperlink hack)
      %if style(1)=='-'
      underlineIdx = (style=='-') | (style=='_');
      if any(underlineIdx)
          underlineFlag = 1;
          %style = style(2:end);
          style = style(~underlineIdx);
      end

      % Check for bold style (only if not underlined)
      boldIdx = (style=='*');
      if any(boldIdx)
          boldFlag = 1;
          style = style(~boldIdx);
      end
      if underlineFlag && boldFlag
          warning('YMA:cprintf:BoldUnderline','Matlab does not support both bold & underline')
      end

      % Check for debug mode (style contains '!')
      debugIdx = (style=='!');
      if any(debugIdx)
          debugFlag = 1;
          style = style(~debugIdx);
      end

      % Check if the remaining style sting is a numeric vector
      %styleNum = str2num(style); %#ok<ST2NM>  % not good because style='text' is evaled!
      %if ~isempty(styleNum)
      if any(style==' ' | style==',' | style==';')
          style = str2num(style); %#ok<ST2NM>
      end
  end

  % Style = valid matlab RGB vector
  if isnumeric(style) && length(style)==3 && all(style<=1) && all(abs(style)>=0)
      if any(style<0)
          underlineFlag = 1;
          style = abs(style);
      end
      style = getColorStyle(style);

  elseif ~ischar(style)
      error('YMA:cprintf:InvalidStyle','Invalid style - see help section for a list of valid style values')

  % Style name
  else
      % Try case-insensitive partial/full match with the accepted style names
      matlabStyles = {'Text','Keywords','Comments','Strings','UnterminatedStrings','SystemCommands','Errors'};
      validStyles  = [matlabStyles, ...
                      'Black','Cyan','Magenta','Blue','Green','Red','Yellow','White', ...
                      'Hyperlinks'];
      matches = find(strncmpi(style,validStyles,length(style)));

      % No match - error
      if isempty(matches)
          error('YMA:cprintf:InvalidStyle','Invalid style - see help section for a list of valid style values')

      % Too many matches (ambiguous) - error
      elseif length(matches) > 1
          error('YMA:cprintf:AmbigStyle','Ambiguous style name - supply extra characters for uniqueness')

      % Regular text
      elseif matches == 1
          style = 'ColorsText';  % fixed by Danilo, 29/8/2011

      % Highlight preference style name
      elseif matches <= length(matlabStyles)
          style = ['Colors_M_' validStyles{matches}];

      % Color name
      elseif matches < length(validStyles)
          colors = [0,0,0; 0,1,1; 1,0,1; 0,0,1; 0,1,0; 1,0,0; 1,1,0; 1,1,1];
          requestedColor = colors(matches-length(matlabStyles),:);
          style = getColorStyle(requestedColor);

      % Hyperlink
      else
          style = 'Colors_HTML_HTMLLinks';  % CWLink
          underlineFlag = 1;
      end
  end

% Convert a Matlab RGB vector into a known style name (e.g., '[255,37,0]')
function styleName = getColorStyle(rgb)
  intColor = int32(rgb*255);
  javaColor = java.awt.Color(intColor(1), intColor(2), intColor(3));
  styleName = sprintf('[%d,%d,%d]',intColor);
  com.mathworks.services.Prefs.setColorPref(styleName,javaColor);

% Fix a bug in some Matlab versions, where the number of URL segments
% is larger than the number of style segments in a doc element
function delta = getUrlsFix(docElement)  %#ok currently unused
  tokens = docElement.getAttribute('SyntaxTokens');
  links  = docElement.getAttribute('LinkStartTokens');
  if length(links) > length(tokens(1))
      delta = length(links) > length(tokens(1));
  else
      delta = 0;
  end

% fprintf(2,str) causes all previous '_'s in the line to become red - fix this
function fixHyperlink(docElement)
  try
      tokens = docElement.getAttribute('SyntaxTokens');
      urls   = docElement.getAttribute('HtmlLink');
      urls   = urls(2);
      links  = docElement.getAttribute('LinkStartTokens');
      offsets = tokens(1);
      styles  = tokens(2);
      doc = docElement.getDocument;

      % Loop over all segments in this docElement
      for idx = 1 : length(offsets)-1
          % If this is a hyperlink with no URL target and starts with ' ' and is collored as an error (red)...
          if strcmp(styles(idx).char,'Colors_M_Errors')
              character = char(doc.getText(offsets(idx)+docElement.getStartOffset,1));
              if strcmp(character,' ')
                  if isempty(urls(idx)) && links(idx)==0
                      % Revert the style color to the CW background color (i.e., hide it!)
                      styles(idx) = java.lang.String('CW_BG_Color');
                  end
              end
          end
      end
  catch
      % never mind...
  end

% Set an element to a particular style (color)
function setElementStyle(docElement,style,specialFlag, majorVersion,minorVersion)
  %global tokens links urls urlTargets  % for debug only
  global oldStyles
  if nargin<3,  specialFlag=0;  end
  % Set the last Element token to the requested style:
  % Colors:
  tokens = docElement.getAttribute('SyntaxTokens');
  try
      styles = tokens(2);
      oldStyles{end+1} = cell(styles);

      % Correct edge case problem
      extraInd = double(majorVersion>7 || (majorVersion==7 && minorVersion>=13));  % =0 for R2011a-, =1 for R2011b+
      %{
      if ~strcmp('CWLink',char(styles(end-hyperlinkFlag))) && ...
          strcmp('CWLink',char(styles(end-hyperlinkFlag-1)))
         extraInd = 0;%1;
      end
      hyperlinkFlag = ~isempty(strmatch('CWLink',tokens(2)));
      hyperlinkFlag = 0 + any(cellfun(@(c)(~isempty(c)&&strcmp(c,'CWLink')),cell(tokens(2))));
      %}

      jStyle = java.lang.String(style);
      if numel(styles)==4 && isempty(char(styles(2)))
          % Attempt to fix discoloration issues - NOT SURE THAT THIS IS OK! - 24/6/2015
          styles(1) = jStyle;
      end
      styles(end-extraInd) = java.lang.String('');
      styles(end-extraInd-specialFlag) = jStyle;  % #ok apparently unused but in reality used by Java
      if extraInd
          styles(end-specialFlag) = jStyle;
      end

      oldStyles{end} = [oldStyles{end} cell(styles)];
  catch
      % never mind for now
  end

  % Underlines (hyperlinks):
  %{
  links = docElement.getAttribute('LinkStartTokens');
  if isempty(links)
      %docElement.addAttribute('LinkStartTokens',repmat(int32(-1),length(tokens(2)),1));
  else
      %TODO: remove hyperlink by setting the value to -1
  end
  %}

  % Correct empty URLs to be un-hyperlinkable (only underlined)
  %{
  urls = docElement.getAttribute('HtmlLink');
  if ~isempty(urls)
      urlTargets = urls(2);
      for urlIdx = 1 : length(urlTargets)
          try
              if urlTargets(urlIdx).length < 1
                  urlTargets(urlIdx) = [];  % '' => []
              end
          catch
              % never mind...
              a=1;  %#ok used for debug breakpoint...
          end
      end
  end
  %}
  % Bold: (currently unused because we cannot modify this immutable int32 numeric array)
  %{
  try
      %hasBold = docElement.isDefined('BoldStartTokens');
      bolds = docElement.getAttribute('BoldStartTokens');
      if ~isempty(bolds)
          %docElement.addAttribute('BoldStartTokens',repmat(int32(1),length(bolds),1));
      end
  catch
      % never mind - ignore...
      a=1;  %#ok used for debug breakpoint...
  end
  %}

  return;  % debug breakpoint

% Display information about element(s)
function dumpElement(docElements)
  %return;
  disp(' ');
  numElements = length(docElements);
  cmdWinDoc = docElements(1).getDocument;
  for elementIdx = 1 : numElements
      if numElements > 1,  fprintf('Element #%d:\n',elementIdx);  end
      docElement = docElements(elementIdx);
      if ~isjava(docElement),  docElement = docElement.java;  end
      %docElement.dump(java.lang.System.out,1)
      disp(docElement)
      tokens = docElement.getAttribute('SyntaxTokens');
      if isempty(tokens),  continue;  end
      links = docElement.getAttribute('LinkStartTokens');
      urls  = docElement.getAttribute('HtmlLink');
      try bolds = docElement.getAttribute('BoldStartTokens'); catch, bolds = []; end
      txt = {};
      tokenLengths = tokens(1);
      for tokenIdx = 1 : length(tokenLengths)-1
          tokenLength = diff(tokenLengths(tokenIdx+[0,1]));
          if (tokenLength < 0)
              tokenLength = docElement.getEndOffset - docElement.getStartOffset - tokenLengths(tokenIdx);
          end
          txt{tokenIdx} = cmdWinDoc.getText(docElement.getStartOffset+tokenLengths(tokenIdx),tokenLength).char;  %#ok
      end
      lastTokenStartOffset = docElement.getStartOffset + tokenLengths(end);
      try
          txt{end+1} = cmdWinDoc.getText(lastTokenStartOffset, docElement.getEndOffset-lastTokenStartOffset).char; %#ok
      catch
          txt{end+1} = ''; %#ok<AGROW>
      end
      %cmdWinDoc.uiinspect
      %docElement.uiinspect
      txt = strrep(txt',sprintf('\n'),'\n');
      try
          data = [cell(tokens(2)) m2c(tokens(1)) m2c(links) m2c(urls(1)) cell(urls(2)) m2c(bolds) txt];
          if elementIdx==1
              disp('    SyntaxTokens(2,1) - LinkStartTokens - HtmlLink(1,2) - BoldStartTokens - txt');
              disp('    ==============================================================================');
          end
      catch
          try
              data = [cell(tokens(2)) m2c(tokens(1)) m2c(links) txt];
          catch
              disp([cell(tokens(2)) m2c(tokens(1)) txt]);
              try
                  data = [m2c(links) m2c(urls(1)) cell(urls(2))];
              catch
                  % Matlab 7.1 only has urls(1)...
                  data = [m2c(links) cell(urls)];
              end
          end
      end
      disp(data)
  end

% Utility function to convert matrix => cell
function cells = m2c(data)
  %datasize = size(data);  cells = mat2cell(data,ones(1,datasize(1)),ones(1,datasize(2)));
  cells = num2cell(data);

% Display the help and demo
function showDemo(majorVersion,minorVersion)
  fprintf('cprintf displays formatted text in the Command Window.\n\n');
  fprintf('Syntax: count = cprintf(style,format,...);  click here for details.\n\n');
  url = 'http://UndocumentedMatlab.com/blog/cprintf/';
  fprintf(['Technical description: ' url '\n\n']);
  fprintf('Demo:\n\n');
  boldFlag = majorVersion>7 || (majorVersion==7 && minorVersion>=13);
  s = ['cprintf(''text'',    ''regular black text'');' 10 ...
       'cprintf(''hyper'',   ''followed %s'',''by'');' 10 ...
       'cprintf(''key'',     ''%d colored'',' num2str(4+boldFlag) ');' 10 ...
       'cprintf(''-comment'',''& underlined'');' 10 ...
       'cprintf(''err'',     ''elements:\n'');' 10 ...
       'cprintf(''cyan'',    ''cyan'');' 10 ...
       'cprintf(''_green'',  ''underlined green'');' 10 ...
       'cprintf(-[1,0,1],  ''underlined magenta'');' 10 ...
       'cprintf([1,0.5,0], ''and multi-\nline orange\n'');' 10];
   if boldFlag
       % In R2011b+ the internal bug that causes the need for an extra space
       % is apparently fixed, so we must insert the sparator spaces manually...
       % On the other hand, 2011b enables *bold* format
       s = [s 'cprintf(''*blue'',   ''and *bold* (R2011b+ only)\n'');' 10];
       s = strrep(s, ''')',' '')');
       s = strrep(s, ''',5)',' '',5)');
       s = strrep(s, '\n ','\n');
   end
   disp(s);
   eval(s);

%%%%%%%%%%%%%%%%%%%%%%%%%% TODO %%%%%%%%%%%%%%%%%%%%%%%%%
% - Fix: Remove leading space char (hidden underline '_')
% - Fix: Find workaround for multi-line quirks/limitations
% - Fix: Non-\n-terminated segments are displayed as black
% - Fix: Check whether the hyperlink fix for 7.1 is also needed on 7.2 etc.
% - Enh: Add font support

author: goGPS-Project


text:
%% element2HexLattice
% Below is a demonstration of the features of the |element2HexLattice| function

%%
clear; close all; clc;

%% Syntax
% |[Es,Vs,Cs]=element2HexLattice(E,V,cPar);|

%% Description
% This function converts an element description (elements and vertices i.e
% nodes into a lattice structure. The lattice structure is returned as a
% hexahederal mesh.

%% Examples
%

%%
% Plot settings
cMap=gjet(4);
fontSize=15;

%% Example 1 Creating a lattice structure on hexahedral meshes

%%
% Creating example geometry.
[V,~]=platonic_solid(2,1); %Vertices of cube
E=1:8; %Element description of the 8-node cube (hexahedral element)
C=(1:size(E,1))'; %color (e.g. material) labels for all elements
[F,~]=element2patch(E,C); %Patch data for plotting

%%
% Create lattice structure
controlParameter.growSteps=0; %0 is normal, positive or negative integers increase or decrease the edge lattice thickness respectively
controlParameter.latticeSide=[]; %Empty outputs both, 1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
[Es,Vs,Cs]=element2HexLattice(E,V,controlParameter); %Get lattice structure

%%
% Visualizing input mesh and lattic structures

% Create patch Data for visualization
[Fs,CsF]=element2patch(Es,Cs); %Patch data for plotting
[Fs1,CsF1]=element2patch(Es(Cs==1,:),Cs(Cs==1,:)); %Patch data for plotting
[Fs2,CsF2]=element2patch(Es(Cs==2,:),Cs(Cs==2,:)); %Patch data for plotting

cFigure;
hs=subplot(2,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
ha=axis; axis off;

subplot(2,2,2);
title('The two complementary lattice structures','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs,Vs,CsF);
colormap(cMap);
cLim=caxis; caxis([1 2]);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

subplot(2,2,3);
title('Lattice side 1','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs1,Vs,cMap(1,:));
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

subplot(2,2,4);
title('Lattice side 2','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs2,Vs,CsF2);
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

drawnow;

%% Example 2 Creating a lattice structure on tetrahedral meshes

%%
% Creating example geometry.

[V,~]=platonic_solid(1,1); %Vertices of tetrahedron
E=1:4; %Element description of the 4-node tetrahedron
C=(1:size(E,1))'; %color (e.g. material) labels for all elements
[F,~]=element2patch(E,C); %Patch data for plotting

%%
% Create lattice structure
controlParameter.growSteps=0; %0 is normal, positive or negative integers increase or decrease the edge lattice thickness respectively
controlParameter.latticeSide=[]; %Empty outputs both, 1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
[Es,Vs,Cs]=element2HexLattice(E,V,controlParameter); %Get lattice structure

%%
% Visualizing input mesh and lattic structures

% Create patch Data for visualization
[Fs,CsF]=element2patch(Es,Cs); %Patch data for plotting
[Fs1,CsF1]=element2patch(Es(Cs==1,:),Cs(Cs==1,:)); %Patch data for plotting
[Fs2,CsF2]=element2patch(Es(Cs==2,:),Cs(Cs==2,:)); %Patch data for plotting

cFigure;
hs=subplot(2,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
ha=axis; axis off;

subplot(2,2,2);
title('The two complementary lattice structures','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs,Vs,CsF);
colormap(cMap);
cLim=caxis; caxis([1 2]);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

subplot(2,2,3);
title('Lattice side 1','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs1,Vs,cMap(1,:));
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

subplot(2,2,4);
title('Lattice side 2','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs2,Vs,CsF2);
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

drawnow;

%% Example 3 Changing lattice structure thickness

%%
% Creating example geometry.

[V,~]=platonic_solid(1,1); %Vertices of tetrahedron
E=1:4; %Element description of the 4-node tetrahedron
C=(1:size(E,1))'; %color (e.g. material) labels for all elements
[F,~]=element2patch(E,C); %Patch data for plotting

%%
% Create lattice structure
controlParameter.latticeSide=[]; %Empty outputs both, 1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice

%%
% lattice structure for 0 growth steps
controlParameter.growSteps=0; %0 is normal, positive or negative integers increase or decrease the edge lattice thickness respectively
[Es_0,Vs_0,Cs_0]=element2HexLattice(E,V,controlParameter); %Get lattice structure

controlParameter.growSteps=1; %0 is normal, positive or negative integers increase or decrease the edge lattice thickness respectively
[Es_1,Vs_1,Cs_1]=element2HexLattice(E,V,controlParameter); %Get lattice structure

controlParameter.growSteps=-1; %0 is normal, positive or negative integers increase or decrease the edge lattice thickness respectively
[Es_n1,Vs_n1,Cs_n1]=element2HexLattice(E,V,controlParameter); %Get lattice structure

%%
% Visualizing input mesh and lattic structures

% Create patch Data for visualization
[Fs_0,CsF_0]=element2patch(Es_0,Cs_0); %Patch data for plotting
[Fs_1,CsF_1]=element2patch(Es_1,Cs_1); %Patch data for plotting
[Fs_n1,CsF_n1]=element2patch(Es_n1,Cs_n1); %Patch data for plotting

cFigure;
subplot(1,3,1);
title('growSteps=1','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs_1,Vs_1,CsF_1);
colormap(cMap);
cLim=caxis; caxis([1 2]);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis off;

subplot(1,3,2);
title('growSteps=0','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs_0,Vs_0,CsF_0);
colormap(cMap);
cLim=caxis; caxis([1 2]);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis off;

subplot(1,3,3);
title('growSteps=-1','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
gpatch(Fs_n1,Vs_n1,CsF_n1);
colormap(cMap);
cLim=caxis; caxis([1 2]);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis off;

drawnow;

%% Example 4 Lattices on multiple elements, adjusting lattice type

%%
% Creating example geometry.
[V,~]=platonic_solid(2,1); %Vertices of cube
E=1:8; %Element description of the 8-node cube (hexahedral element)
[E,V,C]=subHex(E,V); %Subdevide into 8 sub-cubes
[E,V,C]=hex2tet(E,V,C,1); %Convert to tetrahedral elements

[F,~]=element2patch(E,C); %Patch data for plotting

%%
% Create lattice structure
controlParameter.growSteps=-1; %0 is normal, positive or negative integers increase or decrease the edge lattice thickness respectively
controlParameter.latticeSide=1; %Empty outputs both, 1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
[Es,Vs,Cs]=element2HexLattice(E,V,controlParameter); %Get lattice structure

%%
% Visualizing input mesh and lattic structures

% Create patch Data for visualization
[Fs,CsF]=element2patch(Es,Cs); %Patch data for plotting

cFigure;
hs=subplot(1,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

subplot(1,2,2);
title('Lattice side 1','fontSize',fontSize)
hold on;
gpatch(Fs,Vs,cMap(1,:));
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

drawnow;

%%
% Create lattice structure
controlParameter.growSteps=1; %0 is normal, positive or negative integers increase or decrease the edge lattice thickness respectively
controlParameter.latticeSide=2; %Empty outputs both, 1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
[Es,Vs,Cs]=element2HexLattice(E,V,controlParameter); %Get lattice structure

%%
% Visualizing input mesh and lattic structures

% Create patch Data for visualization
[Fs,CsF]=element2patch(Es,Cs); %Patch data for plotting

cFigure;
hs=subplot(1,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

subplot(1,2,2);
title('Lattice side 2','fontSize',fontSize)
hold on;
gpatch(Fs,Vs,cMap(1,:));
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

drawnow;

%% Example 5 Hierarchical lattice structures

%%
% Creating example geometry.

% [V,~]=platonic_solid(2,1); %Vertices of cube
% E=1:8; %Element description of the 8-node cube (hexahedral element)
% C=(1:size(E,1))'; %color (e.g. material) labels for all elements
% [F,~]=element2patch(E,C); %Patch data for plotting

[V,~]=platonic_solid(1,1); %Vertices of tetrahedron
E=1:4; %Element description of the 4-node tetrahedron
C=(1:size(E,1))'; %color (e.g. material) labels for all elements
[F,~]=element2patch(E,C); %Patch data for plotting

%%
% Create first order lattice structure
controlParameter.growSteps=-1; %0 is normal, positive or negative integers increase or decrease the edge lattice thickness respectively
controlParameter.latticeSide=1; %Empty outputs both, 1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
[Es1,Vs1,Cs1]=element2HexLattice(E,V,controlParameter); %Get lattice structure

%%
% Uncomment to create more complex structure
[Es1,Vs1,Cs1]=hex2tet(Es1,Vs1,Cs1,5); %Convert to tetrahedral elements
% [Es1,Vs1,Cs1]=tet2hex(Es1,Vs1,1); %Convert to tetrahedral elements

%%
% Create second order lattice structure
controlParameter.growSteps=-1; %0 is normal, positive or negative integers increase or decrease the edge lattice thickness respectively
controlParameter.latticeSide=1; %Empty outputs both, 1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
[Es2,Vs2,Cs2]=element2HexLattice(Es1,Vs1,controlParameter); %Get lattice structure

%%
% Visualizing input mesh and lattic structures

% Create patch Data for visualization
[Fs1,CsF1]=element2patch(Es1,Cs1); %Patch data for plotting
[Fs2,CsF2]=element2patch(Es2,Cs2); %Patch data for plotting

[indBounary]=tesBoundary(Fs2);
Fs2_b=Fs2(indBounary,:);

[indBounary]=tesBoundary(Fs1);
Fs1_b=Fs1(indBounary,:);

%%

cFigure;
hs=subplot(1,3,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0.25,3);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
ha=axis; axis off;

subplot(1,3,2);
title('Lattice order 1','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0,3);
gpatch(Fs1_b,Vs1,cMap(1,:),'none');
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

subplot(1,3,3);
title('Lattice order 2','fontSize',fontSize)
hold on;
gpatch(F,V,0.5*ones(1,3),'k',0,3);
gpatch(Fs2_b,Vs2,cMap(1,:),'none');
colormap(cMap);
caxis(cLim);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;
axis(ha); axis off;

drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [xi,w]=thirdOrder3DCubPoints()
%%THIRDORDER3DCUBPOINTS Generate third-order cubature points
%               for integration over a 3-dimensional cube with vertices of
%               (1,-1,-1), (-1,-1,-1), (-1,1,-1), (1,1,-1), (1,-1,1),
%               (-1,-1,1), (-1,1,1), and (1,1,1).
%
%INPUTS: None
%
%OUTPUTS: xi This is a 3XnumCubPoints set of points for the standard
%            cube.
%          w A 1XnumCubPoints set of cubature weights. This sums to the
%            volume of the standard cube (8).
%
%This function implements the points given in [1] (6 points).
%
%EXAMPLE:
%We compare a 2nd-order moment computed using these cubature points
%to one computed using monomialIntCube (a 3rd order moment would have
%just been 0). The results are the same within typical finite precision
%limits.
% [xi,w]=thirdOrder3DCubPoints();
% alpha=[2;0;0];
% theMoment=findMomentFromSamp(alpha,xi,w);
% intVal=monomialIntCube(alpha);
% RelErr=(theMoment-intVal)/intVal
%
%REFERENCES:
%[1] F. D. Witherden and P. E. Vincent, "On the identification of symmetric
%    quadrature rules for finite element methods," Computer and Mathematics
%    with Applications, vol. 69, no. 10, pp. 1232-1241, May 2015.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=[-1,  0,  0,  1.3333333333333333333333333333333333333;
    0,  0,  1,  1.3333333333333333333333333333333333333;
    0,  1,  0,  1.3333333333333333333333333333333333333;
    0,  0, -1,  1.3333333333333333333333333333333333333;
    1,  0,  0,  1.3333333333333333333333333333333333333;
    0, -1,  0,  1.3333333333333333333333333333333333333];

w=M(:,4);
xi=M(:,1:3)';

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function joint = xyzhumaneva2joint(pos,transf)

% XYZHUMANEVA2JOINT
%
% COPYRIGHT : Carl Henrik Ek and Neil Lawrence, 2008

% MOCAP

joint(:,1) = pos(1:3:end);
joint(:,2) = pos(2:3:end);
joint(:,3) = pos(3:3:end);

if(exist('transf','var'))
  joint = joint - repmat(transf.center,size(joint,1),1);
  M = rotationMatrix(transf.x,transf.y,transf.z);
  joint = joint*M;
  joint = joint + repmat(transf.center,size(joint,1),1);
end

return

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Fig. 6.6   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

num=1;
den=[1 0];
w=logspace(-2,3);
[m,p]=bode(num,den,w);;
m1=ones(50,1)./m;
den=[1 0 0];
[m2,p]=bode(num,den,w);
m2=ones(50,1)./m2;
loglog(w,m,w,m1,w,m2);
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.6 Magnitude of (j\omega)^n');
bodegrid;

author: Sable


text:
function [tests,pass,perf]=test_mat2wfdb(varargin)

  %Generate 3 different signals and convert them to signed 16 bit in WFDB format
str1=['N=1024;Fs=48000;tm=[0:1/Fs:(N-1)/Fs]'';adu=''V/mV/V'';' ...
    'info=''Example 1'';sig1=double(int16(linspace(-2^15,2^15,N)''));'...
    'sig2=double(int8(sin(2*pi*tm*1000).*(2^7)))./(2^7);' ...
    'sig3=(rand(N,1) > 0.97)*2 -1 + 2^16;sig=[sig1 sig2 sig3];' ...
    'mat2wfdb(sig,''Ex1'',Fs,[],adu,info);'];

cln1=['delete([pwd filesep ''Ex1*'']);'];

inputs={'verbose'};
verbose=0;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};']);
    end
end

%Test the examples
test_string={str1};
clean_up={cln1};

[tests,pass,perf]=test_wrapper(test_string,clean_up,verbose);

author: ikarosilva

repo: wfdb-app-toolbox