text:
% eeg_interp_sph_sline_test - script to test eeg_interp_sph_spline
%

clear

p = eeg_toolbox_defaults;

p.volt.path = 'e:\matlab\eeg_toolbox\eeg_example_data\';
p.volt.file = 'eeg_124ch_simulatedEMSE.txt'; % (400 rows x 125 columns)
p.volt.type = 'ascii';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load electrode co-ordinates

p = elec_open(p);
p.elec.plot = 1;

x = p.elec.data.x;
y = p.elec.data.y;
z = p.elec.data.z;

% % create spherical electrode positions
% r = 10
% fprintf('...generating spherical interpolation points\n');
% [x,y,z] = elec_sphere_points(16,24,r);
%
% p.elec.data.x = x;
% p.elec.data.y = y;
% p.elec.data.z = z;
%
% p.elec.data.Xsp = x;
% p.elec.data.Ysp = y;
% p.elec.data.Zsp = z;
% p.elec.data.Rsp = [r r r];
%
% p.elec.n = 337;
%elec_plot(p)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create simulated potential data

p.volt.sampleHz = 1000;
p.volt.sampleMsec = 1;
p.volt.sampleTime = 1 * p.volt.sampleMsec;
p.volt.samplePoint = 1;
p.volt.epochStart = 0;
p.volt.epochEnd = 1;

p.volt.var = [];
p.volt.timeArray = 0:1:10;
p.volt.points = 10;
p.volt.channels = length(x);
p.volt.epochStart = 0;
p.volt.epochEnd = 10;
p.volt.sweeps = 1;
p.volt.peaks = [];

p.volt.data = repmat(p.elec.data.x',10,1);

V = p.volt.data(p.volt.samplePoint,:);

p.clickTimePoint = 0;
p = eeg_contours_engine(p);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spherical spline interpolation

FV = eeg_interp_sph_spline(V,[x y z]);

p.elec.n = length(FV.vertices);
p.elec.data.x = FV.vertices(:,1);
p.elec.data.y = FV.vertices(:,2);
p.elec.data.z = FV.vertices(:,3);
p.elec.data.Xsp = FV.vertices(:,1);
p.elec.data.Ysp = FV.vertices(:,2);
p.elec.data.Zsp = FV.vertices(:,3);

p.volt.data = repmat(FV.Cdata,10,1);

p.volt.file = 'spherical spline interpolation';

p = eeg_contours_engine(p);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% scd interpolation

FV = eeg_interp_sph_spline_scd(V,[x y z]);

p.elec.n = length(FV.vertices);
p.elec.data.x = FV.vertices(:,1);
p.elec.data.y = FV.vertices(:,2);
p.elec.data.z = FV.vertices(:,3);
p.elec.data.Xsp = FV.vertices(:,1);
p.elec.data.Ysp = FV.vertices(:,2);
p.elec.data.Zsp = FV.vertices(:,3);

p.volt.data = repmat(FV.Cdata,10,1);

p.volt.file = 'spherical spline scd interpolation';

p = eeg_contours_engine(p);

author: PatternRecognition

repo: OpenBMI


text:
function [training_data, mappingScore] = prepareTrainingData(model,printLevel,params)
% given a standard COBRA model, add thermodynamic data to it using
% the Component Contribution method
%
% INPUTS
% model
%
%
% OPTIONAL INPUTS
% params.use_cached_kegg_inchis
% params.use_model_pKas_by_default
% params.uf                          maximum uncertainty
%
%
% OUTPUTS:
% model
% .DfG0                 m x 1 array of component contribution estimated
%                       standard Gibbs energies of formation.
% .covf                 m x m estimated covariance matrix for standard
%                       Gibbs energies of formation.
% .uf                   m x 1 array of uncertainty in estimated standard
%                       Gibbs energies of formation. uf will be large for
%                       metabolites that are not covered by component
%                       contributions.

if ~exist('printLevel','var')
    printLevel = 0;
end
if ~exist('param','var')
    use_cached_kegg_inchis=true;
    use_model_pKas_by_default=true;
else
    if ~isfield(params,'use_cached_kegg_inchis')
        use_cached_kegg_inchis = true;
        % use_cached_kegg_inchis = false;
    else
        use_cached_kegg_inchis=params.use_cached_kegg_inchis;
    end
    if ~isfield(params,'use_model_pKas_by_default')
        use_model_pKas_by_default = true;
    else
        use_model_pKas_by_default=params.use_model_pKas_by_default;
    end
end
% load the training data (from TECRDB, Alberty, etc.)
training_data = loadTrainingData();

% get the InChIs for all the compounds in the training data
% (note that all of them have KEGG IDs)
kegg_inchies = getInchies(training_data.cids, use_cached_kegg_inchis);
inds = ismember(kegg_inchies.cids, training_data.cids);
training_data.std_inchi = kegg_inchies.std_inchi(inds);
training_data.std_inchi_stereo = kegg_inchies.std_inchi_stereo(inds);
training_data.std_inchi_stereo_charge = kegg_inchies.std_inchi_stereo_charge(inds);
training_data.nstd_inchi = kegg_inchies.nstd_inchi(inds);

% use the chemical formulas from the InChIs to verify that each and every
% reaction is balanced.
training_data = balanceReactionsInTrainingData(training_data);

% get the pKas for the compounds in the training data (using ChemAxon)
training_data.kegg_pKa = getTrainingDatapKas(training_data);

% match between the compounds in the model and the KEGG IDs used in the
% training data, and create the group incidence matrix (G) for the
% combined set of all compounds.
[training_data, mappingScore] = createGroupIncidenceMatrix(model, training_data);

% apply the reverse Legendre transform for the relevant training observations (typically
% apparent reaction Keq from TECRDB)
training_data = reverseTransformTrainingData(model, training_data, use_model_pKas_by_default);

author: opencobra

repo: cobratoolbox


text:
function CPD = hhmmF_CPD(bnet, self, Qnodes, d, D, varargin)
% HHMMF_CPD Make the CPD for an F node at depth D of a D-level hierarchical HMM
% CPD = hhmmF_CPD(bnet, self, Qnodes, d, D, ...)
%
%    Q(d-1)
%          \
%           \
%           F(d)
%         /   |
%        /    |
%    Q(d)  F(d+1)
%
% We assume nodes are ordered (numbered) as follows:
% Q(1), ... Q(d), F(d+1), F(d)
%
% F(d)=2 means level d has finished. The prob this happens depends on Q(d)
% and optionally on Q(d-1), Q(d=1), ..., Q(1).
% Also, level d can only finish if the level below has finished
% (hence the F(d+1) -> F(d) arc).
%
% If d=D, there is no F(d+1), so F(d) is just a regular tabular_CPD.
% If all models always finish in the same state (e.g., their last),
% we don't need to condition on the state of parent models (Q(d-1), ...)
%
% optional args [defaults]
%
% termprob - termprob(k,i,2) = prob finishing given Q(d)=i and Q(1:d-1)=k [ finish in last state ]
%
% hhmmF_CPD is a subclass of tabular_CPD so we inherit inference methods like CPD_to_pot, etc.
%
% We create an isolated tabular_CPD with no F parent to learn termprob
% so we can avail of e.g., entropic or Dirichlet priors.
%
% For details, see "Linear-time inference in hierarchical HMMs", Murphy and Paskin, NIPS'01.

ps = parents(bnet.dag, self);
Qps = myintersect(ps, Qnodes);
F = mysetdiff(ps, Qps);
CPD.Q = Qps(end); % Q(d)
assert(CPD.Q == Qnodes(d));
CPD.Qps = Qps(1:end-1); % all Q parents except Q(d), i.e., calling context

ns = bnet.node_sizes(:);
CPD.Qsizes = ns(Qnodes);
CPD.d = d;
CPD.D = D;

Qsz = ns(CPD.Q);
Qpsz = prod(ns(CPD.Qps));

% set default arguments
p = 0.9;
%termprob(k,i,t) Might terminate if i=Qsz; will not terminate if i<Qsz
termprob = zeros(Qpsz, Qsz, 2);
termprob(:, Qsz, 2) = p;
termprob(:, Qsz, 1) = 1-p;
termprob(:, 1:(Qsz-1), 1) = 1;

for i=1:2:length(varargin)
  switch varargin{i},
   case 'termprob', termprob = varargin{i+1};
   otherwise, error(['unrecognized argument ' varargin{i}])
  end
end

ps = [CPD.Qps CPD.Q];
% ns(self) = 2 since this is an F node
CPD.sub_CPD_term = mk_isolated_tabular_CPD(ps, ns([ps self]), {'CPT', termprob});
S = struct(CPD.sub_CPD_term);
CPD.termprob = S.CPT;

CPD = class(CPD, 'hhmmF_CPD', tabular_CPD(bnet, self));

CPD = update_CPT(CPD);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% Author: Guosheng Lin (guosheng.lin@gmail.com)

% this is a simpler demo file for testing on your own images.

function demo_test_simple_voc()

rng('shuffle');
addpath('./my_utils');
dir_matConvNet='../libs/matconvnet/matlab';
run(fullfile(dir_matConvNet, 'vl_setupnn.m'));

run_config=[];

run_config.use_gpu=true;
% run_config.use_gpu=false;
run_config.gpu_idx=1;

% result dir:
result_name=['runner_result_dir' datestr(now, 'YYYYmmDDHHMMSS')];
result_dir=fullfile('../cache_data', 'test_examples_city', result_name);
mkdir_notexist(result_dir);

% the folder that contains testing images:
img_data_dir='../datasets/example_imgs_cityscapes';

% using a trained model which is trained on VOC 2012
run_config.trained_model_path='../model_trained/refinenet_res101_cityscapes.mat';
run_config.class_info=gen_class_info_cityscapes();

% for trained model, control the size of input images
run_config.input_img_short_edge_min=600;
run_config.input_img_short_edge_max=1100;

runner_info=prepare_runner_test_simple(run_config);

img_filenames=my_list_file(img_data_dir);
img_num=length(img_filenames);
for img_idx=1:img_num
    task_info=[];
    task_info.img_dir=img_data_dir;
    task_info.img_filename=img_filenames{img_idx};
    task_result=runner_info.run_task_fn(runner_info, task_info);

    [~, img_name]=fileparts(task_info.img_filename);
    one_cache_file=fullfile(result_dir, [img_name '.png']);
    fprintf('save prediction mask:%s\n', one_cache_file);
    imwrite(task_result.mask_data, run_config.class_info.mask_cmap, one_cache_file);
end

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
%  px = ellipseinteriorpixels(im,ell)
function px = ellipseinteriorpixels(im,ell)

[nr,nc,ncolors] = size(im);

% get a bounding box around the ellipse
bb = [ell.y+ell.a*[-1,1],ell.x+ell.a*[-1,1]];
bb([1,3]) = floor(bb([1,3]));
bb([2,4]) = ceil(bb([2,4]));
bb = max(bb,1);
bb(1:2) = min(bb(1:2),nr);
bb(3:4) = min(bb(1:2),nc);

% crop out this region
imbb = im(bb(1):bb(2),bb(3):bb(4),:);
bw = ellipsepixels([ell.x,ell.y,ell.a*2,ell.b*2,ell.theta],bb);
idx = find(bw);
n = length(idx);
px = zeros(n,ncolors);
for color = 1:ncolors,
  px(:,color) = imbb(idx+(color-1)*nr*nc);
end

author: kristinbranson

repo: JAABA


text:
function val=getoptkey(key,default,varargin)
%
% val=getoptkey(key,default,opt)
%    or
% val=getoptkey(key,default,'key1',val1,'key2',val2, ...)
%
% query the value of a key from a structure or a list of key/value pairs
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%   key: a string name for the target struct field name
%   default: the default value of the key is not found
%   opt: a struct object; the field names will be searched to match the
%        key input, opt can be a list of 'keyname'/value pairs
%
% output:
%   val: val=opt.key if found, otherwise val=default
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

val=default;
if(nargin<=2) return; end
opt=varargin2struct(varargin{:});
if(isstruct(opt) && isfield(opt,key))
    val=getfield(opt,key);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef LaunchVehicleNonSeqEvents < matlab.mixin.SetGet & matlab.mixin.Copyable
    %LaunchVehicleNoneSeqEvents Summary of this class goes here
    %   Detailed explanation goes here

    properties
        nonSeqEvts LaunchVehicleNonSeqEvent

        lvdData LvdData
    end

    properties(Dependent)
        evts LaunchVehicleEvent
    end

    methods
        function obj = LaunchVehicleNonSeqEvents(lvdData)
            obj.lvdData = lvdData;
            obj.nonSeqEvts = LaunchVehicleNonSeqEvent.empty(1,0);
        end

        function values = get.evts(obj)
            values = LaunchVehicleEvent.empty(1,0);

            for(i=1:length(obj.nonSeqEvts))
                values(end+1) = obj.nonSeqEvts(i).evt; %#ok<AGROW>
            end
        end

        function activeNonSeqEvts = getNonSeqEventsForScriptEvent(obj, evt)
            activeNonSeqEvts = LaunchVehicleNonSeqEvent.empty(1,0);
            evtNum = evt.getEventNum();

            for(i=1:length(obj.nonSeqEvts))
                nonSeqEvt = obj.nonSeqEvts(i);

                if(nonSeqEvt.numExecsRemaining > 0)
                    if(not(isempty(nonSeqEvt.lwrBndEvt)))
                        lwrBndEvtNum = nonSeqEvt.lwrBndEvt.getEventNum();
                    else
                        lwrBndEvtNum = 0;
                    end

                    if(not(isempty(nonSeqEvt.uprBndEvt)))
                        uprBndEvtNum = nonSeqEvt.uprBndEvt.getEventNum();
                    else
                        uprBndEvtNum = Inf;
                    end

                    if(evtNum >= lwrBndEvtNum && evtNum <= uprBndEvtNum)
                        activeNonSeqEvts(end+1) = nonSeqEvt; %#ok<AGROW>
                    end
                end
            end
        end

        function addEvent(obj, newEvt)
            obj.nonSeqEvts(end+1) = newEvt;
        end

        function addEventAtInd(obj, newEvt, ind)
            if(not(isempty(obj.nonSeqEvts)))
                if(ind == length(obj.nonSeqEvts))
                    obj.nonSeqEvts(end+1) = newEvt;
                else
                    obj.nonSeqEvts = [obj.nonSeqEvts(1:ind), newEvt, obj.nonSeqEvts(ind+1:end)];
                end
            else
                obj.nonSeqEvts(end+1) = newEvt;
            end
        end

        function removeEvent(obj, nonSeqEvt)
            arguments
                obj(1,1) LaunchVehicleNonSeqEvents
                nonSeqEvt(1,1) LaunchVehicleNonSeqEvent
            end

            evt = nonSeqEvt.evt;

            termCondOptVar = evt.termCond.getExistingOptVar();
            if(not(isempty(termCondOptVar)))
                obj.lvdData.optimizer.vars.removeVariable(termCondOptVar);
            end

            actions = evt.actions;
            for(i=1:length(actions)) %#ok<*NO4LP>
                evt.removeAction(actions(i));
            end

            obj.nonSeqEvts(obj.nonSeqEvts == nonSeqEvt) = [];
        end

        function removeEventFromIndex(obj, ind)
            if(ind >= 1 && ind <= length(obj.nonSeqEvts))
                obj.removeEvent(obj.nonSeqEvts(ind));
            end
        end

        function evtNum = getNumOfEvent(obj, evt)
            evtNum = [];

            if(not(isempty(evt)))
                evtNum = find(obj.nonSeqEvts == evt);
            end
        end

        function evt = getEventForInd(obj, ind)
            evt = LaunchVehicleNonSeqEvent.empty(1,0);

            if(ind >= 1 && ind <= length(obj.nonSeqEvts))
                evt = obj.nonSeqEvts(ind);
            end
        end

        function numEvents = getTotalNumOfEvents(obj)
            numEvents = length(obj.nonSeqEvts);
        end

        function [listboxStr, events] = getListboxStr(obj)
            listboxStr = cell(length(obj.nonSeqEvts),1);

            for(i=1:length(obj.nonSeqEvts))
                listboxStr{i} = obj.nonSeqEvts(i).getListboxStr();
            end

            events = obj.nonSeqEvts;
        end

        function resetAllNumExecsRemaining(obj)
            for(i=1:length(obj.nonSeqEvts))
                obj.nonSeqEvts(i).resetNumExecsRemaining();
            end
        end

        function tf = usesStage(obj, stage)
            tf = false;

            for(i=1:length(obj.evts))
                tf = tf || obj.evts(i).usesStage(stage);
            end
        end

        function tf = usesEngine(obj, engine)
            tf = false;

            for(i=1:length(obj.evts))
                tf = tf || obj.evts(i).usesEngine(engine);
            end
        end

        function tf = usesTank(obj, tank)
            tf = false;

            for(i=1:length(obj.evts))
                tf = tf || obj.evts(i).usesTank(tank);
            end
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;

            for(i=1:length(obj.evts))
                tf = tf || obj.evts(i).usesEngineToTankConn(engineToTank);
            end
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;

            for(i=1:length(obj.evts))
                tf = tf || obj.evts(i).usesStopwatch(stopwatch);
            end
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;

            for(i=1:length(obj.evts))
                tf = tf || obj.evts(i).usesExtremum(extremum);
            end
        end

        function tf = usesCalculusCalc(obj, calculusCalc)
            tf = false;

            for(i=1:length(obj.evts))
                tf = tf || obj.evts(i).usesExtremum(calculusCalc);
            end
        end

        function tf = usesPwrSink(obj, powerSink)
            tf = false;

            for(i=1:length(obj.evts))
                tf = tf || obj.evts(i).usesPwrSink(powerSink);
            end
        end

        function tf = usesPwrSrc(obj, powerSrc)
            tf = false;

            for(i=1:length(obj.evts))
                tf = tf || obj.evts(i).usesPwrSrc(powerSrc);
            end
        end

        function tf = usesPwrStorage(obj, powerStorage)
            tf = false;

            for(i=1:length(obj.evts))
                tf = tf || obj.evts(i).usesPwrStorage(powerStorage);
            end
        end

        function tf = usesSensor(obj, sensor)
            tf = false;

            for(i=1:length(obj.evts))
                tf = tf || obj.evts(i).usesSensor(sensor);
            end
        end

        function tf = usesTankToTankConn(obj, tankToTank)
            tf = false;

            for(i=1:length(obj.nonSeqEvts))
                tf = tf || obj.nonSeqEvts(i).usesTankToTankConn(tankToTank);
            end
        end
    end

	methods(Access = protected)
        function cpObj = copyElement(obj)
            cpObj = copyElement@matlab.mixin.Copyable(obj);

            for(i=1:length(obj.nonSeqEvts))
                cpObj.nonSeqEvts(i) = obj.nonSeqEvts(i).copy();
            end
        end
	end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function value = c1_geg_monomial_integral ( alpha, expon )

%*****************************************************************************80
%
%% C1_GEG_MONOMIAL_INTEGRAL: integral of monomial with Gegenbauer weight on C1.
%
%  Discussion:
%
%    C1_GEG is the interval [-1,+1] with the Gegenbauer weight function
%
%      w(alpha;x) = (1-x^2)^alpha
%
%    with -1.0 < alpha.
%
%    value = integral ( -1 <= x <= +1 ) x^expon (1-x^2)^alpha dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the exponent of (1-X^2).
%    - 1.0 < ALPHA.
%
%    Input, integer EXPON, the exponent.
%    0 <= EXPON.
%
%    Output, real VALUE, the value of the integral.
%
  if ( alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'C1_GEG_MONOMIAL_INTEGRAL - Fatal error!\n' );
    fprintf ( 1, '  ALPHA <= -1.0\n' );
    error ( 'C1_GEG_MONOMIAL_INTEGRAL - Fatal error!' );
  end

  if ( mod ( expon, 2 ) == 1 )
    value = 0.0;
    return
  end

  c = expon;

  arg1 = - alpha;
  arg2 =   1.0 + c;
  arg3 =   2.0 + alpha + c;
  arg4 = - 1.0;

  value1 = r8_hyper_2f1 ( arg1, arg2, arg3, arg4 );

  value = 2.0 * gamma ( 1.0 + c ) * gamma ( 1.0 + alpha ) ...
    * value1 / gamma ( 2.0 + alpha  + c );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ Y ] = classify1( Model, X )

    % Extract Features from the test data
    data = getImages(X);

    weights = Model.weights;
    projection = Model.projection;

    % Project the test data on the new dimension
    xTest = getFeatures(data, projection);

    % Classify the images
    [Y, ~] = LRClassify( weights, xTest );
end

function [Features] =  getImages(data)
    Features = [];
    for i = 1:size(data,1)
        image = reshape(data(i,:),[32,32,3]);
        image = imresize(image,4);
        feat = extract_feature(image);
        Features = horzcat(Features,feat);
    end
    Features = Features';
end

function [Feats] = getFeatures(Features, projection)
    featureSize = size(Features,2);

    % Perform Mean Normalization on the Feature Matrix.
    % Each pixel can have intensity values from 0 - 255
    for i = 1 : featureSize
        Features(:,i) = Features(:,i) - mean(Features(:,i));
    end

    % Project the Features on the new dimension
    Feats = projection * Features';
    Feats = Feats';
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [ n_data, a, b, fx ] = arithmetic_geometric_mean_values ( n_data )

%*****************************************************************************80
%
%% AGM_VALUES returns some values of the AGM.
%
%  Discussion:
%
%    The AGM is defined for nonnegative A and B.
%
%    The AGM of numbers A and B is defined by setting
%
%      A(0) = A,
%      B(0) = B
%
%      A(N+1) = ( A(N) + B(N) ) / 2
%      B(N+1) = sqrt ( A(N) * B(N) )
%
%    The two sequences both converge to AGM(A,B).
%
%    In Mathematica, the AGM can be evaluated by
%
%      ArithmeticGeometricMean [ a, b ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real A, B, the argument ofs the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 14;

  a_vec = [ ...
     22.0, ...
     83.0, ...
     42.0, ...
     26.0, ...
      4.0, ...
      6.0, ...
     40.0, ...
     80.0, ...
     90.0, ...
      9.0, ...
     53.0, ...
      1.0, ...
      1.0, ...
      1.0, ...
      1.5 ];
  b_vec = [ ...
     96.0, ...
     56.0, ...
      7.0, ...
     11.0, ...
     63.0, ...
     45.0, ...
     75.0, ...
      0.0, ...
     35.0, ...
      1.0, ...
     53.0, ...
      2.0, ...
      4.0, ...
      8.0, ...
      8.0 ];
  fx_vec = [ ...
     52.274641198704240049, ...
     68.836530059858524345, ...
     20.659301196734009322, ...
     17.696854873743648823, ...
     23.867049721753300163, ...
     20.717015982805991662, ...
     56.127842255616681863, ...
      0.000000000000000000, ...
     59.269565081229636528, ...
     3.9362355036495554780, ...
     53.000000000000000000, ...
     1.4567910310469068692, ...
     2.2430285802876025701, ...
     3.6157561775973627487, ...
     4.0816924080221632670 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0.0;
    b = 0.0;
    fx = 0.0;
  else
    a = a_vec(n_data);
    b = b_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function model = addInchiToModel(model, molFileDir,  method, printLevel)
% Assigns Inchi to model.inchi structure given set of mol files for each
% metabolite
%
% USAGE:
%
%    model = assignInchiToModel(model, molFileDir, printLevel)
%
% INPUT:
%    model          Model structure with following fields:
%
%                     * .S - `m x n` stoichiometric matrix.
%                     * .mets - `m x 1` array of metabolite identifiers.
%                     * .rxns - `n x 1` array of reaction identifiers.
%                     * .metFormulas - `m x 1` cell array of metabolite formulas. Formulas for
%                       protons should be H, and formulas for water should be H2O.
%                     * .metCharges - `m x 1` numerical array of metabolite charges.
%                     * .metCompartments - optional `m x 1` array of metabolite compartment
%                       assignments. Not required if metabolite
%                       identifiers are strings of the format `ID[*]`
%                       where * is the appropriate compartment identifier.
%
%    molFileDir:    Path to a directory containing molfiles for the
%                   major tautomer of the major microspecies of
%                   each metabolite at pH 7. Molfiles should be
%                   named with the metabolite identifiers in
%                   model.mets (without compartment assignments).
%
% OPTIONAL INPUTS:
%    printLevel:    Verbose level
%
% OUTPUTS:
%    model:          Model structure with following additional fields:
%
%                   * .inchi - Structure containing four `m x 1` cell array's of
%                     IUPAC InChI strings for metabolites, with varying
%                     levels of structural detail.
%
%                   * .inchi.standard: m x 1 cell array of standard inchi
%                   * .inchi.standardWithStereo: m x 1 cell array of standard inchi with stereo
%                   * .inchi.standardWithStereoAndCharge: m x 1 cell array of standard inchi with stereo and charge
%                   * .inchi.nonstandard: m x 1 cell array of non-standard inchi
%
%   inchiBool           m x 1 true if inchi exists
%   molBool             m x 1 true if mol file exists
%   compositeInchiBool  m x 1 true if inchi is composite

% Written output - MetStructures.sdf - An SDF containing all structures input to the
% component contribution method for estimation of standard Gibbs energies.

if ~exist('printLevel','var')
    printLevel = 1;
end
if ~exist('method','var')
    method = 'sdf';
end

switch method
    case 'sdf'
        % Get metabolite structures in molfileDir.
        if printLevel>0
            fprintf('Creating MetStructures.sdf from molfiles.\n')
        end
        sdfFileName = 'MetStructures.sdf';
        includeRs = 0; % Do not include structures with R groups in SDF
        [sdfMetList,noMolMetList] = mol2sdf(model.mets,molFileDir,sdfFileName,includeRs);

        sdfBool = ismember(model.mets,sdfMetList);

        model.molBool = ~ismember(model.mets,noMolMetList);

        if printLevel>0
            fprintf('Converting SDF to InChI strings.\n')
        end
        model.inchi = createInChIStruct(model.mets,[molFileDir filesep sdfFileName]);
    case 'mol'
        %assumes that model.mets without compartment info provides the name of
        %the mol files in molFileDir
        [model.inchi,model.molBool] = createInChIStruct(model.mets,[],molFileDir);
    otherwise
        error('unrecognised method')
end

%identify the metabolites without inchi
model.inchiBool = ~cellfun('isempty',model.inchi.nonstandard);

% Remove InChI for composite compounds as they cause problems later.
model.compositeInchiBool = ~cellfun('isempty',regexp(model.inchi.nonstandard,'\.'));
model.inchi.standard(model.compositeInchiBool) = cell(sum(model.compositeInchiBool),1);
model.inchi.standardWithStereo(model.compositeInchiBool) = cell(sum(model.compositeInchiBool),1);
model.inchi.standardWithStereoAndCharge(model.compositeInchiBool) = cell(sum(model.compositeInchiBool),1);
model.inchi.nonstandard(model.compositeInchiBool) = cell(sum(model.compositeInchiBool),1);

if printLevel>0
    fprintf('%u%s\n',length(model.mets),' = number of model metabolites')
    fprintf('%u%s\n',nnz(model.molBool),' ... with mol files')
    fprintf('%u%s\n',nnz(~model.molBool),' ... without mol files')
    fprintf('%u%s\n',nnz(model.inchiBool),' ... with nonstandard inchi')
    fprintf('%u%s\n',nnz(~model.inchiBool),' ... without nonstandard inchi')
    fprintf('%u%s\n',nnz(model.compositeInchiBool),' ... compositie inchi removed')
end


text:
function options = checkspelling (options)

  potential_options = {
    'K';...
    'onpower';...
    'leida';...
    'order';...
    'maxorder';...
    'covtype';...
    'zeromean';...
    'embeddedlags';...
    'embeddedlags_batched';
    'pca';...
    'pca_spatial';...
    'lowrank';...
    'varimax';...
    'maxFOth'; ...
    'pcamar';...
    'pcapred';...
    'standardise';...
    'standardise_pc';...
    'crosstermsonly';...
    'filter';...
    'detrend';...
    'downsample';...
    'leakagecorr';...
    'sequential';...
    'A';...
    'As';...
    'B';...
    'V';...
    'ndim';...
    'vcomp';...
    'timelag';...
    'exptimelag';...
    'orderoffset';...
    'symmetricprior';...
    'grouping';...
    'priorcov_rate';...
    'AR';...
    'uniqueAR';...
    'id_mixture';...
    'cluster';...
    'Pstructure';...
    'Pistructure';...
    'plotGamma';...
    'tuda';...
    'tudamonitoring';...
    'behaviour';...
    'S';...
    'Sind';...
    'prior';...
    'initial_hmm';...
    'state';...
    'Fs';...
    'cyc';...
    'tol';...
    'stopcriterion';...
    'meancycstop';...
    'cycstogoafterevent';...
    'DirStats';...
    'initTestSmallerK';...
    'initcyc';...
    'initrep';...
    'inittype';...
    'initcriterion';...
    'init_windowsize';...
    'Gamma';...
    'hmm';...
    'fehist';...
    'DirichletDiag';...
    'DirichletDiag_Init';...
    'dropstates';...
    'repetitions';...
    'updateObs';...
    'updateGamma';...
    'updateP';...
    'decodeGamma';...
    'keepS_W';...
    'useParallel';...
    'useMEX';...
    'verbose';...
    'cvfolds';...
    'cvmode';...
    'cvverbose';...
    'BIGNbatch';... % STOCHASTIC
    'BIGNinitbatch';...
    'BIGNinit';...
    'BIGprior';...
    'BIGcyc';...
    'BIGmincyc';...
    'BIGundertol_tostop';...
    'BIGcycnobetter_tostop';...
    'BIGtol';...
    'BIGinitrep';...
    'BIGinitcyc';...
    'BIGdecodeGamma';...
    'BIGforgetrate';...
    'BIGdelay';...
    'BIGbase_weights';...
    'BIGcomputeGamma';...
    'BIGverbose';...
    'BIGinitStrategy';...
    'p';...
    'removezeros';...
    'completelags';...
    'rlowess';...
    'numIterations';...
    'tol';...
    'pad';... % SPECTRA
    'Fs';...
    'fpass';...
    'tapers';...
    'win';...
    'to_do';...
    'loadings';...
    'Nf';...
    'MLestimation';...
    'level';...
    'PriorWeightingP';...
    'PriorWeightingPi';...
    'distribution';...
    'logisticYdim';...
    'balancedata';...
    'regularisation';...
    'classifier';...
    'encodemodel';...
    'Gamma_constraint';...
    'Gamma_constraint_radius';...
    'acrosstrial_constrained';...
    'episodic';...
    'ehmm_priorOFFvsON';...
    'ehmm_regularisation_baseline';...
    'ehmm_init_from_hmm';...
    'ehmm_initfromhmm_crit';...
    'ehmm_baseline_w';...
    'ehmm_baseline_data';...
  };

  current_options = fieldnames(options);
  unrecognized = current_options(~ismember(current_options,potential_options));

  for j = 1:length(unrecognized)
    % Check edit distance
    distances = cellfun(@(x) EditDistance(unrecognized{j},x),potential_options);
    [d,idx] = min(distances);
    if d <= 2 % Be selective, it will only suggest something based on a minor typo
      fprintf(2,sprintf('%s is not a recognized option - did you mean ''%s''?\n',unrecognized{j},potential_options{idx}));
    else
      fprintf(2,sprintf('%s is not a recognized option\n',unrecognized{j}))
    end

    options = rmfield(options,unrecognized{j});
  end

end

function [V,v] = EditDistance(string1,string2)
  % Edit Distance is a standard Dynamic Programming problem. Given two strings s1 and s2,
  % the edit distance between s1 and s2 is the minimum number of operations required to convert string
  % FROM FILE EXCHANGE
  % by : Reza Ahmadzadeh (seyedreza_ahmadzadeh@yahoo.com - reza.ahmadzadeh@iit.it)
  % 14-11-2012

  m=length(string1);
  n=length(string2);
  v=zeros(m+1,n+1);
  for i=1:1:m
      v(i+1,1)=i;
  end
  for j=1:1:n
      v(1,j+1)=j;
  end
  for i=1:m
      for j=1:n
          if (string1(i) == string2(j))
              v(i+1,j+1)=v(i,j);
          else
              v(i+1,j+1)=1+min(min(v(i+1,j),v(i,j+1)),v(i,j));
          end
      end
  end
  V=v(m+1,n+1);
end

author: OHBA-analysis

repo: HMM-MAR


text:
function [tablePathwayRxns] = runFindSparsePathway(model,modelType,start,stop,via,CandidateSubSystem,specificMetabolites2Ignore);
% This function takes in a model (either metabolic model or a whole-body
% model), a start metabolite, and a stop metabolite and finds a shortest
% pathway
ignoreMetabolites = {'atp','h','coa','pi','h2o','o2','nh4','h2o2','nadph','nadh','nadp','nad','adp','udp','na1','hco3',...
    'ppi', 'imp','xmp','amp','ppa','fad','fadh2','oh1','co2','gmp','so4','pap','paps','mg2','zn2',...
    'cl','ca2','cu2','adpcbl','fe3','vitd3','i','caro','mn2','fe2','avite1','btn','ascb_L','phyQ','pnto_R','pydx','k'};

%% do not change hereafter

ignoreMetabolites = [specificMetabolites2Ignore';ignoreMetabolites'];
tablePathwayRxns ='';
s=' ';
osenseStr = 'max';

model = addReaction(model,['start_', start],'reactionFormula',[start,s,'<=>']);
model = addReaction(model,['stop_', stop],'reactionFormula',[stop,s,'<=>']);
printRxnFormula(model,model.rxns(end-1:end))

model2 = model;
compartments = {'e','c','m','n','p','r','l','g','x'};
ignoreMet = [];
for i = 1 : length(compartments)
    ignoreM =strcat(ignoreMetabolites,'[',compartments{i},']');
    ignoreMet = [ignoreMet;ignoreM'];
end
if strcmp(modelType,'WBM')
    sex = model.sex;
    % load list of organs
    OrganLists;
    % expand list of compartments for biofluids
    biofluids = {'bc','bcK','bp','d','u','lu','luLI','luSI','fe','csf','luC','a','aL','sw','luI','bpC','bpL'};

    for i = 1 : length(biofluids)
        ignoreM =strcat(ignoreMetabolites,'[',biofluids{i},']');
        ignoreMet = [ignoreMet;ignoreM'];
    end
    clear ignoreM
    ignoreMet2 =[ignoreMet];
    for i = 1 : length(OrgansListExt)
        ignoreM =strcat(OrgansListExt{i},'_',ignoreMet);
        ignoreMet2 = [ignoreMet2;ignoreM];
    end
    ignoreMet = ignoreMet2;
end

[model2, rxnRemoveList] = removeMetabolites(model2, ignoreMet,true,'legacy');

%Excl_R = [find(contains(model2.rxns,'DM_'));find(contains(model2.rxns,'sink_'));];
%model2 = changeRxnBounds(model2,model2.rxns(Excl_R),0,'b');

% needed for WBM
% set whole
if strcmp(modelType,'WBM')
    model2 = changeRxnBounds(model2,'Whole_body_objective_rxn',0,'l');
    model2 = changeRxnBounds(model2,'Whole_body_objective_rxn',1,'u');
end
%remove all constraints
model2.lb(find(model2.lb<0))=-10000000;
model2.ub(find(model2.ub<0))=0;
model2.lb(find(model2.lb>0))=0;
model2.ub(find(model2.ub>0))=10000000;

% remove coupling constraints if present
if isfield(model2,'C')
    model2 = rmfield(model2,'C');
    model2 = rmfield(model2,'ctrs');
end

% does not work for WBM

%model = changeRxnBounds(model,model.rxns(find(contains(model.rxns,'EX_'))),-1,'l');

model2 = changeObjective(model2,['stop_', stop]);
model2 = changeRxnBounds(model2,['start_', start],-1000,'b');
model2 = changeRxnBounds(model2,['stop_', stop],1000,'l');

if ~isempty(via)
    % enforce flux through via reactions

 %   model2 = changeRxnBounds(model2,via,-1*ones(length(via),1),'u');
    model2 = changeRxnBounds(model2,via,1*ones(length(via),1),'l');
end

FBA = optimizeCbModel(model2);

if FBA.origStat == 1 %problem is feasible
    %%
    param.printLevel = 0;
    param.regularizeOuter = 1;
    rxnPenalty = ones(length(model2.rxns),1);
    % not implemented yet
    %RS = findRxnsFromSubSystem(model2,CandidateSubSystem);
    %rxnPenalty(ismember(model2.rxns,RS)) = 0;
    rxnPenalty(ismember(model2.rxns,['stop_', stop])) = -1;
    rxnPenalty(ismember(model2.rxns,['start_', start])) = -1;
    rxnPenalty(ismember(model2.rxns,via)) = -1;
    [solution,sparseRxnBool2] = findSparsePathway(model2,rxnPenalty,param);
    sparse = model2.rxns(sparseRxnBool2)
    nnz(sparseRxnBool2)
    nnz(solution.v)
    pathwayRxnsAbbr = sparse;
    pathwayRxnsFormulae = printRxnFormula(model,'rxnAbbrList',sparse,'printFlag',0);
    tablePathwayRxns = table(pathwayRxnsAbbr,pathwayRxnsFormulae);
    nnz(solution.v)
    if 1
        % without ignore metabolites
        [involvedMets, deadEnds] = draw_by_rxn(model2, sparse, 'true');
    else
        % with ignored metabolites
        [involvedMets, deadEnds] = draw_by_rxn(model, sparse, 'true');
    end
else
    fprintf('Problem is infeasible.\n');
end


text:
function L = chol_up (A)
%CHOL_UP up-looking Cholesky factorization.
% Example:
%   L = chol_up (A)
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

n = size (A) ;
L = zeros (n) ;
for k = 1:n
    L (k,1:k-1) = (L (1:k-1,1:k-1) \ A (1:k-1,k))' ;
    L (k,k) = sqrt (A (k,k) - L (k,1:k-1) * L (k,1:k-1)') ;
end


text:
function imResult = blendMode_CreateResult(imResult, C, offsetW, offsetH)
%% Combine the result of the blending to obtain the final image. This
%   function is used by every blend function.
%
% Input:
%       A       -       Base Image
%       B       -       Top Image
%   offsetW     -   move picture B horizontally in respect to the top-left
%                   corner of picture A. Default value = 1.
%   offsetH     -   move picture B vertically in respect to the top-left
%                   corner of picture A. Default value = 1.
%
% Output:
%       imResult    -   Result of the blending, having the same size of the
%                       Base Image A.
%

sW = offsetW;
eW = offsetW + size(C, 2) - 1;
sH = offsetH;
eH = offsetH + size(C, 1) - 1;

if offsetW < 1
    sW = 1;
    eW = size(C, 2);
end

if offsetH < 1
    sH = 1;
    eH = size(C, 1);
end

if ((sW <= eW) && (sH <= eH))
    imResult(sH : eH, sW : eW, :) = C;
end

end

author: Sable

repo: mcbench-benchmarks


text:
function [count edges mid loc] = histcn(X, varargin)
% function [count edges mid loc] = histcn(X, edge1, edge2, ..., edgeN)
%
% Purpose: compute n-dimensional histogram
%
% INPUT
%   - X: is (M x N) array, represents M data points in R^N
%   - edgek: are the bin vectors on dimension k, k=1...N.
%     If it is a scalar (Nk), the bins will be the linear subdivision of
%     the data on the range [min(X(:,k)), max(X(:,k))] into Nk
%     sub-intervals
%     If it's empty, a default of 32 subdivions will be used
%
% OUTPUT
%   - count: n-dimensional array count of X on the bins, i.e.,
%         count(i1,i2,...,iN) = cardinal of X such that
%                  edge1(i1) <= X(:,i1) < edge1(i1)+1 and
%                       ...
%                  edgeN(iN) <= X(:,iN) < edgeN(iN)+1
%   - edges: (1 x N) cell, each provides the effective edges used in the
%     respective dimension
%   - mid: (1 x N) cell, provides the mid points of the cellpatch used in
%     the respective dimension
%   - loc: (M x N) array, index location of X in the bins. Points have out
%     of range coordinates will have zero at the corresponding dimension.
%
% DATA ACCUMULATE SYNTAX:
%   [ ... ] = histcn(..., 'AccumData', VAL);
%   where VAL is M x 1 array. Each VAL(k) corresponds to position X(k,:)
%   will be accumulated in the cell containing X. The accumulate result
%   is returned in COUNT.
%   NOTE: Calling without 'AccumData' is similar to having VAL = ones(M,1)
%
%   [ ... ] = histcn(..., 'AccumData', VAL, 'FUN', FUN);
%     applies the function FUN to each subset of elements of VAL.  FUN is
%     a function that accepts a column vector and returns
%     a numeric, logical, or char scalar, or a scalar cell.  A has the same class
%     as the values returned by FUN.  FUN is @SUM by default.  Specify FUN as []
%     for the default behavior.
%
% Usage examples:
%   M = 1e5;
%   N = 3;
%   X = randn(M,N);
%   [N edges mid loc] = histcn(X);
%   imagesc(mid{1:2},N(:,:,ceil(end/2)))
%
% % Compute the mean on rectangular patch from scattered data
%   DataSize = 1e5;
%   Lat = rand(1,DataSize)*180;
%   Lon = rand(1,DataSize)*360;
%   Data = randn(1,DataSize);
%   lat_edge = 0:1:180;
%   lon_edge = 0:1:360;
%   meanData = histcn([Lat(:) Lon(:)], lat_edge, lon_edge, 'AccumData', Data, 'Fun', @mean);
%
% See also: HIST, ACCUMARRAY
%
% Bruno Luong: <brunoluong@yahoo.com>
% Last update: 25/August/2011

if ndims(X)>2
    error('histcn: X requires to be an (M x N) array of M points in R^N');
end
DEFAULT_NBINS = 32;

AccumData = [];
Fun = {};

% Looks for where optional parameters start
% For now only 'AccumData' is valid
split = find(cellfun('isclass', varargin, 'char'), 1, 'first');
if ~isempty(split)
    for k = split:2:length(varargin)
        if strcmpi(varargin{k},'AccumData')
            AccumData = varargin{k+1}(:);
        elseif strcmpi(varargin{k},'Fun')
            Fun = varargin(k+1); % 1x1 cell
        end
    end
    varargin = varargin(1:split-1);
end

% Get the dimension
nd = size(X,2);
edges = varargin;
if nd<length(edges)
    nd = length(edges); % wasting CPU time warranty
else
    edges(end+1:nd) = {DEFAULT_NBINS};
end

% Allocation of array loc: index location of X in the bins
loc = zeros(size(X));
sz = zeros(1,nd);
% Loop in the dimension
for d=1:nd
    ed = edges{d};
    Xd = X(:,d);
    if isempty(ed)
        ed = DEFAULT_NBINS;
    end
    if isscalar(ed) % automatic linear subdivision
        ed = linspace(min(Xd),max(Xd),ed+1);
    end
    edges{d} = ed;
    % Call histc on this dimension
    [dummy loc(:,d)] = histc(Xd, ed, 1);
    % Use sz(d) = length(ed); to create consistent number of bins
    sz(d) = length(ed)-1;
end % for-loop

% Clean
clear dummy

% This is need for seldome points that hit the right border
sz = max([sz; max(loc,[],1)]);

% Compute the mid points
mid = cellfun(@(e) 0.5*(e(1:end-1)+e(2:end)), edges, ...
              'UniformOutput', false);

% Count for points where all coordinates are falling in a corresponding
% bins
if nd==1
    sz = [sz 1]; % Matlab doesn't know what is one-dimensional array!
end

hasdata = all(loc>0, 2);
if ~isempty(AccumData)
    count = accumarray(loc(hasdata,:), AccumData(hasdata), sz, Fun{:});
else
    count = accumarray(loc(hasdata,:), 1, sz);
end

return

end % histcn

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% ASTROTIK by Francesco Santilli
% f2v computes velocity vector of an orbit at different true anomalies
%
% Usage: V = f2v(orbit,f,mi)
%
% where: orbit = [p e i o w] = 3d orbit elements
%        orbit = [p e w]     = 2d orbit elements
%           p = semi-latus rectum [L] (p>0)
%           e = eccentricity [-] (e>=0)
%           i = inclinaion [rad]
%           o = raan [rad]
%           w = argument of perifocus [rad]
%        f(k) = true anomalies [rad]
%        mi = gravitational parameter [L^3*T^-2] (mi>0)
%        V(k,:) = velocity [L*T^-1]

function V = f2v(orbit,f,mi)

    if ~(nargin == 3)
        error('Wrong number of input arguments.')
    end

    D = check(orbit,1);
    K = check(f,1);
    check(mi,0)

    if ~(D==3 || D==5)
        error('Wrong size of input arguments.')
    end

    d3 = (D==5);
    f = f(:);
    cf = cos(f);
    sf = sin(f);

    % rotation matrix
    if d3
        [p,e,i,o,w] = take(orbit);
        R = rotation([o i w]);
        R = R(1:2,:);
    else
        [p,e,w] = take(orbit);
        cw = cos(w);
        sw = sin(w);
        R = [ cw sw
             -sw cw];
    end

    if mi <= 0
        error('mi must be a strictly positive value.')
    end

    if p <= 0
        error('p must be a stricly positive value.')
    end

    if e < 0
        error('e must be a positive value.')
    end

    % velocity
    v = sqrt(mi/p);
    V = v * [-sf e+cf] * R;

end


text:
function [dS,t,f]=mtdspecgramc(data,movingwin,phi,params)
% Multi-taper derivative of the time-frequency spectrum - continuous process
%
% Usage:
%
% [dS,t,f]=mtdspecgramc(data,movingwin,phi,params)
% Input:
%   Note that all times can be in arbitrary units. But the units have to be
%   consistent. So, if E is in secs, win, t have to be in secs, and Fs has to
%   be Hz. If E is in samples, so are win and t, and Fs=1. In case of spike
%   times, the units have to be consistent with the units of data as well.
%
%       data        (in form samples x channels/trials or a single vector) -- required
%       movingwin         (in the form [window winstep] i.e length of moving
%                                                 window and step size.
%                                                 Note that units here have
%                                                 to be consistent with
%                                                 units of Fs - required
%       phi         (angle for evaluation of derivative) -- required
%                       e.g. phi=[0,pi/2] giving the time and frequency
%                       derivatives
%       params: structure with fields tapers, pad, Fs, fpass, trialave
%       -optional
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                    (1) A numeric vector [TW K] where TW is the
%                        time-bandwidth product and K is the number of
%                        tapers to be used (less than or equal to
%                        2TW-1).
%                    (2) A numeric vector [W T p] where W is the
%                        bandwidth, T is the duration of the data and p
%                        is an integer such that 2TW-p tapers are used. In
%                        this form there is no default i.e. to specify
%                        the bandwidth, you have to specify T and p as
%                        well. Note that the units of W and T have to be
%                        consistent: if W is in Hz, T must be in seconds
%                        and vice versa. Note that these units must also
%                        be consistent with the units of params.Fs: W can
%                        be in Hz if and only if params.Fs is in Hz.
%                        The default is to use form 1 with TW=3 and K=5
%                     Note that T has to be equal to movingwin(1).
%
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%           trialave - (average over trials/channels when 1, don't average when 0) - optional. Default 0
% Output:
%       dS      (spectral derivative in form phi x time x frequency x channels/trials if trialave=0;
%               in form phi x time x frequency if trialave=1)
%       t       (times)
%       f       (frequencies)

if nargin < 3; error('Need data, window parameters and angle'); end;
if nargin < 4; params=[]; end;

if length(params.tapers)==3 & movingwin(1)~=params.tapers(2);
    error('Duration of data in params.tapers is inconsistent with movingwin(1), modify params.tapers(2) to proceed')
end

[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear err
data=change_row_to_column(data);
[N,C]=size(data);
Nwin=round(Fs*movingwin(1)); % number of samples in window
Nstep=round(movingwin(2)*Fs); % number of samples to step through
nfft=max(2^(nextpow2(Nwin)+pad),Nwin);
f=getfgrid(Fs,nfft,fpass); Nf=length(f);
params.tapers=dpsschk(tapers,Nwin,Fs); % check tapers
params.tapers=tapers;
winstart=1:Nstep:N-Nwin+1;
nw=length(winstart);
if trialave==0; dS=zeros(length(phi),nw,Nf,C); else dS=zeros(length(phi),nw,Nf); end;
for n=1:nw;
   indx=winstart(n):winstart(n)+Nwin-1;
   datawin=data(indx,:);
   [ds,f]=mtdspectrumc(datawin,phi,params);
   dS(:,n,:,:)=ds;
end;
dS=squeeze(dS);
sz=size(dS);
% if length(sz)==3;
%    dS=permute(dS,[2 1 3 4]);
% elseif length(phi)>1
%    dS=permute(dS,[2 1 3]);
% end;
winmid=winstart+round(Nwin/2);
t=winmid/Fs;


text:
function A00 = synA00Qmin(A11, sizeA00, cmax)
%-----------------------------------------------------------------------------
%
% For each point of colour 00 this function assigns the minimum value at the
% neighbouring gridpoints of colour 11.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: December 12, 2001.
% (c) 1998-2002 Stichting CWI, Amsterdam
%-----------------------------------------------------------------------------
[n11, m11]=size(A11);
if nargin == 3
  o=[0 0];
  if ~all(size(o) == size(sizeA00))
    error(' synA00Qmin - unexpected dimensions of sizeA00 ')
  else
    clear o;
    n00=sizeA00(1);
    m00=sizeA00(2);
  end
elseif nargin == 2
  n00=n11+1;
  m00=m11+1;
else
  error(' synA00Qmin - number of arguments should be either 2 or 3 ')
end
%[n00, m00]=size(A00);
if m00 == m11+1
  if n00 == n11+1
    S = min(extL(A11, cmax), extR(A11, cmax));
    A00=min(extD(S, cmax), extU(S, cmax));
  elseif n00 == n11
    S = min(extL(A11, cmax), extR(A11, cmax));
    A00=min(S, stripD(extU(S, cmax)));
  else
    disp([' size A11 = ' int2str(size(A11)) ' size A00 = ' int2str([n00 m00])]);
    error(' synA00Qmin - A11 and target A00 do not match ');
  end
elseif m00 == m11
  if n00 == n11+1
    S = min(stripR(extL(A11, cmax)), A11);
    A00=min(extD(S, cmax), extU(S, cmax));
  elseif n00 == n11
    S = min(stripR(extL(A11, cmax)), A11);
    A00=min(S, stripD(extU(S, cmax)));
  else
    disp([' size A11 = ' int2str(size(A11)) ' size A00 = ' int2str([n00 m00])]);
    error(' synA00Qmin - A11 and target A00 do not match ');
  end
else
  disp([' size A11 = ' int2str(size(A11)) ' size A00 = ' int2str([n00 m00])]);
  error(' synA00Qmin - A11 and target A00 do not match ');
end
%------------------------------------------------------------------------------


text:
function varargout = information(R, opts, varargin)

%INFORMATION Computes mutual information using different methods and
% different bias correction procedures.
%
%   ------
%   SYNTAX
%   ------
%       [...] = information(R, opts, output list...)
%
%   ---------
%   ARGUMENTS
%   ---------
%   R           - Response matrix.
%   opts        - Options structure.
%   output list - List of strings specifying what to compute.
%
%   -------------------
%   THE RESPONSE MATRIX
%   -------------------
%   L-dimensional responses to S distinct stimuli are stored in a response
%   matrix R of size L-by-T-by-S, T being the maximum number of trials
%   available for any of the stimuli. Emtpy trials, i.e., elements of R not
%   corresponding to a recorded response, can take any value.
%
%   ---------------------
%   THE OPTIONS STRUCTURE
%   ---------------------
%   The options structure can include any the following fields:
%
%   opts.nt
%   -------
%       This field specifies the number of trials (responses) recorded for
%       each stimulus. It can be either a scalar (for constant number of
%       trials per stimulus) or an array of length S.
%
%       NT must satisfy the following two conditions:
%       - max(nt) = T
%       - length(nt) = S (if nt is an array)
%
%   opts.method
%   -----------
%       This field specifies which estimation method to use and can be one
%       of the following strings:
%
%       --------------------------
%       | 'dr' | Direct method   |
%       | 'gs' | Gaussian method |
%       --------------------------
%
%       IMPORTANT!
%       ==========
%       The direct method requires the response values to be discretized
%       into non-negative integer values (this is meant only in a numerical
%       sense, the MATLAB variable still needs to be of type double). See
%       function BINR.M for instruction on how to discretize the responses.
%       Failing to properly discretizing the response will result in Matlab
%       crashing.
%
%   opts.bias
%   ---------
%       This field specifies the bias correction procedure. It can be one
%       of the following strings:
%
%       -------------------------------------
%       | 'qe'    | Quadratic EXtrapolation |
%       | 'pt'    | Panzeri & Treves 1996   |
%       | 'gsb'   | Gaussian bias           |
%       | 'naive' | Biased naive estimates  |
%       -------------------------------------
%
%   opts.btsp (optional)
%   --------------------
%       This field must be a (non-negative) scalar specifying how many
%       bootstrap estimates to compute.
%
%       Bootstrap estimates are performed by means of pairing stimuli and
%       responses at random and computing the entropy quantities for these
%       random pairings; each estimate corresponds to a different random
%       pairing configuration.
%
%       See the examples below for additional information on how to use
%       this option.
%
%       DEFAULT: 0.
%
%   opts.verbose (optional)
%   -----------------------
%       If this field exists and is set to true a summary of the selected
%       options is displayed and additional checks are performed on the
%       input variables. No warnings are displayed unless this options is
%       enabled.
%
%       This feature is useful to check whether INFORMATION is being called
%       correctly. It is therefore highly reccomended for new users or when
%       first running of the program with new input options. However, keep
%       in mind that these checks drammatically increases computation time
%       and are thus not reccommended for computationally intensive
%       session.
%
%       DEFAULT: false.
%
%   ---------------
%   THE OUTPUT LIST
%   ---------------
%   To specify which IT quantities need to compute, one or more of the
%   following strings has to be specified:
%
%       =================================================================
%       | Option  | Description       | Expression (in terms of ENTROPY |
%       |         |                   | output options)                 |
%       =================================================================
%       | 'I'     | I(R;S)            | I     = HR - HRS                |
%       | 'Ish'   | I(R;S) shuffle    | Ish   = HR - HiRS + HshRS - HRS |
%       |---------------------------------------------------------------|
%       | 'IX'    | I(R)              | IX    = HlR - HR                |
%       |---------------------------------------------------------------|
%       | 'ILIN'  | I_lin(S;R)        | ILIN  = HlR - HiRS              |
%       |---------------------------------------------------------------|
%       | 'SYN'   | Syn               | SYN   = HR - HRS - HlR + HiRS   |
%       | 'SYNsh' | Syn shuffle       | SYNsh = HR + HshRS - HRS - HlR  |
%       |---------------------------------------------------------------|
%       | 'ISS'   | I_sig_sim         | ISS   = HiR - HlR               |
%       |---------------------------------------------------------------|
%       | 'IC'    | I_cor             | IC    = HR - HRS + HiRS - HiR   |
%       | 'ICsh'  | I_cor shuffle     | ICsh  = HR + HshRS - HRS - HiR  |
%       |---------------------------------------------------------------|
%       | 'ICI'   | I_cor_ind         | ICI   = ChiR - HiR              |
%       |---------------------------------------------------------------|
%       | 'ICD'   | I_cor_dep         | ICD   = HR - HRS - ChiR + HiRS  |
%       | 'ICDsh' | I_cor_dep shuffle | ICDsh = HR + HshRS - HRS - ChiR |
%       |---------------------------------------------------------------|
%       | 'ILB1'  | I_LB1             | ILB1  = HR - HiRS               |
%       |---------------------------------------------------------------|
%       | 'ILB2'  | I_LB2             | ILB2  = ChiR - HiRS             |
%       =================================================================
%
%   Outputs are returned IN THE SAME ORDER as that specified in the output
%   list.
%
%   IMPORTANT: Not all combinations of method, bias and output options are
%   possible. For example, bias correction 'pt' can only be used together
%   with method 'dr'. The allowed combinations of method, bias and output
%   options are summarized in the following tables:
%
%       =============================================
%       | DIRECT METHOD                             |
%       =============================================
%       |         | 'naive' | 'qe'  | 'pt'  | 'gsb' |
%       |-------------------------------------------|
%       | 'I'     |    X    |   X   |   X   |   -   |
%       | 'Ish'   |    X    |   X   |   X   |   -   |
%       |---------|---------|-------|-------|-------|
%       | 'IX'    |    X    |   X   |   X   |   -   |
%       |---------|---------|-------|-------|-------|
%       | 'ILIN'  |    X    |   X   |   X   |   -   |
%       |---------|---------|-------|-------|-------|
%       | 'SYN'   |    X    |   X   |   X   |   -   |
%       | 'SYNsh' |    X    |   X   |   X   |   -   |
%       |---------|---------|-------|-------|-------|
%       | 'ISS'   |    X    |   X   |   -   |   -   |
%       |---------|---------|-------|-------|-------|
%       | 'IC'    |    X    |   X   |   -   |   -   |
%       | 'ICsh'  |    X    |   X   |   -   |   -   |
%       |---------|---------|-------|-------|-------|
%       | 'ICI'   |    X    |   X   |   -   |   -   |
%       |---------|---------|-------|-------|-------|
%       | 'ICD'   |    X    |   X   |   -   |   -   |
%       | 'ICDsh' |    X    |   X   |   -   |   -   |
%       |---------|---------|-------|-------|-------|
%       | 'ILB1'  |    X    |   X   |   X   |   -   |
%       |---------|---------|-------|-------|-------|
%       | 'ILB2'  |    X    |   X   |   -   |   -   |
%       =============================================
%
%   Legend: X: combination available
%           -: combination NOT permitted
%
%
%       =============================================
%       | GAUSSIAN METHOD                           |
%       =============================================
%       |         | 'naive' | 'qe'  | 'pt'  | 'gsb' |
%       |-------------------------------------------|
%       | 'I'     |    X    |  n.r. |   -   |   X   |
%       | 'Ish'   |    X    |  n.r. |   -   |   X   |
%       |---------|---------|-------|-------|-------|
%       | 'IX'    |    X    |  n.r. |   -   |   X   |
%       |---------|---------|-------|-------|-------|
%       | 'ILIN'  |    X    |  n.r. |   -   |   X   |
%       |---------|---------|-------|-------|-------|
%       | 'SYN'   |    X    |  n.r. |   -   |   X   |
%       | 'SYNsh' |    X    |  n.r. |   -   |   X   |
%       |---------|---------|-------|-------|-------|
%       | 'ISS'   |   NaN   |  NaN  |   -   |  NaN  |
%       |---------|---------|-------|-------|-------|
%       | 'IC'    |   NaN   |  NaN  |   -   |  NaN  |
%       | 'ICsh'  |   NaN   |  NaN  |   -   |  NaN  |
%       |---------|---------|-------|-------|-------|
%       | 'ICI'   |   NaN   |  NaN  |   -   |  NaN  |
%       |---------|---------|-------|-------|-------|
%       | 'ICD'   |   NaN   |  NaN  |   -   |  NaN  |
%       | 'ICDsh' |   NaN   |  NaN  |   -   |  NaN  |
%       |---------|---------|-------|-------|-------|
%       | 'ILB1'  |    X    |  n.r. |   -   |   X   |
%       |---------|---------|-------|-------|-------|
%       | 'ILB2'  |   NaN   |  NaN  |   -   |  NaN  |
%       =============================================
%
%   Legend: X   : combination available
%           -   : combination NOT permitted
%           n.r.: combination available but not recommended
%           NaN : NaN returned
%
%   --------
%   EXAMPLES
%   --------
%   In the following examples, we assume R to be a 2-by-10-by-3 matrix
%   i.e., R stores 2-dimensional responses to 3 different stimuli. We also
%   assume that, while 10 trials are available for stimulus 1 and 2, only 7
%   trials have been recorded for stimulus 3.
%
%   - Estimate I(S;R) using the direct method and no bias corrections
%
%       opts.nt = [10 10 7];
%       opts.method = 'dr';
%       opts.bias = 'naive';
%       X = information(R, opts, 'I');
%
%   - Estimate I(S;R) and I_shuffle(S;R) using direct method and the
%     quadratic extrapolation bias correction
%
%       opts.nt = [10 10 7];
%       opts.method = 'dr';
%       opts.bias = 'gsb';
%       [X, Y] = information(R, opts, 'Ish', 'I');
%
%     where the estimate of I_shuffle(S;R) is stored in the X and that of
%     I(S;R) in Y.
%
%   - Compute gaussian naive estimate of I(S;R) together with 20 bootstrap
%     estimates:
%
%       opts.nt = [10 10 7];
%       opts.method = 'gs';
%       opts.bias = 'naive';
%       opts.btsp = 20;
%       X = information(R, opts, 'I');
%
%     Note that, in this case, Y is an array of size 21-by-1: Y(1) gives
%     the estimate for I(S;R) computed using the input matrix R; Y(2:21)
%     are are 20 distinct bootstrap estimates of I(S;R).
%
%   -------
%   REMARKS
%   -------
%   - Field-names in the option structure are case-sensitive
%
%   - Ouput options are case INsensitive
%
%   - It is more efficient to call INFORMATION with several output options
%     rather than calling the function repeatedly. For example:
%
%         [X, Y] = information(R, opts, 'I', 'Ish');
%
%     is faster than
%
%         X = information(R, opts, 'I');
%         Y = information(R, opts, 'Ish');
%
%   - Some MEX files in the toolbox create static arrays which are used
%     to store computations performed in previous calls to the routines.
%     This memory is freed automatically when Matlab is quitted. However,
%     consider using
%
%         clear mex;
%
%     when needing to free all of Matlab's available memory.

% NOTE:
% The function also computes IXS and IXSsh according to the following
% equations:
%
%       =================================================================
%       | Option  | Description       | Expression (in terms of ENTROPY |
%       |         |                   | output options)                 |
%       =================================================================
%       | 'IXS'   | I(R|S)            | IXS   = HiRS - HRS              |
%       | 'IXSsh' | I(R|S) shuffle    | IXSsh = HshRS - HRS             |
%       =================================================================
%
% However, since this quantity is meaningful only for L=2, this feature
% is kept hidden. The combination tables for this quantity are
% as follows:
%
%       =============================================
%       | DIRECT METHOD                             |
%       =============================================
%       |         | 'naive' | 'qe'  | 'pt'  | 'gsb' |
%       |-------------------------------------------|
%       | 'IXS'   |    X    |   X   |   X   |   -   |
%       | 'IXSsh' |    X    |   X   |   X   |   -   |
%       =============================================
%
%   Legend: X: combination available
%           -: combination NOT available
%
%
%       =============================================
%       | GAUSSIAN METHOD                           |
%       =============================================
%       |         | 'naive' | 'qe'  | 'pt'  | 'gsb' |
%       |-------------------------------------------|
%       | 'IXS'   |    X    |  n.r. |   -   |   X   |
%       | 'IXSsh' |    X    |  n.r. |   -   |   X   |
%       =============================================
%
%   Legend: X   : combination available
%           -   : combination NOT available
%           n.r.: combination available but not recommended

%   Copyright (C) 2009 Cesare Magri
%   Version: 1.0.5

% -------
% LICENSE
% -------
% This software is distributed free under the condition that:
%
% 1. it shall not be incorporated in software that is subsequently sold;
%
% 2. the authorship of the software shall be acknowledged and the following
%    article shall be properly cited in any publication that uses results
%    generated by the software:
%
%      Magri C, Whittingstall K, Singh V, Logothetis NK, Panzeri S: A
%      toolbox for the fast information analysis of multiple-site LFP, EEG
%      and spike train recordings. BMC Neuroscience 2009 10(1):81;
%
% 3.  this notice shall remain in place in each source file.

% NOTE:
% In this function HiRS is indeed output option 'HiRS' of function
% entropy.m, thus corresponding to the variable HlRS in the same function.

whereI     = strcmpi(varargin, 'i');
whereIsh   = strcmpi(varargin, 'ish');
whereIX    = strcmpi(varargin, 'ix');
whereIXS   = strcmpi(varargin, 'ixs');
whereIXSsh = strcmpi(varargin, 'ixssh');
whereILIN  = strcmpi(varargin, 'ilin');
whereSYN   = strcmpi(varargin, 'syn');
whereSYNsh = strcmpi(varargin, 'synsh');
whereISS   = strcmpi(varargin, 'iss');
whereIC    = strcmpi(varargin, 'ic');
whereICsh  = strcmpi(varargin, 'icsh');
whereICI   = strcmpi(varargin, 'ici');
whereICD   = strcmpi(varargin, 'icd');
whereICDsh = strcmpi(varargin, 'icdsh');
whereILB1  = strcmpi(varargin, 'ilb1');
whereILB2  = strcmpi(varargin, 'ilb2');

whereuvar  = strcmpi(varargin, 'uvar');
if sum(whereuvar)
  Runit = varargin{find(whereuvar)+1};
  sel = find(whereuvar);
  sel = [sel sel+1];
  varargin(sel) = [];
end

doI     = any(whereI);
doIsh   = any(whereIsh);
doIX    = any(whereIX);
doIXS   = any(whereIXS);
doIXSsh = any(whereIXSsh);
doILIN  = any(whereILIN);
doSYN   = any(whereSYN);
doSYNsh = any(whereSYNsh);
doISS   = any(whereISS);
doIC    = any(whereIC);
doICsh  = any(whereICsh);
doICI   = any(whereICI);
doICD   = any(whereICD);
doICDsh = any(whereICDsh);
doILB1  = any(whereILB1);
doILB2  = any(whereILB2);

% Checks ------------------------------------------------------------------
specifiedOutputOptsVec = ...
    [doI doIsh doIX doIXS doIXSsh doILIN doSYN doSYNsh doISS doIC doICsh doICI doICD doICDsh doILB1 doILB2];
NspecifiedOutputOpts = sum(specifiedOutputOptsVec);
lengthVarargin = length(varargin);
if NspecifiedOutputOpts~=lengthVarargin
    msg = 'Unknown selection or repeated option in output list.';
    error('information:unknownOutputOpt', msg);
end

% Restrictions on possible combinantions ----------------------------------
if strcmpi(opts.method, 'dr')
    % Can't apply bias-correction gsb with method dr:
    if strcmpi(opts.bias, 'gsb')
        msg = 'Bias correction ''gsb'' can only be used in conjunction with method ''gs''.';
        error('Information:drMethodAndGsbBias', msg);
    end

    % Can't compute ISS, IC, ICsh, ICI, ICD, ICDsh or ILB2 for
    % bias-correction pt
    if strcmpi(opts.bias, 'pt') && (doISS || doIC || doICsh || doICI || doICD || doICDsh || doILB2)
        msg = 'One or more of the selected output options are not available for bias correction ''pt''.';
        error('information:ptBiasAndNonAvailableOutputOpt', msg);
    end
end

% Default verbose value:
if ~isfield(opts, 'verbose')
    opts.verbose = false;
end;

isGaussianMethod = false;
if strcmpi(opts.method, 'gs')
    isGaussianMethod = true;

    % Gaussian and QE is not recommended:
    if opts.verbose && strcmpi(opts.bias, 'qe')
        msg = 'Usage of bias correction ''qe'' in conjunction with gaussian method is not recommended.';
        warning('Information:gsMethodAndQeBias', msg);
    end

    % Gaussian and PT is not allowed:
    if strcmpi(opts.bias, 'pt')
        msg = 'Bias correction ''pt'' can only be used in conjunction with method ''dr''.';
        error('Information:gsMethodAndPtBias', msg);
    end
end

allOutputOpts = {'HR' 'HRS' 'HlR' 'HiR' 'HiRS' 'ChiR' 'HshRS'};
positionInOuputOptsList = 0;

% What needs to be computed by ENTTROPY -----------------------------------
% Need to compute H(R)?
doHR = false;
if doI || doIsh || doIX || doSYN || doSYNsh || doIC || doICsh || doICD || doICDsh || doILB1
    positionInOuputOptsList = positionInOuputOptsList + 1;
    whereHR = positionInOuputOptsList;
    doHR = true;
end

% Need to compute H(R|S)?
doHRS = false;
if doI || doIsh || doIXS || doIXSsh || doSYN || doSYNsh || doIC || doICsh || doICD || doICDsh
    positionInOuputOptsList = positionInOuputOptsList + 1;
    whereHRS = positionInOuputOptsList;
    doHRS = true;
end

% Need to compute H_lin(R)?
doHlR = false;
if doIX || doSYN || doILIN || doSYNsh || doISS
    positionInOuputOptsList = positionInOuputOptsList + 1;
    whereHlR = positionInOuputOptsList;
    doHlR = true;
end

% Need to compute H_ind(R)?
doHiR = false;
if (doISS || doIC || doICsh || doICI) && ~isGaussianMethod
    positionInOuputOptsList = positionInOuputOptsList + 1;
    whereHiR = positionInOuputOptsList;
    doHiR = true;
end

% Need to compute H_ind(R|S)?
doHiRS = false;
if doIsh || doIXS || doILIN || doSYN || doIC || doICD || doILB1 || doILB2
    positionInOuputOptsList = positionInOuputOptsList + 1;
    whereHiRS = positionInOuputOptsList;
    doHiRS = true;
end

% Need to compute Chi(R)?
doChiR = false;
if (doICI || doICD || doICDsh || doILB2) && ~isGaussianMethod
    positionInOuputOptsList = positionInOuputOptsList + 1;
    whereChiR = positionInOuputOptsList;
    doChiR = true;
end

% Need to compute H_sh(R|S)?
doHshRS = false;
if doIsh || doIXSsh || doSYNsh || doICsh || doICDsh
    positionInOuputOptsList = positionInOuputOptsList + 1;
    whereHshRS = positionInOuputOptsList;
    doHshRS = true;
end

% Computing information theoretic quantities ------------------------------
outputOptsList = allOutputOpts([doHR doHRS doHlR doHiR doHiRS doChiR doHshRS]);

H = cell(positionInOuputOptsList, 1);
if ~exist('Runit', 'var')
  [H{:}] = entropy(R, opts, outputOptsList{:});
else
  [H{:}] = entropy(R, opts, outputOptsList{:}, 'uvar', Runit);
end

% Assigning output --------------------------------------------------------
varargout = cell(length(varargin),1);

% I = HR - HRS
if doI
    varargout(whereI) = {H{whereHR} - H{whereHRS}};
end

% Ish = HR - HiRS + HshRS - HRS
if doIsh
    varargout(whereIsh) = {H{whereHR} - H{whereHiRS} + H{whereHshRS} - H{whereHRS}};
end

% IX = HlR - HR
if doIX
    varargout(whereIX) = {H{whereHlR} - H{whereHR}};
end

% IXS = HiRS - HRS
if doIXS
    varargout(whereIXS) = {H{whereHiRS} - H{whereHRS}};
end

% IXSsh = HshRS - HRS
if doIXSsh
    varargout(whereIXSsh) = {H{whereHshRS} - H{whereHRS}};
end

% ILIN = HlR - HiRS
if doILIN
    varargout(whereILIN) = {H{whereHlR} - H{whereHiRS}};
end

% SYN = HR - HRS - HlR + HiRS
if doSYN
    varargout(whereSYN) = {H{whereHR} - H{whereHRS} - H{whereHlR} + H{whereHiRS}};
end

% SYNsh = HR + HshRS - HRS - HlR
if doSYNsh
    varargout(whereSYNsh) = {H{whereHR} + H{whereHshRS} - H{whereHRS} - H{whereHlR}};
end

% ISS = HiR - HlR
if doISS
    if ~isGaussianMethod
        varargout(whereISS) = {H{whereHiR} - H{whereHlR}};
    else
        varargout(whereISS) = {NaN};
    end
end

% IC = HR - HRS + HiRS - HiR
if doIC
    if ~isGaussianMethod
        varargout(whereIC) = {H{whereHR} - H{whereHRS} + H{whereHiRS} - H{whereHiR}};
    else
        varargout(whereIC) = {NaN};
    end
end

% ICsh  = HR + HshRS - HRS - HiR
if doICsh
    if ~isGaussianMethod
        varargout(whereICsh) = {H{whereHR} + H{whereHshRS} - H{whereHRS} - H{whereHiR}};
    else
        varargout(whereICsh) = {NaN};
    end
end

% ICI= ChiR - HiR
if doICI
    if ~isGaussianMethod
        varargout(whereICI) = {H{whereChiR} - H{whereHiR}};
    else
        varargout(whereICI) = {NaN};
    end
end

% ICD = HR - HRS - ChiR + HiRS
if doICD
    if ~isGaussianMethod
        varargout(whereICD) = {H{whereHR} - H{whereHRS} - H{whereChiR} + H{whereHiRS}};
    else
        varargout(whereICD) = {NaN};
    end
end

% ICDsh = HR + HshRS - HRS - ChiR
if doICDsh
    if ~isGaussianMethod
        varargout(whereICDsh) = {H{whereHR} + H{whereHshRS} - H{whereHRS} - H{whereChiR}};
    else
        varargout(whereICDsh) = {NaN};
    end
end

% ILB1 = HR - HiRS
if doILB1
    varargout(whereILB1) = {H{whereHR} - H{whereHiRS}};
end

% ILB2 = ChiR - HiRS
if doILB2
    if ~isGaussianMethod
        varargout(whereILB2) = {H{whereChiR} - H{whereHiRS}};
    else
        varargout(whereILB2) = {NaN};
    end
end


text:
function value = p21_i4 ( action, name, value )

%*****************************************************************************80
%
%% P21_I4 sets or gets I4 parameters for problem 21.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'D' to set a parameter to its default value.
%    'S' to set a parameter,
%    'G' to get a parameter,
%    'I' to increment a parameter.
%
%    Input, string NAME, the name of the variable.
%    '#' is the number of calls to the integrand routine.
%
%    Input/output, integer VALUE.
%    * If ACTION = 'I', then VALUE is an input quantity, and is the
%      new value to be added to NAME.
%    * If ACTION = 'S', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'G', then VALUE is an output quantity, and is the
%      current value of NAME.
%
  persistent calls;

  if ( action == 'D' | action == 'd' )

    if ( name == '#' | name == '*' )
      calls = 0;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == '#' )
      value = calls;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P21_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P21_I4 - Fatal error!' );
    end

  elseif ( action == 'I' | action == 'i' )

    if ( name == '#' )
      calls = calls + value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P21_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P21_I4 - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == '#' )
      calls = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P21_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P21_I4 - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P21_I4 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P21_I4 - Fatal error!' );

  end

  return
end


text:
% 3D Lattice Boltzmann (BGK) model of a fluid.
% D3Q19 model. At each timestep, particle densities propagate
% outwards in the directions indicated in the figure. An
% equivalent 'equilibrium' density is found, and the densities
% relax towards that state, in a proportion governed by omega.
%               Iain Haslam, March 2006.
nx=12;ny=nx;nz=nx; omega=1.0; density=1.0;t1=1/3; t2=1/18; t3=1/36;
F=repmat(density/19,[nx ny nz 19]); FEQ=F; matsize=nx*ny*nz;
CI=[0:matsize:matsize*19];
BOUND=zeros(nx,ny,nz);
for i=1:nx, for j=1:ny, for k=1:nz
	BOUND(i,j,k)=((i-5)^2+(j-6)^2+(k-7)^2)<6;
end, end, end
BOUND(:,:,1)=1;BOUND(:,1,:)=1;
ON=find(BOUND); %matrix offset of each Occupied Node
TO_REFLECT=[ON+CI(2) ON+CI(3) ON+CI(4) ON+CI(5)	ON+CI(6) ON+CI(7) ON+CI(8) ...
  ON+CI(9) ON+CI(10) ON+CI(11) ON+CI(12) ON+CI(13) ON+CI(14) ON+CI(15) ...
  ON+CI(16) ON+CI(17) ON+CI(18) ON+CI(19)];
REFLECTED=[ON+CI(3) ON+CI(2) ON+CI(5) ON+CI(4) ON+CI(7) ON+CI(6) ON+CI(11) ...
  ON+CI(10) ON+CI(9) ON+CI(8) ON+CI(15) ON+CI(14) ON+CI(13) ON+CI(12) ...
  ON+CI(19) ON+CI(18) ON+CI(17) ON+CI(16)];
avu=1; prevavu=1; ts=0; deltaU=1e-7; numactivenodes=sum(sum(sum(1-BOUND)));
while (ts<4000 & 1e-10<abs((prevavu-avu)/avu)) | ts<100
	% Propagate
	%nearest-neighbours
	F(:,:,:,2)=F(:,:,[nz 1:nz-1],2);
	F(:,:,:,3)=F(:,:,[2:nz 1],3);
	F(:,:,:,4)=F(:,[ny 1:ny-1],:,4);
	F(:,:,:,5)=F(:,[2:ny 1],:,5);
	F(:,:,:,6)=F([nx 1:nx-1],:,:,6);
	F(:,:,:,7)=F([2:nx 1],:,:,7);
	%next-nearest neighbours
	F(:,:,:,8)= F([nx 1:nx-1],[ny 1:ny-1],:,8);
	F(:,:,:,9)= F([nx 1:nx-1],[2:ny 1],:,9);
	F(:,:,:,10)=F([2:nx 1],[ny 1:ny-1],:,10);
	F(:,:,:,11)=F([2:nx 1],[2:ny 1],:,11);
	F(:,:,:,12)=F([nx 1:nx-1],:,[nz 1:nz-1],12);
	F(:,:,:,13)=F([nx 1:nx-1],:,[2:nz 1],13);
	F(:,:,:,14)=F([2:nx 1],:,[nz 1:nz-1],14);
	F(:,:,:,15)=F([2:nx 1],:,[2:nz 1],15);
	F(:,:,:,16)=F(:,[ny 1:ny-1],[nz 1:nz-1],16);
	F(:,:,:,17)=F(:,[ny 1:ny-1],[2:nz 1],17);
	F(:,:,:,18)=F(:,[2:ny 1],[nz 1:nz-1],18);
	F(:,:,:,19)=F(:,[2:ny 1],[2:nz 1],19);
	BOUNCEDBACK=F(TO_REFLECT); %Densities bouncing back at next timestep
	% Relax; calculate equilibrium state (FEQ) with equivalent speed and density to F
	DENSITY = sum(F,4);
	UX=(sum(F(:,:,:,[6 8 9 12 13]),4)-sum(F(:,:,:,[7 10 11 14 15]),4))./DENSITY;
	UY=(sum(F(:,:,:,[4 8 10 16 17]),4)-sum(F(:,:,:,[5 9 11 18 19]),4))./DENSITY;
	UZ=(sum(F(:,:,:,[2 12 14 16 18]),4)-sum(F(:,:,:,[3 13 15 17 19]),4))./DENSITY;
	UX(1,:,:)=UX(1,:,:)+deltaU; %Increase inlet pressure
	UX(ON)=0; UY(ON)=0; UZ(ON)=0; DENSITY(ON)=0; U_SQU=UX.^2+UY.^2+UZ.^2;
	U8=UX+UY;U9=UX-UY;U10=-UX+UY;U11=-U8;U12=UX+UZ;U13=UX-UZ;
	U14=-U13;U15=-U12;U16=UY+UZ;U17=UY-UZ;U18=-U17;U19=-U16;
	% Calculate equilibrium distribution: stationary
	FEQ(:,:,:,1)=t1*DENSITY.*(1-3*U_SQU/2);
	% nearest-neighbours
	FEQ(:,:,:,2)=t2*DENSITY.*(1 + 3*UZ + 9/2*UZ.^2 - 3/2*U_SQU);
	FEQ(:,:,:,3)=t2*DENSITY.*(1 - 3*UZ + 9/2*UZ.^2 - 3/2*U_SQU);
	FEQ(:,:,:,4)=t2*DENSITY.*(1 + 3*UY + 9/2*UY.^2 - 3/2*U_SQU);
	FEQ(:,:,:,5)=t2*DENSITY.*(1 - 3*UY + 9/2*UY.^2 - 3/2*U_SQU);
	FEQ(:,:,:,6)=t2*DENSITY.*(1 + 3*UX + 9/2*UX.^2 - 3/2*U_SQU);
	FEQ(:,:,:,7)=t2*DENSITY.*(1 - 3*UX + 9/2*UX.^2 - 3/2*U_SQU);
	% next-nearest neighbours
	FEQ(:,:,:,8) =t3*DENSITY.*(1 + 3*U8  + 9/2*(U8).^2  - 3*U_SQU/2);
	FEQ(:,:,:,9) =t3*DENSITY.*(1 + 3*U9  + 9/2*(U9).^2  - 3*U_SQU/2);
	FEQ(:,:,:,10)=t3*DENSITY.*(1 + 3*U10 + 9/2*(U10).^2 - 3*U_SQU/2);
	FEQ(:,:,:,11)=t3*DENSITY.*(1 + 3*U11 + 9/2*(U11).^2 - 3*U_SQU/2);
	FEQ(:,:,:,12)=t3*DENSITY.*(1 + 3*U12 + 9/2*(U12).^2 - 3*U_SQU/2);
	FEQ(:,:,:,13)=t3*DENSITY.*(1 + 3*U13 + 9/2*(U13).^2 - 3*U_SQU/2);
	FEQ(:,:,:,14)=t3*DENSITY.*(1 + 3*U14 + 9/2*(U14).^2 - 3*U_SQU/2);
	FEQ(:,:,:,15)=t3*DENSITY.*(1 + 3*U15 + 9/2*(U15).^2 - 3*U_SQU/2);
	FEQ(:,:,:,16)=t3*DENSITY.*(1 + 3*U16 + 9/2*(U16).^2 - 3*U_SQU/2);
	FEQ(:,:,:,17)=t3*DENSITY.*(1 + 3*U17 + 9/2*(U17).^2 - 3*U_SQU/2);
	FEQ(:,:,:,18)=t3*DENSITY.*(1 + 3*U18 + 9/2*(U18).^2 - 3*U_SQU/2);
	FEQ(:,:,:,19)=t3*DENSITY.*(1 + 3*U19 + 9/2*(U19).^2 - 3*U_SQU/2);
	F=omega*FEQ+(1-omega)*F;
	F(REFLECTED)=BOUNCEDBACK;
	prevavu=avu;avu=sum(sum(sum(UX)))/numactivenodes; ts=ts+1;
end
figure;zcut=5;colormap(gray(2));image(2-BOUND(:,:,5));hold on;
%Thanks to Thomas Wagner for correcting the transposed results plots
quiver(UY(:,:,zcut),UX(:,:,zcut));xlabel('y');ylabel('x');
title(['Flow field at z=',num2str(zcut),', after ',num2str(ts),'\deltat']);
figure;ycut=5;colormap(gray(2));image(2-squeeze(BOUND(:,ycut,:)));hold on;
quiver(squeeze(UZ(:,ycut,:)),squeeze(UX(:,ycut,:)));xlabel('z');ylabel('x');
title(['Flow field at y=',num2str(ycut),', after ',num2str(ts),'\deltat']);

repo: Aero-matlab


text:
function [] =  TestKShortestPath(case_number)
% Tests the function kShortestPath(case_number) for one of the 3 basic test cases
% specified by case_number. Then it displays the reults.
%==============================================================
% Meral Shirazipour
% DATE :           December 9 decembre 2009
% Last Updated:    April 2 2010 ; August 2 2011
%==============================================================
switch case_number
    case 1
        netCostMatrix = [inf 1 inf 1 ; 1 inf 1 1 ;inf 1 inf inf ;inf inf inf inf ];
        source=3;
        destination=4;
        k = 5;
    case 2
        netCostMatrix = [inf 5 10 15 20 inf; 5 inf 10 inf inf 5;10 10 inf 15 inf 10; 15 inf 15 inf 20 15;20 inf inf 20 inf 20; inf 5 10 15 20 inf];
        source=1;
        destination=6;
        k = 20;
    case 3
        netCostMatrix = [inf 10 50 15 inf; 10 inf 10 inf 10; 50 10 inf 5 5; 15 inf 5 inf 15; inf 10 5 15 inf];
        source=1;
        destination=5;
        k = 20;
    otherwise
        error('The only case options available are 1, 2 or 3');
end

%------Show case selected:------:
fprintf('You selected case #%d, a %dx%d network:\n',case_number,size(netCostMatrix,1),size(netCostMatrix,1));
disp(netCostMatrix);
fprintf('The path request is from source node %d to destination node %d, with K = %d \n',source,destination, k);

%------Call kShortestPath------:
[shortestPaths, totalCosts] = kShortestPath(netCostMatrix, source, destination, k);

%------Display results------:
fprintf('\nResult of Function call: kShortestPath(netCostMatrix, source, destination, k)  = \n\n');

if isempty(shortestPaths)
    fprintf('No path available between these nodes\n\n');
else
    for i = 1: length(shortestPaths)
        fprintf('Path # %d:\n',i);
        disp(shortestPaths{i})
        fprintf('Cost of path %d is %5.2f\n\n',i,totalCosts(i));
    end
end

end


text:
function dy = eulerODE01(x,y)
%
% Differential Equation:
%
% y'' + 6 * y' + 9 * y = 0
%
% With the initial conditions,
% y(1) = 5
% y'(1) = -1
% it has the solution:
%
% y = x^2 + 4/sqrt(x)
%
% Adams, Calculus of Several Variables, pp. 394-395
%
dy = zeros(2,1) ;
dy(1) = y(2) ;
dy(2) = -( 6*y(2) + 9*y(1) ) ;
%
% END
%


text:
function dt6FileName = dtiRawPreprocessSimple(dwRawFileName, bvalsFile, bvecsFile, t1FileName, dt6BaseName, numBootStrapSamples, dwOutMm)
%
% dtiRawPreprocessSimple([dwRawFileName], [bvalFile], [bvecFile],
%                        [t1FileName], [dt6BaseName='dti'],
%                        [numBootStrapSamples=500], [dwOutMm=[2 2 2]])
%
% A function for running a simple version of the mrDiffusion pre-processing
% steps on raw DWI data. This is similar to dtiRawPreprocess, but skips
% alignment to a t1 anatomical and assumes that you already have
% appropriate bvals/bvecs files.
%
% This assumes that you've already generated the raw NIFTI data file with
% something like dinifti.
%
% HISTORY:
% 2010.11.29 RFD wrote it.
%

% 1 = always clobber existing files, -1 = never replace existing files, just
% load them and use them, 0 = ask the user what to do for each file
clobber = -1;

outSuffix = '_aligned';
excludeVols = [];
bsplineInterpFlag = false;

if(~exist('dwOutMm','var')||isempty(dwOutMm))
    % The eddy-current-corrected, aligned diffusion images will be resampled at
    % 2mm isotropic voxels by default.
    dwOutMm = [2 2 2];
end

if(~exist('dt6BaseName','var')||isempty(dt6BaseName))
    dt6BaseName = '';
end

if(~exist('numBootStrapSamples','var')||isempty(numBootStrapSamples))
  numBootStrapSamples = 500;
end

mrDiffusionDir = fileparts(which('mrDiffusion.m'));

if(~exist('dwRawFileName','var')||isempty(dwRawFileName))
    [f,p] = uigetfile({'*.nii.gz';'*.*'},'Select a raw NIFTI file for input...');
    if(isnumeric(f)), disp('User canceled.'); return; end
    dwRawFileName = fullfile(p,f);
end

if(~exist(dwRawFileName,'file'))
  error(['Couldn''t find diffusion data file "' dwRawFileName '".']);
end

[dataDir,inBaseName] = fileparts(dwRawFileName);
[junk,inBaseName] = fileparts(inBaseName);
if(isempty(dataDir)), dataDir = pwd; end
mnB0Name = fullfile(dataDir,[inBaseName '_b0.nii.gz']);

outBaseName = [inBaseName outSuffix];
outBaseDir = fullfile(dataDir,outBaseName);
inBaseDir = fullfile(dataDir,inBaseName);
% Default output dir is one level above the dataDir. We assume that
% this is the 'subjectDir', which contains the 'raw' data dir.
subjectDir = fileparts(dataDir);
if(isempty(subjectDir)), subjectDir = pwd; end

if(~exist('t1FileName','var')||isempty(t1FileName))
    t1FileName = fullfile(subjectDir,'t1');
    [f,p] = uigetfile('*.nii.gz','Select the t1-file for alignment...',t1FileName);
    if(isnumeric(f)), error('user canceled.'); end
    t1FileName = fullfile(p,f);
elseif(strcmpi(t1FileName,'MNI'))
    t1FileName = fullfile(mrDiffusionDir,'templates','MNI_EPI.nii.gz');
    disp('The MNI EPI template will be used for alignment.');
end
if(~exist(t1FileName,'file'))
    error(['Couldn''t find t1 data file "' t1FileName '".']);
end

%bvalsFile = [inBaseDir '.bval'];
%bvecsFile = [inBaseDir '.bvec'];
if(~exist('bvalsFile','var')||isempty(bvalsFile))
    [f,p] = uigetfile({'*.bval';'*.bvals';'*.*'},'Select the bvals file...');
    if(isnumeric(f)), disp('User canceled.'); return; end
    bvalsFile = fullfile(p,f);
end
if(~exist('bvecsFile','var')||isempty(bvecsFile))
    [f,p] = uigetfile({'*.bvec';'*.bvecs';'*.*'},'Select the bvecs file...');
    if(isnumeric(f)), disp('User canceled.'); return; end
    bvecsFile = fullfile(p,f);
end

ecFile = [inBaseDir '_ecXform.mat'];
acpcFile = [inBaseDir '_acpcXform.mat'];
alignedBvecsFile = [outBaseDir '.bvecs'];
alignedBvalsFile = [outBaseDir '.bvals'];
dwAlignedRawFile = [outBaseDir '.nii.gz'];

disp('loading raw data...');
dwRaw = niftiRead(dwRawFileName);
oldPhaseDim = dwRaw.phase_dim;
% *** FIXME: the canXform might need to be applied to bvecs.
[dwRaw,canXform] = niftiApplyCannonicalXform(dwRaw);
% Make sure there is a valid phase-encode dir, as this is crucial for
% eddy-current correction.

if notDefined('phaseEncodeDir')
    while(dwRaw.phase_dim<1 || dwRaw.phase_dim>3)
        prompt = sprintf('Phase-encode dir is currently %d, but must be 1, 2 or 3. New value:',dwRaw.phase_dim);
        resp = inputdlg(prompt,'Set phase encode direction',1,{'2'});
        if(~isempty(resp))
            dwRaw.phase_dim = round(str2double(resp{1}));
        end
    end
else
    dwRaw.phase_dim=phaseEncodeDir;
end

doResamp = false;

bvecs = dlmread(bvecsFile);
bvals = dlmread(bvalsFile);
goodVols = squeeze(max(max(max(dwRaw.data))))~=0;
% Negative bvals are used to indicate bad volumes that should be skipped
if(any(bvals<0))
    goodVols = goodVols&bvals>0;
end
if(~isempty(excludeVols))
    goodVols(excludeVols) = false;
end
if(~all(goodVols))
  fprintf('Found %d bad volumes in data- removing them from analysis...\n',sum(~goodVols));
  dwRaw.data = dwRaw.data(:,:,:,goodVols);
  bvecs = bvecs(:,goodVols);
  bvals = bvals(goodVols);
  doResamp = true;
else
  if(length(goodVols)<size(bvecs,2))
	warning('mrDiffusion:dimMismatch', 'More bvecs than vols- ignoring some bvecs...');
    bvecs = bvecs(:,goodVols);
  end
  if(length(goodVols)<size(bvals,2))
	warning('mrDiffusion:dimMismatch', 'More bvals than vols- ignoring some bvals...');
    bvals = bvals(goodVols);
  end
end

% Get a mean b0 to be used for eddy-current correction and
% alignment to a structural scan.
if(clobber==1||~exist(mnB0Name,'file'))
  dtiRawComputeMeanB0(dwRaw, bvals, mnB0Name);
else
  if(clobber==0)
    resp = questdlg([mnB0Name ' exists- would you like to overwrite it?'], 'Clobber mnB0', 'Overwrite','Use Existing File','Abort','Use Existing File');
    if(strcmpi(resp,'Abort')), error('User aborted.'); end
    if(strcmpi(resp,'Overwrite'))
      dtiRawComputeMeanB0(dwRaw, bvals, mnB0Name);
    end
  end
end

% Compute the eddy-current correction for all the DWIs (SLOW!)
if(clobber==1||~exist(ecFile,'file'))
    dtiRawRohdeEstimateEddyMotion(dwRaw, mnB0Name, bvals, ecFile);
    doResamp = true;
else
    if(clobber==0)
        resp = questdlg([ecFile ' exists- would you like to overwrite it?'], 'Clobber EddyCorrect', ...
            'Overwrite','Use Existing File','Abort','Use Existing File');
        if(strcmpi(resp,'Abort')), error('User aborted.'); end
        if(strcmpi(resp,'Overwrite'))
            dtiRawRohdeEstimateEddyMotion(dwRaw, mnB0Name, bvals, ecFile);
            doResamp = true;
        end
    end
end

% Compute the dti-structural alignment
if(clobber==1||~exist(acpcFile,'file'))
  dtiRawAlignToT1(mnB0Name, t1FileName, acpcFile);
  doResamp = true;
else
  if(clobber==0)
    resp = questdlg([acpcFile ' exists- would you like to overwrite it?'], 'Clobber AcPc', ...
		    'Overwrite','Use Existing File','Abort','Use Existing File');
    if(strcmpi(resp,'Abort')), error('User aborted.'); end
    if(strcmpi(resp,'Overwrite'))
      dtiRawAlignToT1(mnB0Name, t1FileName, acpcFile);
      doResamp = true;
    end
  end
end

% Resample all the DWIs, applying the dti-to structural xform and
% the eddy-current correction xforms.
if((doResamp&&clobber==-1)||clobber==1||~exist(dwAlignedRawFile,'file'))
  dtiRawResample(dwRaw, ecFile, acpcFile, dwAlignedRawFile, bsplineInterpFlag, dwOutMm);
else
  if(clobber==0)
    resp = questdlg([dwAlignedRawFile ' exists- would you like to overwrite it?'], 'Clobber Resampled Data', ...
		    'Overwrite','Use Existing File','Abort','Use Existing File');
    if(strcmpi(resp,'Abort')), error('User aborted.'); end
    if(strcmpi(resp,'Overwrite'))
      dtiRawResample(dwRaw, ecFile, acpcFile, dwAlignedRawFile, bsplineInterpFlag, dwOutMm);
    end
  end
end

if((doResamp&&clobber==-1)||clobber==1||~exist(alignedBvecsFile,'file')||~exist(alignedBvalsFile,'file'))
    dtiRawReorientBvecs(bvecs, ecFile, acpcFile, alignedBvecsFile);
    %copyfile(bvalsFile, alignedBvalsFile);
	dlmwrite(alignedBvalsFile,bvals,' ');
else
    if(clobber==0)
        resp = questdlg([alignedBvecsFile ' exists- would you like to overwrite it?'], 'Clobber Bvecs/bvals', ...
            'Overwrite','Use Existing File','Abort','Use Existing File');
        if(strcmpi(resp,'Abort')), error('User aborted.'); end
        if(strcmpi(resp,'Overwrite'))
            dtiRawReorientBvecs(bvecsFile, ecFile, acpcFile, alignedBvecsFile);
            dlmwrite(alignedBvalsFile,bvals,' ');
        end
    end
end

clear dwRaw;
dwRawAligned = niftiRead(dwAlignedRawFile);

bs.n = numBootStrapSamples;
% We'll use the non-realigned bvecs since we want to count bvecs that are
% only a little differnt due to motion correction as 'repeats'. Also, we
% can count a direction with just a sign-flip as a 'repeat' since it will
% contain essentially the same diffusion info.
% Note that this code is now used just to compute the number of unique
% diffusion directions. We now use a residual bootstrap, so the repetion
% pattern is no long important for the bootstrap.
[bs.permuteMatrix, nNonRepeats, nUniqueDirs] = dtiBootGetPermMatrix(dlmread(bvecsFile), dlmread(bvalsFile));
% We still need an over-determined tensor fit to do residual bootstrap.
% We'll skip the bootstrap for datasets with fewer than 14 measurements
% (7 is the minimum for tensor fitting).
if(size(dwRawAligned.data,4)<14)
    warning('mrDiffusion:bootstrap','Not enough redundancy in the data- skipping bootstrap.');
    bs.n = 0;
end
bs.showProgress = false;
if(isempty(dt6BaseName))
    dt6BaseName = fullfile(subjectDir,sprintf('dti%02d',nUniqueDirs));
    if(bsplineInterpFlag)
        dt6BaseName = [dt6BaseName 'bs'];
    end
else
    if(isempty(fileparts(dt6BaseName)))
        % Don't allow a relative path
        dt6BaseName = fullfile(subjectDir,dt6BaseName);
    end
end

dt6FileName = dtiRawFitTensorMex(dwRawAligned, alignedBvecsFile, alignedBvalsFile, dt6BaseName, bs);
% Make sure paths are correct
%dtiRawFixDt6File(dt6FileName);

% snCheckIm = dtiComputeSpatialNorm(dt6FileName);
% imwrite(makeMontage(snCheckIm),fullfile(dt6BaseName,'MNIalignmentCheck.png'),'CreationTime',now,'Author','mrDiffusion from Stanford University','Description','b0 with MNI overlay');

dt6 = load(dt6FileName,'files');
files = dt6.files;
% Make the t1 filename relative
k = strfind(t1FileName,dtiGetSubjDirInDT6(dt6FileName));
if isempty(k) || ~exist(t1FileName(k:end),'file')
  [p,t1RelFileName,e] = fileparts(t1FileName);
  t1RelFileName = [t1RelFileName e];
  if(exist(fullfile(subjectDir,'t1'),'dir'))
	t1RelFileName = fullfile('t1',t1RelFileName);
  end
else
  t1RelFileName = t1FileName(k+1+length(dtiGetSubjDirInDT6(dt6FileName)):end);
end
files.t1 = t1RelFileName;
% add the raw data file names (these will be full paths)
files.alignedDwRaw = dwAlignedRawFile;
file.alignedDwBvecs = alignedBvecsFile;
file.alignedDwBvals = alignedBvalsFile;
save(dt6FileName,'files','-APPEND');

[pddT1,xform,mm] = dtiRawCheckTensors(fullfile(dt6BaseName,'bin','tensors.nii.gz'),t1FileName);
pddT1 = flipdim(permute(pddT1,[2 1 3 4]),1);
imSlices = [1:2:size(pddT1,3)];
m = makeMontage3(pddT1,imSlices,mm(1),0,[],[],0);
%imshow(m)
imwrite(m,fullfile(dt6BaseName,'t1pdd.png'),'CreationTime',now,'Author','mrDiffusion from Stanford University','Description','T1 with PDD overlay');

% Setup conTrack directory
fiberDir = fullfile(subjectDir,'fibers');
if(~exist(fiberDir,'dir')); mkdir(fiberDir); end;
conTrackDir = fullfile(fiberDir,'conTrack');
if(~exist(conTrackDir,'dir')); mkdir(conTrackDir); end;
paramsName = fullfile(conTrackDir,'met_params.txt');

% Create conTrack options file if it doesn't exist or we need to clobber
bWriteParams = 0;
if(clobber==1||~exist(paramsName,'file'))
    bWriteParams = 1;
else
    if(clobber==0&&exist(paramsName,'file'))
        resp = questdlg([paramsName ' exists- would you like to overwrite it?'], 'Clobber conTrack params', 'Overwrite','Use Existing File','Abort','Use Existing File');
        if(strcmpi(resp,'Abort')); error('User aborted.'); end
        if(strcmpi(resp,'Overwrite'))
            bWriteParams = 1;
        end
    end
end
if( bWriteParams )
    mtr = mtrCreate();
    mtr = mtrSet(mtr, 'tensors_filename', fullfile(subjectDir,'bin','tensors.nii.gz'));
    mtr = mtrSet(mtr, 'fa_filename', fullfile(subjectDir,'bin','wmMask.nii.gz'));
    mtr = mtrSet(mtr, 'pdf_filename', fullfile(subjectDir,'bin','pddDispersion.nii.gz'));
    mtrSave(mtr,paramsName);
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function plotMicArray(mic_dirs_deg, R)
% PLOTMICARRAY plots the arrangement of the microphones in a spherical array
%
%   mic_dirs:   the directions of the microphones in [azi1 elev1; ...]
%       convention
%   R:          radius of the array
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PLOTMICARRAY.M - 11/7/2013
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mic_dirs_rad = mic_dirs_deg*pi/180;
Nmic = size(mic_dirs_deg,1);

hold on
% set up unit sphere information
numSphereFaces = 20;
[unitSphereX, unitSphereY, unitSphereZ] = sphere(numSphereFaces);
% radius of each sphere
spheresRadius = ones(Nmic,1)*0.05*R;
%plot 3d axes
line([0 2*R],[0 0], [0 0],'color','r');
text(2*R,0,0,'x','Color','r','FontSize',24);
line([0 0],[0 2*R], [0 0],'color','g');
text(0,2*R,0,'y','Color','g','FontSize',24);
line([0 0],[0 0], [0 2*R],'color','b');
text(0,0,2*R,'z','Color','b','FontSize',24);

spheresX = R*cos(mic_dirs_rad(:,1)).*cos(mic_dirs_rad(:,2));
spheresY = R*sin(mic_dirs_rad(:,1)).*cos(mic_dirs_rad(:,2));
spheresZ = R*sin(mic_dirs_rad(:,2));
% for each given sphere, shift the scaled unit sphere by the
% location of the sphere and plot
for i=1:Nmic
    sphereX = spheresX(i) + unitSphereX*spheresRadius(i);
    sphereY = spheresY(i) + unitSphereY*spheresRadius(i);
    sphereZ = spheresZ(i) + unitSphereZ*spheresRadius(i);
    h = surf(sphereX, sphereY, sphereZ);
    set(h, 'FaceColor','b', 'FaceLighting', 'gouraud')

    text(1.1*spheresX(i),1.1*spheresY(i),1.1*spheresZ(i), num2str(i),'color','w','FontSize',24)
end
sphereX = unitSphereX*R;
sphereY = unitSphereY*R;
sphereZ = unitSphereZ*R;
h = surf(sphereX, sphereY, sphereZ);
set(h, 'FaceColor','c', 'FaceLighting', 'gouraud')
light('Position',[0 0 1],'Style','infinite');
material shiny
hold off
axis equal
grid on
view(-30,15)
set(gca,'visible','off')

author: polarch

repo: Spherical-Array-Processing


text:
%% RATE OF CONVERGENCE OF LINEAR ELEMENT FOR POISSON EQUATION
%
% This example is to show the rate of convergence of linear finite element
% approximation of the Poisson equation on the unit square with the
% following boundary conditions:
%
% - Non-empty Dirichlet boundary condition.
% - Pure Neumann boundary condition.
% - Robin boundary condition.
%
% The basis, data structure and numerical test is summarized in <a
% href="matlab:ifem Poissonfemrate">Poissonfemrate</a>.
%
% See also PoissonP2femrate, Poissonafemrate
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

clear variables
%% Setting
[node,elem] = squaremesh([0,1,0,1],0.25);
mesh = struct('node',node,'elem',elem);
option.L0 = 3;
option.maxIt = 4;
option.printlevel = 1;

%% Non-empty Dirichlet boundary condition.
option.plotflag = 1;
pde = sincosdata;
mesh.bdFlag = setboundary(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
femPoisson(mesh,pde,option);

%% Pure Neumann boundary condition.
option.plotflag = 0;
% pde = sincosNeumanndata;
pde = sincosdata;
mesh.bdFlag = setboundary(node,elem,'Neumann');
femPoisson(mesh,pde,option);

%% Pure Robin boundary condition.
option.plotflag = 0;
pde = sincosRobindata;
mesh.bdFlag = setboundary(node,elem,'Robin');
femPoisson(mesh,pde,option);

%% Conclusion
%
% The optimal rate of convergence of the H1-norm (1st order) and L2-norm
% (2nd order) is observed. The 2nd order convergent rate between two
% discrete functions ||DuI-Duh|| is known as superconvergence.
%
% MGCG converges uniformly in all cases.

repo: ifem


text:
function crcCheck = checkCRC(response)
receivedMsg = response;
challengeMsg = append_crc(response(1:length(response)-2));
if all(challengeMsg==receivedMsg)
    crcCheck = 1;
else
    crcCheck = 0;
end


text:
function f = sinh( f )
%SINH   Hyperbolic sine of a SEPARABLEAPPROX.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty( f ) )
    return
end

f = compose( f, @sinh );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function r8vec_index_delete_dupes_test ( )

%*****************************************************************************80
%
%% R8VEC_INDEX_DELETE_DUPES_TEST tests R8VEC_INDEX_DELETE_DUPES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  n_max = 25;
  n = 0;
  x = [];
  indx = [];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_INDEX_DELETE_DUPES_TEST\n' );
  fprintf ( 1, '  R8VEC_INDEX_DELETE_DUPES deletes duplicates.\n' );
  fprintf ( 1, '\n' );

  xval = 8.0;
  [ n, x, indx ] = r8vec_index_insert ( n, x, indx, xval );

  xval = 7.0;
  [ n, x, indx ] = r8vec_index_insert ( n, x, indx, xval );

  seed = 123456789;

  for i = 1 : 20
    [ xval, seed ] = r8_uniform_ab ( 0.0, 20.0, seed );
    xval = round ( xval );
    fprintf ( 1, '  %f\n', xval );
    [ n, x, indx ] = r8vec_index_insert ( n, x, indx, xval );
  end

  xval = 7.0;
  [ n, x, indx ] = r8vec_index_insert ( n, x, indx, xval );

  xval = 8.0;
  [ n, x, indx ] = r8vec_index_insert ( n, x, indx, xval );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Indexed list of entries:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  INDX(I)  X(I)  X(INDX(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %3d  %3d  %9f  %9f\n', i, indx(i), x(i), x(indx(i)) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Call R8VEC_INDEX_DELETE_DUPES to delete duplicates:\n' );

  [ n, x, indx ] = r8vec_index_delete_dupes ( n, x, indx );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Indexed list of unique entries:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  INDX(I)  X(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %3d  %3d  %9f\n', i, indx(i), x(i) );
  end

  return
end