text:
% By NG Hoi-Fung, Ivan
% Discription: Convert the NMEA sentence to time, llh data
% Output CSV format: gps time, latitude, longitude, MSL(m), Geoid(m), utc time(yyyymmddHHMMSS.FFF)
% csv file will output to the original path that nmea file directory

clc;
clear;
close all;

GPS_TIME_START = datetime(1980,1,6,0,0,0,0);
leap_second = 18;

nmea_path = 'E:\sparse and urban jiulongtang midundao Huawei VINS test\20200314\COM55_200314_093715.ubx'; % input nmea sentence directory
[nmea_dir, exp_name] = fileparts(nmea_path);

support_type = {'$G*GGA','$G*RMC','$G*GLL','$G*BWC'};

% GGA: UTC, Lat, NE, Lon, NE, Fix Quality, No Sat, HDOP, MSL(m), Geoid(m)
% RMC: UTC, Status, Lat, NE, Lon, NE, Speed, Track angle, Date, Magnetic, checksum
% GLL: Lat, NE, Lon, NE, UTC, Date active, checksum
% BWC: UTC, Lat, NE, Lon, NE,bearing(true), bearing(magnetic), dist to waypoint(NM), waypoint id, checksum
% RMB: data status, Cross-track error(NM), Origin waypoint ID, Destination waypoint ID

date_str = {};
S = textread(nmea_path,'%s');
nmea_pos = double.empty(0,4);

for i = 1:length(S)
    type = regexp(S{i}, regexptranslate('wildcard', support_type));
    sup = 0;
    sup = find(~cellfun(@isempty,type));
    k = strfind(S{i},'$GNRMC');
    if k
       disp('');
    end
    if ~isempty(sup)
        utc = nan; lat = nan; lon = nan; msl = nan; geoid = nan;
        lat_ = nan; lon_ = nan; date_ = nan;
        lat_factor = 1; lon_factor = 1;

        str = S{i};
        str = str(type{sup}:end)
        C = strsplit(str,',', 'CollapseDelimiters',false);
        %         disp('');
        switch(sup)
            case 1 % GGA
                utc = str2double(C{2});
                utc_ = C{2};
                try
                    lat_ = str2double(C{3});
                    if C{4}~='N'; lat_factor = -1; end
                catch
                end
                try
                    lon_ = str2double(C{5});
                    if C{6}~='E'; lon_factor = -1; end
                catch
                end
                try msl = str2double(C{10}); catch, end
                try geoid = str2double(C{12}); catch, end
            case 2 % RMC
                utc = str2double(C{2});
                utc_ = C{2};
                try
                    lat_ = str2double(C{4});
                    if C{5}~='N'; lat_factor = -1; end
                catch
                end
                try
                    lon_ = str2double(C{6});
                    if C{7}~='E'; lon_factor = -1; end
                catch
                end
                try
                    datetime(C{10},'InputFormat','ddMMyy');
                    if isnan(date_); date_ = C{10}; end
                catch
                end
            case 3 % GLL
                utc = str2double(C{6});
                utc_ = C{6};
                try
                    lat_ = str2double(C{2});
                    if C{3}~='N'; lat_factor = -1; end
                catch
                end
                try
                    lon_ = str2double(C{4});
                    if C{5}~='E'; lon_factor = -1; end
                catch
                end
            case 4 % BWC
                utc = str2double(C{2});
                utc_ = C{2};
                try
                    lat_ = str2double(C{3});
                    if C{4}~='N'; lat_factor = -1; end
                catch
                end
                try
                    lon_ = str2double(C{5});
                    if C{6}~='E'; lon_factor = -1; end
                catch
                end
        end

        try
            datetime(utc_,'InputFormat','HHmmss.S');
        catch
            utc_ = nan;
        end

        if ~isnan(lat_)
            degree = floor(lat_/100);
            minute = floor(lat_ - degree*100)/60;
            second = (lat_ - (degree*100 + floor(lat_ - degree*100)))*60/3600;
            lat = degree + minute + second;
        end

        if ~isnan(lon_)
            degree = floor(lon_/100);
            minute = floor(lon_ - degree*100)/60;
            second = (lon_ - (degree*100 + floor(lon_ - degree*100)))*60/3600;
            lon = degree + minute + second;
        end

        if isempty(nmea_pos)
            nmea_pos = [utc,nan,nan,nan,nan];
        end
        utc_idx = find(nmea_pos(:,1)==utc);
        if isempty(utc_idx)
            utc_idx = size(nmea_pos,1)+1;
            nmea_pos(utc_idx,:) = [utc,nan,nan,nan,nan];
        end
        if isnan(nmea_pos(utc_idx,2)) && ~isnan(lat); nmea_pos(utc_idx,2)=lat; end
        if isnan(nmea_pos(utc_idx,3)) && ~isnan(lon); nmea_pos(utc_idx,3)=lon; end
        if isnan(nmea_pos(utc_idx,4)) && ~isnan(msl); nmea_pos(utc_idx,4)=msl; end
        if isnan(nmea_pos(utc_idx,5)) && ~isnan(geoid); nmea_pos(utc_idx,5)=geoid; end
        if ~isnan(date_); date_str(utc_idx,1) = {date_}; end
        if ~isnan(utc_); date_str(utc_idx,2) = {utc_}; end
    end
end

empty_date_idx = find(cellfun(@isempty,date_str(:,1)));
date_idx = find(~cellfun(@isempty,date_str(:,1)));
for i = 1:size(empty_date_idx,1)
    [~, store_idx] = min(date_idx - empty_date_idx(i));
    date_str(empty_date_idx(i),1) = date_str(date_idx(store_idx),1);
end
validIdx = find(~isnan(nmea_pos(:,1)));
nmea_pos = nmea_pos(validIdx,:);
date_str = date_str(validIdx,:);
gps_time = zeros(size(date_str,1),1);
utc_time = zeros(size(date_str,1),1);
for i = 1:size(gps_time,1)
    tempt = datetime(strcat(date_str{i,:}),'InputFormat','ddMMyyHHmmss.S');
    gps_time(i) = seconds(tempt - GPS_TIME_START) + leap_second;
    utc_time(i) = str2double(datestr(tempt, 'yyyymmddHHMMSS.FFF'));
end
nmea_pos(:,1) = gps_time;
nmea_pos = [nmea_pos utc_time];

% if exist('plot_openstreetmap'); plot_openstreetmap(nmea_pos(:,[2:3])); end
if exist('plot_osm'); plot_osm(nmea_pos(:,[2:3])); end

% Output format: gps time, latitude, longitude, MSL(m), Geoid(m), utc time(yyyymmddHHMMSS.FFF)
nmea_output = [nmea_dir,'\',exp_name,'_nmea.csv'];
fid = fopen(nmea_output,'w+');
for i = 1:size(nmea_pos,1)
    fprintf(fid,'%d,%.7f,%.7f,%.1f,%.1f,%.2f\n',nmea_pos(i,:));
end
fclose(fid);

ublox_data = csvread('COM55_200314_093715_nmea.csv');
%
kmlwrite('COM55_200314_093715_nmea.kml',ublox_data(:,2),ublox_data(:,3),'Icon',...
    'http://maps.google.com/mapfiles/kml/shapes/shaded_dot.png','IconScale',0.5,'Color',[1,0,0],'Name','  ');

author: weisongwen

repo: researchTools


text:
function i4vec_write ( output_filename, n, x )

%*****************************************************************************80
%
%% I4VEC_WRITE writes an I4VEC file.
%
%  Discussion:
%
%    An I4VEC is a vector of I4's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer N, the number of points.
%
%    Input, integer X(N), the vector.
%

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

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4VEC_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4VEC_WRITE - Error!' );
  end

  for j = 1 : n
    fprintf ( output_unit, '%d\n', x(j) );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [B,E,S,info] = DECOLOR(D,opt)
% DEteting Contiguous Outliers in the LOw-rank Representation
% http://arxiv.org/PS_cache/arxiv/pdf/1109/1109.0882v1.pdf
% eexwzhou@ust.hk
% Syntex: [B,S] = DECOLOR(D); or [B,S] = DECOLOR(D,opt);
% Input:
%   D -- 2D matrix
%   opt -- options. Usually, default setting is good. No need to specify.
%   opt.K: desired rank of the estimated low-rank component.
%          Default: \sqrt(min(size(D))) is good generally.
%   opt.lambda: a constant controls the strength of smoothness regularize
%               lambda ~ [1 5] is recommended. Default: 1
%   opt.sigma: STD of noise in the image. If not specified, computed online
%   opt.tol: convergence precision. Default: 1e-4
% Output:
%   B -- Low-rank component
%   S -- Outlier support
%   info -- other information

disp('^_^^_^^_^^_^^_^^_^ DECOLOR ^_^^_^^_^^_^^_^');
tic;

%% default parameter setting
if ~exist('opt','var'); opt = []; end
if ~isfield(opt,'tol'); opt.tol = 1e-4; end
if ~isfield(opt,'K'); opt.K = floor(sqrt(min(size(D)))); end
if ~isfield(opt,'lambda'); opt.lambda = 1; end % gamma = opt.lambda * beta;
if ~isfield(opt,'sigma'); opt.sigma = []; end % sigma can be estimated online

%% variable initialize
D = double(D);
B = D; % the low-rank matrix
S = false(size(D)); % background support
alpha = []; % Default setting by soft-impute
beta = 0.5*(std(D(:)))^2; % Start from a big value
minbeta = 0.5*(3*std(D(:))/20)^2; % lower bound: suppose SNR <= 20
sigma = opt.sigma; % if empty, will be estimated online
card = sum(S(:)); % record mumber of outliers
minCard = numel(D)/1e4; % minimum number of outliers
maxOuterIts = 50; % max number of iteration

% graph cuts initialization
% GCO toolbox is called
if opt.lambda > 0
    hMRF = GCO_Create(numel(D),2);
    GCO_SetSmoothCost( hMRF, [0 1;1 0] );
    AdjMatrix = getAdj(size(D));
    amplify = 10 * opt.lambda;
    GCO_SetNeighbors( hMRF, amplify * AdjMatrix );
end

%% outer loop
energy_old = inf; % total energy
for outerIts = 1:maxOuterIts
    disp(['---------------- Outer Loop:  ' num2str(outerIts) ' ----------------']);

    %% update B
    disp('*** Estimate Low-rank Matrix *** ');
    [B,Bnorm,alpha] = softImpute(D,B,~S,alpha,opt.K);
    E = D - B;

    %% estimate sigma
    if isempty(opt.sigma)
        sigma_old = sigma;
        sigma = std(E(~S(:)));
        if abs(sigma_old-sigma)/abs(sigma_old) < 0.01
            sigma = sigma_old; % if the change is not too large
        end
    end
    % update beta
    if card < minCard
        beta = beta/2;
    else
        beta = min(max([beta/2,0.5*(3*sigma)^2 minbeta]),beta);
    end
    gamma = opt.lambda * beta;

    %% estimate S
    disp('*** Estimate Outlier Support *** ');
    disp(['$$$ beta = ' num2str(beta) '; gamma = ' num2str(gamma) '; sigma = ' num2str(sigma)]);
    if opt.lambda > 0
        % call GCO to run graph cuts
        GCO_SetDataCost( hMRF, (amplify/gamma)*[ 0.5*(E(:)).^2, ones(numel(E),1)*beta]' );
        GCO_Expansion(hMRF);
        S = reshape(GCO_GetLabeling(hMRF)==2,size(S));
        card = sum(S(:)==1);
        energy_cut = (gamma/amplify)*double(GCO_ComputeEnergy(hMRF));
    else
        % direct hard thresholding if no smoothness
        S = 0.5*E.^2 > beta;
        card = sum(S(:));
        energy_cut = 0.5*norm(D-B-E,2)^2 + beta*card;
    end

    %% display energy
    energy = energy_cut + alpha * Bnorm;
    disp(['>>> the number of outliers is ' num2str(card)]);
    disp(['>>> the objectvive energy is ' num2str(energy)]);

    %% check termination condition
    if card > 0 && abs(energy_old-energy)/energy < opt.tol; break; end
    energy_old = energy;

end

info.opt = opt;
info.time = toc;
info.outerIts = outerIts;
info.energy = energy;
info.rank = rank(B);
info.alpha = alpha;
info.beta = beta;
info.sigma = sigma;

if opt.lambda > 0
    GCO_Delete(hMRF);
end

end

%% function to get the adjcent matirx of the graph
function W = getAdj(sizeData)
numSites = prod(sizeData);
id1 = [1:numSites, 1:numSites, 1:numSites];
id2 = [ 1+1:numSites+1,...
        1+sizeData(1):numSites+sizeData(1),...
        1+sizeData(1)*sizeData(2):numSites+sizeData(1)*sizeData(2)];
value = ones(1,3*numSites);
W = sparse(id1,id2,value);
W = W(1:numSites,1:numSites);
end

%% function for soft-impute
function [Z,Znorm,alpha] = softImpute(X,Z,Omega,alpha0,maxRank)
%
% This program implements the soft-impute algorithm followed by
% postprocessing in the Matrix completion paper Mazumder'10 IJML
% min || Z - X ||_Omega + \alpha || Z ||_Nulear
% \alpha is decrease from alpha0 to the minima value that makes rank(Z) <= maxRank

% X is the incomplete matrix
% maxRank is the desired rank in the constraint
% Omega is the mask with value 1 for data and 0 for missing part
if isempty(Z)
    Z = X;
end
if isempty(Omega)
    Omega = true(size(X));
end
if isempty(alpha0)
    [dummy,D] = svd(X,'econ');
    alpha0 = D(2,2);
end
if isempty(maxRank)
    maxRank = -1;
end
% parameters
eta = 0.707;
epsilon = 1e-4;
maxInnerIts = 20;
%% trivial
% no rank constraint
if maxRank >= min(size(X))
    Z = X;
    [dummy,D] = svd(Z,'econ');
    Znorm = sum(diag(D));
    alpha = 0;
    return;
end
% no observation
if sum(Omega(:)) == 0
    % no data
    Z = zeros(size(X));
    Znorm = 0;
    alpha = alpha0;
    return;
end
%% soft-impute
% 1. initialize
outIts = 0;
alpha = alpha0;
% 2. Do for alpha = alpha0 > alpha_1 > alpha_2 > ... > alpha_maxRank
disp('begin soft-impute iterations');
while 1
    outIts = outIts + 1;
    energy = inf;
    for innerIts = 1:maxInnerIts
        % (a)i
        C = X.*Omega + Z.*(1-Omega);
        [U,D,V] = svd(C,'econ');
        VT = V';
        % soft impute
        d = diag(D);
        idx = find(d > alpha);
        Z = U(:,idx) * diag( d(idx) - alpha ) * VT(idx,:);
        % (a)ii
        Znorm = sum(d(idx)-alpha);
        energy_old = energy;
        energy = alpha*Znorm + norm(Z(Omega(:))-X(Omega(:)),'fro')/2;
        if abs(energy - energy_old) / energy_old < epsilon
            break
        end
    end
    % check termination condition of alpha
    k = length(idx); % rank of Z
    disp(['alpha = ' num2str(alpha) ';    rank = ' num2str(k) ';  number of iteration: ' num2str(innerIts)]);
    if k <= maxRank && alpha > 1e-3
        alpha = alpha*eta;
    else
        break;
    end
end
end

author: andrewssobral

repo: lrslibrary


text:
function ptt = pulsetransit(ecgann, pulseann)
%   ptt = pulsetransit(ecgann, pulseann)
%
%   OVERVIEW:   Calculate PTT
%               match pulse annotations with corresponding ecg annotations
%
%   INPUT:
%
%   OUTPUT:
%
%   DEPENDENCIES & LIBRARIES:
%   REFERENCE:
%	REPO:
%       https://github.com/cliffordlab/hrv_toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Main script written by Adriana N. Vest
%       Dependent scripts written by various authors
%       (see functions for details)
%	COPYRIGHT (C) 2016
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%%

ptt = NaN((length(ecgann)-1),3);
for i = 1:(length(ecgann)-1)
    x = find(pulseann < ecgann(i+1) & pulseann > ecgann(i));
    ptt(i,1) = ecgann(i);
    if x > 0
        if length(x)>1
            ptt(i,2) = pulseann(x(1));
        else
            ptt(i,2) = pulseann(x);
        end
    end
    ptt(i,3) = ptt(i,2) - ptt(i,1);
end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function whelp(classpath)
% Display help about the available options to the console.
%
% classpath - the full path to the class omiting 'weka'.
%
% Example: whelp classifiers.bayes.BayesNet

    eval(['weka.',classpath,'.main(''-h'')']);
end

author: Sable

repo: mcbench-benchmarks


text:
function [r,xp] = spm_dcm_erp_bma (BMS_file,stats,params)
% Compute posterior over connections, or modulatory gains in them, from BMA
% FORMAT [r,xp] = spm_dcm_erp_bma (BMS_file,stats,params)
%
% BMS_file      Name of Bayesian Model Selection .mat file
% stats         'ffx' or 'rfx'
% params        Parameter data structure
%               .type          'A' (connection) or 'B' (gain in connection)
%               .hier          'forward', 'backward' or 'lateral' (for conn_type='A')
%               .ip            eg 1, 2, 3 indexes modulatory input (for conn_type='B')
%               .to            to region eg 3
%               .from          from region eg 1
%               .xt            exceedance threshold (typically set to 1)
%               .C             [nr x nr] contrast matrix where nr is the number of regions
%
%
% r             posterior samples
% xp            exceedance probability
%               This is the posterior probability that the connection is
%               larger than params.xt. Alternatively, if you are looking at
%               a contrast of connections, its the posterior probability
%               that the contrast is greater than zero.
%
% The parameters returned by Bayesian Model Averaging (BMA) are the 'latent'
% variables A and B which are Gaussian (and consequently can be positive or
% negative).
%
% The corresponding connection strengths (rA) or gains in connection
% strength (rB) are an exponential function of these latent variables.
% These are the values we are interested in and want to make an inference
% about.
%
% This routine computes the posterior distribution over rA or rB by
% generating samples from the latent variables, and exponentiating each
% sample.
%
% The probability that the rA or RB values are greater than some threshold
% xt (such as unity) is then just the proportion of posterior samples that
% are greater than xt.
%
% If a contrast matrix (C) is not specifed this function looks at a single
% connection or gain. To look at relative sizes of connection/gain values
% enter a C matrix. eg. to test, in a 3 region DCM, is connection from 3
% to 2 bigger than 2 to 3 ? set C=[0 0 0; 0 0 1; 0 -1 0].
%
%--------------------------------------------------------------------------
%
% Example usage:
%
% 1. Look at a single connection value:
%
% params.type='A'; params.hier='forward';
% params.to=3; params.from=1; params.xt=1;
% spm_dcm_erp_bma([],'ffx',params);
%
% 2. Look at a single gain value:
%
% params.type='B'; params.ip=1;
% params.to=1; params.from=1; params.xt=1;
% spm_dcm_erp_bma([],'ffx',params);
%
% 3. Look at a contrast of connection values:
%
% params.type='B'; params.ip=1;
% params.C=[0 0 0; 0 0 1; 0 -1 0];
% spm_dcm_erp_bma([],'ffx',params);
%__________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_dcm_erp_bma.m 4529 2011-10-18 14:58:20Z guillaume $

%-Get parameters
%--------------------------------------------------------------------------
try
    load (BMS_file);
catch
    [BMS_file,sts] = spm_select(1,'^BMS\.mat$','Select BMS.mat file');
    if ~sts, return; end
    load (BMS_file);
end

if isfield(params,'C')
    con = 1;
else
    con = 0;
end

switch lower(stats)
    case 'ffx'
        bma = BMS.DCM.ffx.bma;
    case 'rfx'
        bma = BMS.DCM.rfx.bma;
    otherwise
        error('The stats field shout be set to ''ffx'' or ''rfx''');
end

%-Get mean and SD of effects
%--------------------------------------------------------------------------
switch params.type
    case 'A'
        hier_type = {'forward','backward','lateral'};
        h = find(ismember(hier_type,params.hier));
        M = bma.mEp.A{h};
        S = bma.sEp.A{h};
        if ~con
            str=sprintf('%s connection (rA) from region %d to %d:',hier_type{h},params.from,params.to);
        end
    case 'B'
        M = bma.mEp.B{params.ip};
        S = bma.sEp.B{params.ip};
        if ~con
            str=sprintf('Effect (rB) of modulatory input %d on connection from region %d to %d:',params.ip,params.from,params.to);
        end
    otherwise
        error('conn_type must be ''A'' or ''B''');
end

%-Generate N samples in latent space, then exponentiate to get estimate
% of effects
%--------------------------------------------------------------------------
N=10000;

switch con
    case 1
        % Look at contrast of connections
        nr = size(M,1);
        m  = M(:);
        s  = S(:);
        c  = params.C(:);

        % Full posterior covariances are not stored
        % only univariate SD's so use these
        x  = (s*ones(1,N)).*randn(nr*nr,N)+m*ones(1,N);
        r  = c'*exp(x);
        rm = mean(r);
        xp = length(find(r>0))/length(r);
        switch params.type
            case 'A'
                str=sprintf('Contrast of %s connections (rA)',hier_type{h});
            case 'B'
                str=sprintf('Contrast of effects (rB) of modulatory input %d ',params.ip);
        end
    case 0
        % Just look at a single connection
        m  = M(params.to,params.from);
        s  = S(params.to,params.from);
        x  = s*randn(N,1)+m;
        r  = exp(x);
        rm = mean(r);
        xp = length(find(r>params.xt))/length(r);
end

%-Display
%--------------------------------------------------------------------------
figure
hist(r,20);
set(gca,'FontSize',18);
set(gca,'YTickLabel',[]);
ylabel(sprintf('p(r%s|Y)',params.type));
xlabel(sprintf('r%s',params.type));

%-Report
%--------------------------------------------------------------------------
fprintf('%s\n',str);
fprintf('Posterior mean = %1.2f\n',rm);
fprintf('Exceedance probability = %1.6f\n', xp);

author: spm


text:
function pred_value = xval_featureselect_nmdscluster(X, Y, pthreshold, holdout_method, xyzlist)
%
% pred_value = xval_featureselect_nmdscluster(X, Y, p-value selection for univariate feature selection, holdout_method)
%
% pred_value: cross-validated predictions of outcome data
% X: n x variables matrix of predictors
% Y: n x 1 vector of outcomes
%
% Tor Wager, June 2010
%
% Go to any LASSO output directory and run this:
%
% maskInfo = iimg_read_img(fullfile(pwd, 'mask.img'), 2);
% dat{1} =  iimg_get_data(maskInfo, anticimages);
% pred_value = xval_simple_ols_loo(X, Y)

[N, k] = size(X);
pred_value = zeros(N, 1);

holdout_set = nested_select_holdout_set;

create_figure('test', 1, 2);

fprintf('Fold: %03d', 0);
for wh_fold = 1:length(holdout_set)

    fprintf('\b\b\b%3.0f ', wh_fold);

    % select training data
    Xi = X;
    Yi = Y;

    Yi(holdout_set{wh_fold}) = [];
    Xi(holdout_set{wh_fold}, :) = [];              % leave out the missing observation(s)

    Xtest = X(holdout_set{wh_fold}, :);
    Ytest = Y(holdout_set{wh_fold}, :);          % only used later when we test prediction

    ntrain = size(Xi, 1);
    nholdout = size(Xtest, 1);

    % select features based on univariate correlations
    [r, p, Tstat] = correlation_fast_series(Xi, Yi);
    wh_features = p <= pthreshold;
    nfeatures(wh_fold) = sum(wh_features);
    if nfeatures(wh_fold) == 0
        disp('Warning: no features pass threshold');
        wh_features = p <= prctile(p, 10);
    end

    Xi = Xi(:, wh_features);

    c = []; c.outcome = Yi; c.covs_nointerest = [];
    % get average data from contiguous blobs (??)
    xyz = xyzlist(wh_features, :);
    contig = spm_clusters(xyzlist(wh_features, :)');
    for ii = 1:max(contig)
        if sum(contig == ii) < 3, contig(contig == ii) = 0; end
    end
    u = unique(contig);    u(u == 0) = [];
    for ii = 1:length(u)
        c.dat(:, ii) = nanmean(Xi(:, contig == u(ii))')';
    end
    c = nmdsfig_tools('get_correlations',c);
    c.GroupSpace = mdscale(c.D, 8, 'Replicates', 10,'Options',statset('MaxIter',500),'Criterion','sstress');
    c.ClusterSolution.classes = clusterdata(c.GroupSpace, 'maxclust', 5, 'linkage', 'average');

    c.APPLY_CLUSTER = apply_cluster_solution(c.ClusterSolution.classes,...
        c.r,...
        'names',[],'bcov',c.outcome, 'n', size(c.dat, 1), 'dointeractive', 0);
    class_avg_dat = [];
     for i = 1:max(c.ClusterSolution.classes)
        wh = find(c.ClusterSolution.classes == i);
        if length(wh) > 1
            class_avg_dat = [class_avg_dat nanmean(c.dat(:,wh)')'];
        else
            class_avg_dat = [class_avg_dat c.dat(:,wh)];
        end
     end
    c.class_avg_dat = class_avg_dat;

    c.class_STEPWISE = stepwise_tor(c.class_avg_dat,c.outcome,c.APPLY_CLUSTER.class_names);

    Xi = [ones(ntrain, 1) c.class_avg_dat(:, c.class_STEPWISE.inmodel)]; % add intercept

    % apply to Xtest
    c.testdat = zeros(size(Xtest, 1), size(c.dat, 2));
    for ii = 1:length(u)
        c.testdat(:, ii) = nanmean(Xtest(:, contig == u(ii))')';
    end
     class_avg_dat = [];
     for i = 1:max(c.ClusterSolution.classes)
        wh = find(c.ClusterSolution.classes == i);
        if length(wh) > 1
            class_avg_dat = [class_avg_dat nanmean(c.testdat(:,wh)')'];
        else
            class_avg_dat = [class_avg_dat c.testdat(:,wh)];
        end
     end
    c.class_avg_testdat = class_avg_dat;
    Xtest = c.class_avg_testdat(:, c.class_STEPWISE.inmodel);

    % Make prediction
    b = pinv(Xi) * Yi;
    pred_value(holdout_set{wh_fold}, 1) = [ones(nholdout, 1) Xtest] * b;

    create_figure('test', 1, 2, 1); subplot(1, 2, 1); plot(Xi*b, Yi, 'ko');
    plot(pred_value(holdout_set{wh_fold}, 1), Ytest, 'ro', 'MarkerFaceColor', 'r');
    subplot(1, 2, 2);
   title('Black circles = training set; Red = holdout obs');
    drawnow;

end

cm = colormap(jet(N));
figure; hold on;
for i = 1:N
    plot(pred_value(i), Y(i), 'ko', 'MarkerFaceColor', cm(i, :));
end
xlabel('Predicted outcome (xval)'); ylabel('Outcome');
title('Color = order in data series');

function holdout_set = nested_select_holdout_set
        % purpose: return holdout_set variable

        nobs_s = length(Y);

        switch lower(holdout_method)
            case 'loo'
                holdout_set = cell(1, nobs_s);
                for i = 1:nobs_s, holdout_set{i} = i; end

            case 'l4o_covbalanced'
                disp('Selecting holdout sets: Balancing on covariates and outcome, and also trying to ensure that each obs is selected equally often.');
                holdout_proportion = 4 ./ nobs_s; % prop for leave-4-out
                nfolds = 40;
                if isempty(cov_val)
                    wh_holdout = xval_select_holdout_set(Y, [], nfolds, holdout_proportion, verbose);
                else
                    wh_holdout = xval_select_holdout_set(Y, cov_val, nfolds, holdout_proportion, verbose);
                end
                holdout_set = cell(1, nfolds);
                for k = 1:nfolds
                    holdout_set{k} = wh_holdout(:, k);
                end

            case 'categorical_covs'

                holdout_set = xval_select_holdout_set_categoricalcovs(cov_val);

            case 'balanced4'
                nfolds = 4;
                holdout_set = cell(1, nfolds);
                [ys, indx] = sort(Y);
                for k = 1:nfolds

                    holdout_set{k} = indx(k:nfolds:end);
                    if isempty(holdout_set{k}), error('Holdout set construction error: Dataset too small?'); end

                end

            otherwise error('Unknown holdout method.  See help.');
        end
end

end % main function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function varargout = virtualsignallab(varargin)
global choice
% VIRTUALSIGNALLAB M-file for virtualsignallab.fig
%      VIRTUALSIGNALLAB, by itself, creates a new VIRTUALSIGNALLAB or raises the existing
%      singleton*.
%
%      H = VIRTUALSIGNALLAB returns the handle to a new VIRTUALSIGNALLAB or the handle to
%      the existing singleton*.
%
%      VIRTUALSIGNALLAB('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in VIRTUALSIGNALLAB.M with the given input arguments.
%
%      VIRTUALSIGNALLAB('Property','Value',...) creates a new VIRTUALSIGNALLAB or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before virtualsignallab_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to virtualsignallab_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES
% Edit the above text to modify the response to help virtualsignallab
% Last Modified by GUIDE v2.5 02-Feb-2005 16:15:56
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @virtualsignallab_OpeningFcn, ...
                   'gui_OutputFcn',  @virtualsignallab_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
% --- Executes just before virtualsignallab is made visible.
function virtualsignallab_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to virtualsignallab (see VARARGIN)
% Choose default command line output for virtualsignallab
load train;
tr=y;
handles.tr=tr;
load handel;
ha=y;
handles.ha=ha;
load laughter;
la=y;
handles.la=la;
load chirp;
ch=y;
handles.ch=ch;
d=size(handles.tr,1);
h_f=findobj(0,'tag','Fs');
m=get(h_f,'string');
n=str2num(m);
handles.n=n;
handles.x=0;
handles.y=0;
handles.current_data=0;
handles.col='  blanc   ';
handles.dia=' no domain ';
handles.string=' no ';
choice=[1 1 1]
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes virtualsignallab wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = virtualsignallab_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
varargout{1} = handles.output;
% --- Executes during object creation, after setting all properties.
function soundmenu_CreateFcn(hObject, eventdata, handles)
% hObject    handle to soundmenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on selection change in soundmenu.
function soundmenu_Callback(hObject, eventdata, handles)
val=get(hObject,'value');
str=get(hObject,'string');
switch str{val}
    case 'please select'
        errordlg('you must select a field entery','bad input','modal');
    case 'train'
        handles.current_data=handles.tr;
        handles.string='train';
    case 'handel'
        handles.current_data=handles.ha;
        handles.string='handel';
    case 'laughter'
         handles.current_data=handles.la;
         handles.string='laughter';
    case 'chirp'
        handles.current_data=handles.ch;
        handles.string='chirp';
    otherwise
    handles.current_data=0;
end
    guidata(hObject,handles);
% hObject    handle to soundmenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Hints: contents = get(hObject,'String') returns soundmenu contents as cell array
%        contents{get(hObject,'Value')} returns selected item from
%        soundmenu
% --- Executes during object creation, after setting all properties.
function diagrammenu_CreateFcn(hObject, eventdata, handles)
% hObject    handle to diagrammenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on selection change in diagrammenu.
function diagrammenu_Callback(hObject, eventdata, handles)
% hObject    handle to diagrammenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns diagrammenu contents as cell array
%        contents{get(hObject,'Value')} returns selected item from
%        diagrammenu
n=handles.n;
val=get(hObject,'value');
str=get(hObject,'string');
switch str{val}
    case 'please select'
        y=handles.current_data;
        x=0;
        errordlg('you must select a field entery','bad input','modal');
    case 'time domain'
      y=handles.current_data;
      d=size(y,1);
      T=d/n;
      x=linspace(0,T,d);
      handles.dia='time domain';
    case 'spectrum'
      y=handles.current_data;
      d=size(y,1);
      q=1/d*fft(y);
      real=abs(q);
      y=fftshift(real);
      x=linspace(-n/2,n/2,d);
       handles.dia='frequency domain';
    case 'periodogram'
      y=handles.current_data;
      d=size(y,1);
      q=1/d*fft(y);
      real=abs(q);
      y=fftshift(real.^2);
      x=linspace(-n/2,n/2,d);
       handles.dia='frequency domain';
  otherwise
      y=handles.current_data;
      x=0;
end
      handles.y=y;
      handles.x=x;
     guidata(hObject,handles);
% --- Executes during object creation, after setting all properties.
function colormenu_CreateFcn(hObject, eventdata, handles)
% hObject    handle to colormenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on selection change in colormenu.
function colormenu_Callback(hObject, eventdata, handles)
global choice;
h_c=findobj(0,'tag','colormenu');
val=get(h_c,'value');
str=get(h_c,'string');
switch str{val}
    case 'please select'
        errordlg('you must select a field entery','bad input','modal');
    case 'blue'
        choice=[0 0 1];
        handles.col='blue';
    case 'red'
        choice=[1 0 0];
        handles.col='red';
    case 'green'
        choice=[0 1 0];
        handles.col='green';
    otherwise
        choice=[1 1 1];
end
guidata(hObject,handles);
% hObject    handle to colormenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Hints: contents = get(hObject,'String') returns colormenu contents as cell array
%        contents{get(hObject,'Value')} returns selected item from
%        colormenu
% --- Executes on button press in mute.
function mute_Callback(hObject, eventdata, handles)
% hObject    handle to mute (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of mute
% --- Executes on button press in play.
function play_Callback(hObject, eventdata, handles)
% hObject    handle to play (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Hint: get(hObject,'Value') returns toggle state of play
% --- Executes during object creation, after setting all properties.
sound(handles.current_data,handles.n);
guidata(hObject,handles);
function Fs_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Fs (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
function Fs_Callback(hObject, eventdata, handles)
h_f=findobj(0,'tag','Fs');
m=get(h_f,'string');
n=str2num(m);
handles.n=n;
guidata(hObject,handles);
% hObject    handle to Fs (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of Fs as text
%        str2double(get(hObject,'String')) returns contents of Fs as a
%        double
% --- Executes on button press in ok.
function ok_Callback(hObject, eventdata, handles)
global choice
x=handles.x;
y=handles.y;
h_p=plot(x,y);
title([' A plot of the  ',handles.string,' signal with the  ',handles.dia,' in ',handles.col ,'color']);
set(h_p,'color',choice);
sro=findobj(0,'tag','soundmenu');
set(sro,'value',1);
sra=findobj(0,'tag','diagrammenu');
set(sra,'value',1);
srb=findobj(0,'tag','colormenu');
set(srb,'value',1);
guidata(hObject,handles);
% hObject    handle to ok (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% --- Executes on button press in PLAY.


text:
X=randn(64,200);
Y=randn(200,20000)';

tic
XYt=mexCalcXYt(X,Y);
t=toc;
fprintf('mex-file time: %fs\n',t);

tic
XYt2=X*Y';
t=toc;
fprintf('matlab-file time: %fs\n',t);

sum((XYt(:)-XYt2(:)).^2)

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
N = 15; % Set N as an odd number;
w0 = 0.4*pi;
n0 = floor(N/2);
hn = zeros(1,N);
for i = 1:length(hn)
    if i == n0
        hn(i) = w0/pi;
    else
        hn(i) = sin(w0*(i-n0))/(pi*(i-n0));
    end;
end;
subplot(3,1,1);stem(hn);title('h[n]');

Np = 1000;
hn = [hn,zeros(1, Np-length(hn))];%zero padding
[w,y] = calculateDiscreteFourierTransform(hn);

subplot(3,1,2);plot(w,abs(y)),title('幅度谱');
subplot(3,1,3);plot(w,angle(y)),title('相位谱');

author: VipaiLab

repo: Signals-and-Systems-course


text:
function Write_PhaseMap(FileWritingSetup)
% write_PhaseMap(FileWritingSetup)
%
% FUNCTION:
%
%   This function takes the symbolically derived impact map equations and
%   writes them to a function file.
%
%
% ARGUMENTS: FileWritingSetup.
%
%   SolnOne = A struct with a field containing an expression for the rates
%             immediately after a sticking collision with foot one
%
%   SolnTwo = A struct with a field containing an expression for the rates
%             immediately after a sticking collision with foot two
%
%   States = [2N x 2] cell array where the first column is the name of
%       the states, and the second column is comments about that states.
%       The index of a state in this array will match the index in the
%       dynamics function to be written. It is assumed that the second
%       block of N states are the time derivatives of the first block of N
%       states. Furthermore, it is assumed that the first N states are
%       positions and are not affected by impact.
%
%   Parameters = [K x 2] cell array with all of the input parameters. The
%       first column is for names, the second for comments
%
%   Directory = the name of the directory to write the files in.
%
%
% RETURNS:
%   This function will write a file for the dynamics equations in each
%   field of Dyn.
%
%
% NOTE:
%
%   1)  It is assumed that the symbols in all of the inputs are self
%       consistent; that is - if Contacts{1,1} = H1, then there must be a field
%       in Dyn.(FileNames{i,1}) called H1.
%
%
% Written by Matthew Kelly
% December 6, 2013
% Cornell University
%
% See also DERIVE_EOM

States = FileWritingSetup.States;
SolnOne = FileWritingSetup.SolnOne;
SolnTwo = FileWritingSetup.SolnTwo;
Directory = FileWritingSetup.Directory;

CurrDir = cd;
if ~exist(Directory,'dir')
    mkdir(Directory);
end
cd(Directory);

fid = fopen('phaseMap.m','w');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                           Function Header                               %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

fprintf(fid,['function StateAfter = phaseMap(StateBefore,PhaseBefore,PhaseAfter)\n']);
fprintf(fid,'%% function StateAfter = phaseMap(StateBefore,PhaseBefore,PhaseAfter)\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% Computer Generated File -- DO NOT EDIT\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% This function was created by the function Write_PhaseMap()\n');
fprintf(fid,['%% ' datestr(now) '\n']);
fprintf(fid,'%%\n');
fprintf(fid,'%% Dymanics Model: retractable double pendulum biped\n');
fprintf(fid,'%% Phases: Flight (F) == Both feet in the air\n');
fprintf(fid,'%%         Double Stance (D) == Both feet on the ground\n');
fprintf(fid,'%%         Single Stance One (S1) == Foot One -> ground, Foot Two -> air\n');
fprintf(fid,'%%         Single Stance Two (S2) == Foot One -> air, Foot Two -> ground\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% FUNCTION:\n');
fprintf(fid,'%%   This function will compute the impact map that occurs during\n');
fprintf(fid,'%%   transitions between any two phases of motion. Detailed assumptions can\n');
fprintf(fid,'%%   be found in Derive_EoM.m\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% ARGUMENTS:\n');
fprintf(fid,'%%   StateBefore = [Nx1] state vector\n');
fprintf(fid,'%%   PhaseBefore = element of set {''F'',''D'',''S1'',''S2''} = phase before impact\n');
fprintf(fid,'%%   PhaseAfter = element of set {''F'',''D'',''S1'',''S2''} = phase after impact\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% RETURNS:\n');
fprintf(fid,'%%   StateAfter = [Nx1] state vector\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% Matthew Kelly\n');
fprintf(fid,'%% Cornell University\n');
fprintf(fid,'%%\n');
fprintf(fid,'%%\n');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                            Select Impact Map                            %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

fprintf(fid,'switch PhaseBefore    \n');
fprintf(fid,'    case ''F''    \n');
fprintf(fid,'        switch PhaseAfter    \n');
fprintf(fid,'            case ''F''    \n');
fprintf(fid,'                StateAfter = StateBefore;    \n');
fprintf(fid,'            case ''D'' %%Assume that collision happens on Foot One then Two    \n');
fprintf(fid,'                StateInt = MapOne(StateBefore);    \n');
fprintf(fid,'                StateAfter = MapTwo(StateInt);    \n');
fprintf(fid,'            case ''S1''    \n');
fprintf(fid,'                StateAfter = MapOne(StateBefore);    \n');
fprintf(fid,'            case ''S2''    \n');
fprintf(fid,'                StateAfter = MapTwo(StateBefore);    \n');
fprintf(fid,'            otherwise    \n');
fprintf(fid,'                error(''Invalid phase after collision. Must be element of set: {''''F'''',''''D'''',''''S1'''',''''S2''''}'');    \n');
fprintf(fid,'        end    \n');
fprintf(fid,'    case ''D''    \n');
fprintf(fid,'         StateAfter = StateBefore;    \n');
fprintf(fid,'    case ''S1''    \n');
fprintf(fid,'        switch PhaseAfter    \n');
fprintf(fid,'            case {''F'',''S1''}    \n');
fprintf(fid,'                StateAfter = StateBefore;    \n');
fprintf(fid,'            case {''D'',''S2''}    \n');
fprintf(fid,'                StateAfter = MapTwo(StateBefore);    \n');
fprintf(fid,'            otherwise    \n');
fprintf(fid,'                error(''Invalid phase after collision. Must be element of set: {''''F'''',''''D'''',''''S1'''',''''S2''''}'');    \n');
fprintf(fid,'        end    \n');
fprintf(fid,'    case ''S2''    \n');
fprintf(fid,'        switch PhaseAfter    \n');
fprintf(fid,'            case {''F'',''S2''}    \n');
fprintf(fid,'                StateAfter = StateBefore;    \n');
fprintf(fid,'            case {''D'',''S1''}    \n');
fprintf(fid,'                StateAfter = MapOne(StateBefore);    \n');
fprintf(fid,'            otherwise    \n');
fprintf(fid,'                error(''Invalid phase after collision. Must be element of set: {''''F'''',''''D'''',''''S1'''',''''S2''''}'');    \n');
fprintf(fid,'        end    \n');
fprintf(fid,'    otherwise    \n');
fprintf(fid,'                error(''Invalid phase before collision. Must be element of set: {''''F'''',''''D'''',''''S1'''',''''S2''''}'');    \n');
fprintf(fid,'end    \n');
fprintf(fid,'\n');
fprintf(fid,'end\n');
fprintf(fid,'\n\n');
fprintf(fid,'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n');
fprintf(fid,'%%%%%%%%  ~~~~~~~~~~~~~~~~~~~ Sub Functions  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%%%%%\n');
fprintf(fid,'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n');
fprintf(fid,'\n\n');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                        Write MapOne subfunction                         %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

    fprintf(fid,'function StateAfter = MapOne(StateBefore)    \n\n');

    for idxState=1:size(States,1)
        fprintf(fid,[States{idxState,1} ' = StateBefore(' num2str(idxState) '); %% ' ...
            States{idxState,2} '\n']);
    end %idxState
    fprintf(fid,'\n');

    Ndof = size(States,1)/2;
    N = num2str(Ndof);
    fprintf(fid,'StateAfter = zeros(size(StateBefore));\n');
    fprintf(fid,['StateAfter(1:' N ',:) = '...
        'StateBefore(1:' N ',:); %% Configuration not changed by collision\n']);
    for idx=1:Ndof
        writeIndex = idx + Ndof;
        writeName = [States{writeIndex,1} '_plus'];
        vectorizedExpression = char(SolnOne.(writeName));
        fprintf(fid,['StateAfter(' num2str(writeIndex) ') = ' ...
            vectorizedExpression ';\n']);
    end %idxAcc
    fprintf(fid,'\n');

    fprintf(fid,'end\n\n\n');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                        Write MapTwo subfunction                         %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

    fprintf(fid,'function StateAfter = MapTwo(StateBefore)    \n\n');

    for idxState=1:size(States,1)
        fprintf(fid,[States{idxState,1} ' = StateBefore(' num2str(idxState) '); %% ' ...
            States{idxState,2} '\n']);
    end %idxState
    fprintf(fid,'\n');

    Ndof = size(States,1)/2;
    N = num2str(Ndof);
    fprintf(fid,'StateAfter = zeros(size(StateBefore));\n');
    fprintf(fid,['StateAfter(1:' N ',:) = '...
        'StateBefore(1:' N ',:); %% Configuration not changed by collision\n']);
    for idx=1:Ndof
        writeIndex = idx + Ndof;
        writeName = [States{writeIndex,1} '_plus'];
        vectorizedExpression = char(SolnTwo.(writeName));
        fprintf(fid,['StateAfter(' num2str(writeIndex) ') = ' ...
            vectorizedExpression ';\n']);
    end %idxAcc
    fprintf(fid,'\n');

    fprintf(fid,'end\n');

fclose(fid);

cd(CurrDir);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
% tHis subroutine assigns creates a grid with
% spacing dx,dy (in degreees). The size will
% be selected interactiVELY. The bvalue in each
% volume around a grid point containing ni earthquakes
% will be calculated as well as the magnitude
% of completness
%   Stefan Wiemer 1/95

report_this_filefun(mfilename('fullpath'));

global no1 bo1 inb1 inb2

if sel == 'in'
    % get the grid parameter
    % initial values
    %
    dd = 1.00;
    dx = 1.00 ;
    ni = 100;
    ra = 5;

    % make the interface
    %
    figure_w_normalized_uicontrolunits(...
        'Name','Grid Input Parameter',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'units','points',...
        'Visible','off', ...
        'Position',[ wex+200 wey-200 550 300]);
    axis off

    labelList2=['Weighted LS - automatic Mcomp | Weighted LS - no automatic Mcomp '];
    labelPos = [0.2 0.7  0.6  0.08];
    hndl2=uicontrol(...
        'Style','popup',...
        'Position',labelPos,...
        'Units','normalized',...
        'String',labelList2,...
        'Callback','inb2 =get(hndl2,''Value''); ');

    labelList=['Maximum likelihood - automatic Mcomp | Maximum likelihood  - no automatic Mcomp '];
    labelPos = [0.2 0.8  0.6  0.08];
    hndl1=uicontrol(...
        'Style','popup',...
        'Position',labelPos,...
        'Units','normalized',...
        'String',labelList,...
        'Callback','inb1 =get(hndl1,''Value''); ');

    % creates a dialog box to input grid parameters
    %
    freq_field=uicontrol('Style','edit',...
        'Position',[.60 .50 .22 .10],...
        'Units','normalized','String',num2str(ra),...
        'Callback','ra=str2double(get(freq_field,''String'')); set(freq_field,''String'',num2str(ra));');

    freq_field2=uicontrol('Style','edit',...
        'Position',[.60 .40 .22 .10],...
        'Units','normalized','String',num2str(dx),...
        'Callback','dx=str2double(get(freq_field2,''String'')); set(freq_field2,''String'',num2str(dx));');

    freq_field3=uicontrol('Style','edit',...
        'Position',[.60 .30 .22 .10],...
        'Units','normalized','String',num2str(dd),...
        'Callback','dd=str2double(get(freq_field3,''String'')); set(freq_field3,''String'',num2str(dd));');

    close_button=uicontrol('Style','Pushbutton',...
        'Position',[.60 .05 .15 .12 ],...
        'Units','normalized','Callback','close;done','String','Cancel');

    go_button1=uicontrol('Style','Pushbutton',...
        'Position',[.20 .05 .15 .12 ],...
        'Units','normalized',...
        'Callback','inb1 =get(hndl1,''Value'');inb2 =get(hndl2,''Value'');close,sel =''ca'', bcrossVt',...
        'String','Go');

    text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.20 1.0 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.l ,...
        'FontWeight','bold',...
        'String','Automatically estimate magn. of completeness?   ');

    txt3 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.30 0.65 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.l ,...
        'FontWeight','bold',...
        'String',' Grid Parameter');
    txt5 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0. 0.42 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing along projection [km]');

    txt6 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0. 0.32 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Spacing in depth in km:');

    txt1 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0. 0.53 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m,...
        'FontWeight','bold',...
        'String','Radius in km');

    set(gcf,'visible','on');
    watchoff

end   % if sel == in

% get the grid-size interactively and
% calculate the b-value in the grid by sorting
% thge seimicity and selectiong the ni neighbors
% to each grid point

if sel == 'ca'

    figure_w_normalized_uicontrolunits(xsec_fig)
    hold on

    messtext=...
        ['To select a polygon for a grid.       '
        'Please use the LEFT mouse button of   '
        'or the cursor to the select the poly- '
        'gon. Use the RIGTH mouse button for   '
        'the final point.                      '
        'Mac Users: Use the keyboard "p" more  '
        'point to select, "l" last point.      '
        '                                      '];

    zmap_message_center.set_message('Select Polygon for a grid',messtext);

    hold on
ax = findobj('Tag','main_map_ax');
[x,y, mouse_points_overlay] = select_polygon(ax);
    zmap_message_center.set_info('Message',' Thank you .... ')

    plos2 = plot(x,y,'b-','era','xor');        % plot outline
    sum3 = 0.;
    pause(0.3)

    %create a rectangular grid
    xvect=[min(x):dx:max(x)];
    yvect=[min(y):dd:max(y)];
    gx = xvect;gy = yvect;
    tmpgri=zeros((length(xvect)*length(yvect)),2);
    n=0;
    for i=1:length(xvect)
        for j=1:length(yvect)
            n=n+1;
            tmpgri(n,:)=[xvect(i) yvect(j)];
        end
    end
    %extract all gridpoints in chosen polygon
    XI=tmpgri(:,1);
    YI=tmpgri(:,2);

    ll = polygon_filter(x,y, XI, YI, 'inside');
    %grid points in polygon
    newgri=tmpgri(ll,:);

    % Plot all grid points
    plot(newgri(:,1),newgri(:,2),'+k')

    if length(xvect) < 2 || length(yvect) < 2
        errordlg('Selection too small! (not a matrix)');
        return
    end

    itotal = length(newgri(:,1));

    zmap_message_center.set_info(' ','Running... ');think
    %  make grid, calculate start- endtime etc.  ...
    %
    t0b = min(newa.Date)  ;
    n = newa.Count;
    teb = newa(n,3) ;
    tdiff = round((teb - t0b)*365/par1);
    loc = zeros(3, length(gx)*length(gy));

    % loop over  all points
    %
    i2 = 0.;
    i1 = 0.;
    bvg = [];
    allcount = 0.;
    wai = waitbar(0,' Please Wait ...  ');
    set(wai,'NumberTitle','off','Name','b-value grid - percent done');;
    drawnow
    %
    % loop

    % overall b-value
    [bv magco stan av me mer me2,  pr] =  bvalca3(newa,inb1,inb2);
    bo1 = bv; no1 = newa.Count;
    %
    for i= 1:length(newgri(:,1))
        x = newgri(i,1);y = newgri(i,2);
        allcount = allcount + 1.;
        i2 = i2+1;

        % calculate distance from center point and sort wrt distance
        l = sqrt(((xsecx' - x)).^2 + ((xsecy + y)).^2) ;
        %[s,is] = sort(l);
        %b = newa(is(:,1),:) ;       % re-orders matrix to agree row-wise

        % take first ni points
        l = l <= ra;
        b = newa.subset(l);      % new data per grid point (b) is sorted in distance

        if isempty(b) == 1; b = newa.subset(1); end
        if length(b(:,1)) >= 50;
            % call the b-value function
            lt =  b(:,3) >= t1 &  b(:,3) <t2 ;
            if  length(b(lt,1)) > 20;
                [bv magco stan av me mer me2,  pr] =  bvalca3(b(lt,:),inb1,inb2);
                bo1 = bv; no1 = newa.Count;
            else
                bv = NaN; pr = 50;
            end
            lt = b(:,3) >= t3 &  b(:,3) < t4 ;
            if  length(b(lt,1)) > 20;
                [bv2 magco stan av me mer me2,  pr] =  bvalca3(b(lt,:),inb1,inb2);
            else
                bv2 = NaN; pr = 50;
            end
            lt = b(:,3) >= t4 &  b(:,3) < t5 ;
            if  length(b(lt,1)) > 20;
                [bv3 magco stan av me mer me2,  pr2] =  bvalca3(b(lt,:),inb1,inb2);
            else
                bv3 = NaN; pr = 50;
            end

            l2 = sort(l);

            l2 = sort(l);
            b2 = b;
            if inb2 ==  1
                l = b(:,6) >= magco;
                b2 = b(l,:);
            end
            [av2 bv2 stan2 ] =  bmemag(b2);
            if pr >=90
                bvg = [bvg ; bv magco x y length(b(:,1)) bv2 pr av stan  max(b(:,6)) bv-bv2 bv-bv3 bv2-bv3 pr];
            else
                bvg = [bvg ; NaN NaN x y NaN NaN NaN NaN NaN  NaN NaN NaN NaN NaN];
            end
        else
            bvg = [bvg ; NaN NaN x y NaN NaN NaN NaN NaN  NaN NaN NaN NaN NaN];
        end
        waitbar(allcount/itotal)
    end  % for  newgri

    % save data
    %
    %  set(txt1,'String', 'Saving data...')
    drawnow
    gx = xvect;gy = yvect;

    catSave3 =...
        [ 'zmap_message_center.set_info(''Save Grid'',''  '');think;',...
        '[file1,path1] = uiputfile(fullfile(hodi, ''eq_data'', ''*.mat''), ''Grid Datafile Name?'') ;',...
        ' sapa2 = [''save '' path1 file1 '' ll tmpgri bvg xvect yvect gx gy ni dx dd par1 ni newa maex maey maix maiy ''];',...
        ' if length(file1) > 1, eval(sapa2),end , done']; eval(catSave3)

    close(wai)
    watchoff

    % reshape a few matrices
    %
    normlap2=ones(length(tmpgri(:,1)),1)*nan;
    normlap2(ll)= bvg(:,1);
    re3=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,5);
    r=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,6);
    meg=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,2);
    old1 =reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,7);
    pro=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,8);
    avm=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,9);
    stanm=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,10);
    maxm=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,11);
    db12=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,12);
    db13=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,13);
    db23=reshape(normlap2,length(yvect),length(xvect));

    normlap2(ll)= bvg(:,14);
    dp=reshape(normlap2,length(yvect),length(xvect));

    old = re3;

    % View the b-value map
    view_bvt

end   %  if sel = ca

% Load exist b-grid
if sel == 'lo'
    load_existing_bgrid_version_A
end

author: CelsoReyes


text:
filename = 'test_rectangular';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'circleInclusion'; %% Parámetro a ir variando
cost = {'enforceCh_CCstar_L2';'perimeter'};
%cost = {'enforceCh_CCstar_L2';'perimeter'};
weights = [1 1];
%weights = 1;
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'MMA';
%optimizer = 'AlternatingPrimalDual';

incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';
fracRadius = 0.5;

nsteps = 1;
Vfrac_final = 0.5;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 0.5;
optimality_initial = 1e-3;
constr_initial = 1e-3;
Perimeter_target = 5;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

 %Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
selectiveC_Cstar = 'IsotropyHexagon';

% For all tests
plotting = false;
printing = false;
printing_physics = false;
monitoring = false;
maxiter = 30;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function fx = p56_fun ( n, x )

%*****************************************************************************80
%
%% P56_FUN evaluates the integrand for problem 56.
%
%  Interval:
%
%    -1 <= x <= 1
%
%  Integrand:
%
%    1 / ( x^6 + 0.9 )
%
%  Approximate Integral (20 digits):
%
%    1.9922524079504000171...
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Kahaner,
%    Comparison of Numerical Quadrature Formulas,
%    in Mathematical Software,
%    edited by John R Rice,
%    Academic Press, 1971.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  fx(1:n) = 1.0 ./ ( x(1:n).^6 + 0.9 );

  return
end


text:
function pred = ml_predictsvmperf(trials, model)
% Prediction function for the Support Vector Machine (SVMperf).
% Prediction = ml_predictsvmperf(Trials, Model)
%
% In:
%   Trials  : the data a matrix, as in ml_predict
%
%   Model   : predictive model as produced by ml_trainsvmlight
%
% Out:
%   Prediction  : discrete probability distribution, formatted as
%                 {'disc' [NxC] [Cx1]}, with element #2 being the per-class probability and
%                 element #3 the original target values per class
%                 thus, the expected target values are Prediction{2}*Prediction{3}
%
% See also:
%   ml_trainsvmperf
%
%                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                           2010-04-04

if isfield(model,'voted')
    pred = ml_predictvote(trials,model);
else
    trials = hlp_applyscaling(trials,model.sc_info);
    raw_pred = min(+1,max(-1,svmperfclassify(trials, ones(size(trials,1),1), model, '-v 0')));
    pred = {'disc', [(1-raw_pred)/2 1-(1-raw_pred)/2], model.classes};
end

repo: REST


text:
function remove = spm_cfg_eeg_remove_bad_trials
% configuration file for copying
%__________________________________________________________________________
% Copyright (C) 2009-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_cfg_eeg_remove_bad_trials.m 5377 2013-04-02 17:07:57Z vladimir $

%--------------------------------------------------------------------------
% D
%--------------------------------------------------------------------------
D        = cfg_files;
D.tag    = 'D';
D.name   = 'File Name';
D.filter = 'mat';
D.num    = [1 1];
D.help   = {'Select the M/EEG mat file.'};

%--------------------------------------------------------------------------
% prefix
%--------------------------------------------------------------------------
prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the output dataset. Default prefix is ''r''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.val     = {'r'};

%--------------------------------------------------------------------------
% remove
%--------------------------------------------------------------------------
remove          = cfg_exbranch;
remove.tag      = 'remove';
remove.name     = 'Remove bad trials';
remove.val      = {D, prefix};
remove.help     = {'Removes bad trials and re-orders trials to conform to condlist'}';
remove.prog     = @eeg_remove;
remove.vout     = @vout_eeg_remove;
remove.modality = {'EEG'};

%==========================================================================
function out = eeg_remove(job)
% construct the S struct
S           = job;
S.D         = S.D{1};
out.D       = spm_eeg_remove_bad_trials(S);
out.Dfname  = {fullfile(out.D)};

%==========================================================================
function dep = vout_eeg_remove(job)
% return dependencies
dep(1)            = cfg_dep;
dep(1).sname      = 'Cleaned M/EEG data';
dep(1).src_output = substruct('.','D');
dep(1).tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2)            = cfg_dep;
dep(2).sname      = 'Cleaned M/EEG datafile';
dep(2).src_output = substruct('.','Dfname');
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});


text:
% Example of using CaptureFigVid
% Cheers, Dr. Alan Jennings, Research assistant professor,
% Department of Aeronautics and Astronautics, Air Force Institute of Technology

%% Set up 3D plot to record
figure(171);clf;
surf(peaks,'EdgeColor','none','FaceColor','interp','FaceLighting','phong')
daspect([1,1,.3]);axis tight;

%% Set up recording parameters (optional), and record
OptionZ.FrameRate=15;OptionZ.Duration=5.5;OptionZ.Periodic=true;
CaptureFigVid([-20,10;-110,10;-190,80;-290,10;-380,10], 'WellMadeVid',OptionZ)


function newstr = m2tstrjoin(cellstr, delimiter, floatFormat)
% This function joins a cell of strings to a single string (with a
% given delimiter in between two strings, if desired).
    if ~exist('delimiter','var') || isempty(delimiter)
        delimiter = '';
    end
    if ~exist('floatFormat','var') || isempty(floatFormat)
        floatFormat = '%g';
    end
    if isempty(cellstr)
        newstr = '';
        return
    end

    % convert all values to strings first
    nElem = numel(cellstr);
    for k = 1:nElem
        if isnumeric(cellstr{k})
            cellstr{k} = sprintf(floatFormat, cellstr{k});
        elseif iscell(cellstr{k})
            cellstr{k} = m2tstrjoin(cellstr{k}, delimiter, floatFormat);
            % this will fail for heavily nested cells
        elseif ~ischar(cellstr{k})
            error('matlab2tikz:join:NotCellstrOrNumeric',...
                'Expected cellstr or numeric.');
        end
    end

    % inspired by strjoin of recent versions of MATLAB
    newstr = cell(2,nElem);
    newstr(1,:)         = reshape(cellstr, 1, nElem);
    newstr(2,1:nElem-1) = {delimiter}; % put delimiters in-between the elements
    newstr(2,end)       = {''}; % for Octave 4 compatibility
    newstr = [newstr{:}];
end


text:
function newp = eleminterpolate(p,tree,weight)
%% ELEMINTERPOLATE interpolate a piecewise constant function.
%
% newp = eleminterpolate(p,tree) interpolate a piecewise constant function
% p from a coarse grid to a fine grid or the other way around. tree(:,1:3)
% stores the binary tree of the coarsening. tree(:,1) is the index of
% parent element in the coarsened mesh and tree(:,2:3) are two children
% indices in the original mesh.
%
% For uniform refinement, the prolongation from the coarse grid to the fine
% grid is simply |newp = repmat(p,4,1)| and from the fine grid to the
% coarse one is |newp = p(1:NT/4)|.
%
% Example
%
%   [node,elem] = squaremesh([0 1 0 1],1/2);
%   p = 1:size(elem,1);
%   figure(1);
%   subplot(1,4,1); showsolution(node,elem,p,'EdgeColor','k'); view(2);
%   [node,elem,~,~,tree] = bisect(node,elem,[1 2]);
%   p = eleminterpolate(p,tree);
%   subplot(1,4,2); showsolution(node,elem,p,'EdgeColor','k'); view(2);
%   [node,elem,~,~,tree] = bisect(node,elem,[1 2]);
%   p = eleminterpolate(p,tree);
%   subplot(1,4,3); showsolution(node,elem,p,'EdgeColor','k'); view(2);
%   [node,elem,~,~,tree] = coarsen(node,elem,'all');
%   p = eleminterpolate(p,tree);
%   subplot(1,4,4); showsolution(node,elem,p,'EdgeColor','k'); view(2);
%
% See also bisect, coarsen, nodeinterpolate
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('weight','var'), weight = 1; end
%%
newp = p;
if (size(tree,1)==0), return; end % no change
NTin = length(p);
NTf = max(tree(:,3));
if NTin < NTf 		   % coarse grid to fine grid
    NTc = NTin;
    newp = zeros(NTf,1);
    isNew = tree(:,3); % Right child is a new element
    inCoarse = true(NTf,1);
    inCoarse(isNew) = false;
    % Case 0: still in coarse element
    newp(inCoarse) = p;
    % Case 1: The parent element is in the coarse mesh. A triangle t could
    % be bisected twice and a child could be a parent
    idx = (tree(:,1) <= NTc);
    if weight == 1
        newp(tree(idx,2)) = p(tree(idx,1));
        newp(tree(idx,3)) = p(tree(idx,1));
    else
        newp(tree(idx,2)) = weight*p(tree(idx,1));
        newp(tree(idx,3)) = weight*p(tree(idx,1));
    end
    % Case 2: The parent is in the intermediate mesh
    idx = (tree(:,1) > NTc);
    if weight == 1
        newp(tree(idx,2)) = newp(tree(idx,1));
        newp(tree(idx,3)) = newp(tree(idx,1));
    else
        newp(tree(idx,2)) = weight*newp(tree(idx,1));
        newp(tree(idx,3)) = weight*newp(tree(idx,1));
    end
else      % fine grid to coarse grid
    newp = p;
    newp(tree(:,3)) = [];
end

repo: ifem


text:
function r8vec_uniform_01_test ( )

%*****************************************************************************80
%
%% R8VEC_UNIFORM_01_TEST tests R8VEC_UNIFORM_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_UNIFORM_01\n' );
  fprintf ( 1, '  R8VEC_UNIFORM_01 returns a random R8VEC\n' );
  fprintf ( 1, '  with entries in [ 0.0, 1.0 ]\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  for j = 1 : 3

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Input SEED = %12d\n', seed );
    fprintf ( 1, '\n' );

    [ r, seed ] = r8vec_uniform_01 ( n, seed );

    r8vec_print ( n, r, '  Random R8VEC:' );

  end

  return
end


text:
function coord_w = minc_voxel2world(coord_v,mat,opt);
% Convert coordinates in the voxel space into coordinates in the world
% space.
%
% SYNTAX:
% COORD_W = MINC_VOXEL2WORLD(COORD_V,MAT,OPT)
%
% _________________________________________________________________________
% INPUTS:
%
% COORD_V
%       (matrix N*3) each row is a vector of 3D coordinates in voxel space.
%
% MAT
%       (matrix 4*4) an affine transformation from voxel to world
%       coordinates. See the help of NIAK_READ_VOL for more infos. It is
%       generally the HDR.INFO.MAT field of the header of a volume file.
%
% OPT
%       (structure, optional) with the following fields :
%
%       FLAG_ZERO
%           (boolean, default false) if FLAG_ZERO is true, voxel
%           coordinates start from 1 (default behaviour in matlab),
%           otherwise they start from 0 (default behaviour in C/C++ or
%           MINC).
%
% _________________________________________________________________________
% OUTPUTS:
%
% COORD_W
%       (matrix N*3) each row is a vector of 3D coordinates in world space.
%
% _________________________________________________________________________
% SEE ALSO:
% MINC_READ, MINC_WRITE, MINC_VOXEL2WORLD, MINC_WORLD2VOXEL
%
% _________________________________________________________________________
% COMMENTS:
%
% Copyright (c) Pierre Bellec, Centre de recherche de l'institut de
% gériatrie de Montréal, Département d'informatique et de recherche
% opérationnelle, Université de Montréal, 2013.
% See licensing information in the code.
% Keywords : affine transformation, coordinates

% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
if nargin < 3
    flag_zero = false;
else
    if isfield(opt,'flag_zero')
        flag_zero = opt.flag_zero;
    else
        flag_zero = false;
    end
end
if flag_zero
    coord_w = [coord_v ones([size(coord_v,1) 1])]*(mat');
else
    coord_w = [coord_v-1 ones([size(coord_v,1) 1])]*(mat');
end
coord_w = coord_w(:,1:3);

author: brainstorm-tools


text:
function [m2] = yd22m2(yd2)
% Convert area from square yards to square meters.
% Chad A. Greene 2012
m2 = yd2*0.83612736;


text:
filename='Bridge_hexahedra_coarse';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;


text:
% Data File DTXY1
% Dynamics of a particle, thrown
% with initial velocity v0 under
% angle "alfa' to the Horizon
 m    =  'm';    % mass of the particle
 Fx   = '-k*xt'; % Projections of forces
 Fy   = '-k*yt - m*9.81'; % on axis x and y
 x0   = '0';     % Initial
 y0   = '0';     % coordinates
 v0   = 'v0';    % Initial velocity
 alfa = 'alfa';  % angle between v0 and horizon
 Tend = 20;      % upper bound of integration
 eps  = 1e-10;   % desirable accuracy
 np   = 2;       % number of parameters
 P{1} = 'm';     % mass of the particle
 P{2} = 'k';     % coefiicient of resistance


text:
function plotEmptySphere(varargin)
% plots white sphere

if check_option(varargin,'parent')
  ax = get_option(varargin,'parent');
else
  ax = gca;
end

[x, y, z] = sphere;
x = 0.99*x;
y = 0.99*y;
z = 0.99*z;
surface(x,y,z,'FaceColor', 'w','EdgeColor','none','parent',ax,'handlevisibility','off')
hold(ax,'on')

dth = 15*degree;
th = -pi/2+dth:dth:pi/2-dth;
rh = linspace(0,2*pi,100);

[th,rh] = meshgrid(th,rh);

[x,y,z] = sph2cart(rh, th, 1);

line(x,y,z,'color',[1 1 1] * 0.8,'parent',ax,'handlevisibility','off');

drh = 15*degree;
rh = 0:drh:2*pi-drh;
th = linspace(-pi/2,pi/2,50);

[th,rh] = meshgrid(th,rh);

[x,y,z] = sph2cart(rh, th, 1);

line(x.',y.',z.','color',[1 1 1] * 0.8,'parent',ax,'handlevisibility','off')

axis(ax,'equal','vis3d','off');
set(ax,'XDir','rev','YDir','rev',...
  'XLim',[-1.02,1.02],'YLim',[-1.02,1.02],'ZLim',[-1.02,1.02]);
view(3);

if nargout == 0, clear h;end

end

author: mtex-toolbox

repo: mtex

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


text:
function C = myunion(A,B)
% MYUNION Union of two sets of positive integers (much faster than built-in union)
% C = myunion(A,B)

if isempty(A)
  ma = 0;
else
  ma = max(A);
end

if isempty(B)
  mb = 0;
else
  mb = max(B);
end

if ma==0 & mb==0
  C = [];
elseif ma==0 & mb>0
  C = B;
elseif ma>0 & mb==0
  C = A;
else
  %bits = sparse(1, max(ma,mb));
  bits = zeros(1, max(ma,mb));
  bits(A) = 1;
  bits(B) = 1;
  C = find(bits);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [split,numpieces]=explode(string,delimiters)
% EXPLODE    Splits string into pieces.
%
% EXPLODE(STRING,DELIMITERS) returns a cell array with the pieces
% of STRING found between any of the characters in DELIMITERS.
%
% :Usage:
% ::
%
%     [SPLIT,NUMPIECES] = EXPLODE(STRING,DELIMITERS)
%
% also returns the number of pieces found in STRING.
%
% :Inputs:
%
%   **STRING:**
%        the string to split (string)
%
%   **DELIMITERS:**
%        the delimiter characters (string)
%
% :Outputs:
%
%   **SPLIT:**
%        the split string (cell array), each cell is a piece
%
%   **NUMPIECES:**
%        the number of pieces found (integer)
%
% :Examples:
% ::
%
%    STRING = 'ab_c,d,e fgh'
%    DELIMITERS = '_,'
%    [SPLIT,NUMPIECES] = EXPLODE(STRING,DELIMITERS)
%    SPLIT = 'ab'    'c'    'd'    'e fgh'
%    NUMPIECES = 4
%
% :See also:
% IMPLODE, STRTOK
%
% ..
%   Created: Sara Silva (sara@itqb.unl.pt) - 2002.04.30
%
%   Modified: Matthew Davidson (matthew@psych.columbia.edu) - 2005.10.11
% ..

if isempty(string) % empty string, return empty and 0 pieces
    split{1}='';
    numpieces=0;

elseif isempty(delimiters) % no delimiters, return whole string in 1 piece
    split{1}=string;
    numpieces=1;

else % non-empty string and delimiters, the correct case

    remainder=string;
    i=0;

    % If the first letter is in the delimiters, add a blank piece at the
    % beginning
    firstletteridx = strfind(delimiters,string(1));
    if ~isempty(firstletteridx)
        i=i+1;
        split{i}='';
        remainder = string(2:end);
    end

    while ~isempty(remainder)
        [piece,remainder]=strtok(remainder,delimiters);
        i=i+1;
        split{i}=piece;
    end
    numpieces=i;

end


text:
function flag=deletemeshfile(fname)
%
% flag=deletemeshfile(fname)
%
% delete a given work mesh file under the working directory
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
%
% input:
%     fname: specified file name (without path)
%
% output:
%     flag: not used
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

try
    if(exist(fname))
	delete(fname);
    end
catch
    error(['You do not have permission to delete temporary files. If you are working in a multi-user ',...
         'environment, such as Unix/Linux and there are other users using iso2mesh, ',...
         'you may need to define ISO2MESH_SESSION=''yourstring'' to make your output ',...
         'files different from others; if you do not have permission to ',mwpath(''),...
         ' as the temporary directory, you have to define ISO2MESH_TEMP=''/path/you/have/write/permission'' ',...
         'in matlab/octave base workspace.']);
end


text:
function json=jsonset(fname,mmap,varargin)
%
% json=jsonset(fname,mmap,'$.jsonpath1',newval1,'$.jsonpath2','newval2',...)
%
% Fast writing of JSON data records to stream or disk using memory-map
% (mmap) returned by loadjson/loadbj and JSONPath-like keys
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
% initially created on 2022/02/02
%
% input:
%      fname: a JSON/BJData/UBJSON string or stream, or a file name
%      mmap: memory-map returned by loadjson/loadbj of the same data
%            important: mmap must be produced from the same file/string,
%            otherwise calling this function may cause data corruption
%      '$.jsonpath1,2,3,...':  a series of strings in the form of JSONPath
%            as the key to each of the record to be written
%
% output:
%      json: the modified JSON string or, in the case fname is a filename,
%            the cell string made of jsonpaths that are successfully
%            written
%
% examples:
%      % create test data
%       d.arr={[1,2],'a',struct('c',2)}; d.obj=struct('k','test')
%      % convert to json string
%       str=savejson('',d,'compact',1)
%      % parse and return mmap
%       [dat, mmap]=loadjson(str);
%      % display mmap entries
%       savejson('',mmap)
%      % replace value using mmap
%       json=jsonset(str,mmap,'$.arr.[2].c','5')
%      % save same json string to file (must set savebinary 1)
%       savejson('',d,'filename','file.json','compact',1,'savebinary',1);
%      % fast write to file
%       json=jsonset('file.json',mmap,'$.arr.[2].c','5')
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(regexp(fname,'^\s*(?:\[.*\])|(?:\{.*\})\s*$','once'))
    inputstr=fname;
else
    if(~exist('memmapfile','file'))
        fid=fopen(fname,'r+b');
    end
end

mmap=[mmap{:}];
keylist=mmap(1:2:end);

opt=struct;
for i=1:2:length(varargin)
    if(isempty(regexp(varargin{i},'^\$', 'once')))
        opt.(encodevarname(varargin{i}))=varargin{i+1};
    end
end

json={};
for i=1:2:length(varargin)
    if(regexp(varargin{i},'^\$'))
        [tf,loc]=ismember(varargin{i},keylist);
        if(tf)
            bmap=mmap{loc*2};
            if(ischar(varargin{i+1}))
                val=varargin{i+1};
            else
                val=savejson('',varargin{i+1},'compact',1);
            end
            if(length(val)<=bmap(2))
                val=[val repmat(' ',[1,bmap(2)-length(val)])];
                if(exist('inputstr','var'))
                    inputstr(bmap(1):bmap(1)+bmap(2)-1)=val;
                else
                    if(exist('memmapfile','file'))
                        rec={'uint8', [1 bmap(2)],  'x'};
                        fmap=memmapfile(fname,'writable',true,'offset', bmap(1)-1, 'format', rec, 'repeat',1);
                        fmap.Data.x=uint8(val);
                    else
                        fseek(fid,bmap(1)-1,'bof');
                        fwrite(fid,val);
                    end
                    json{end+1}={varargin{i},val};
                end
            end
        end
    end
end

if(exist('fid','var') && fid>=0)
    fclose(fid);
end

if(exist('inputstr','var'))
    json=inputstr;
end