text:
function ConvertPhyKilo2Neurosuite(basepath,basename)

% Converts Phy kilosort output files into klusters-compatible
% fet,res,clu,spk files.  Works on a single shank of a recording, assumes a
% 16bit .dat, .xml file is present in "basepath" (home folder) and
% that they are named basename.dat and basename.xml.
% It assumes that your directory contains:
% spike_clusters.npy, spikes_times.npy, cluster_group.tsv, pcfeatures.npy,
% templates.npy, cluster_ids.npy and shanks.npy (these last two are
% generated by a phy plugin -Export Shanks- found in the lab wiki)
%
% Inputs:
%   basepath -  directory path to the main recording folder with .dat and .xml
%               as well as shank folders made by makeProbeMapKlusta2.m (default is
%               current directory matlab is pointed to)
%   basename -  shared file name of .dat and .xml (default is last part of
%               current directory path, ie most immediate folder name)
%
% Eliezyer de Oliveira, 2018

if ~exist('basepath','var')
    [~,basename] = fileparts(cd);
    basepath = cd;
end
savepath = basepath;

if ~exist('rez','var');
    load(fullfile(basepath,'rez.mat'))
end

Nchan = rez.ops.Nchan;
connected    = rez.connected;
xcoords      = rez.xc;
ycoords      = rez.yc;
% Nchan = rez.ops.Nchan;
% connected    = ones(Nchan, 1);
% xcoords      = ones(Nchan, 1);
% ycoords      = (1:Nchan)';

par = LoadParameters(fullfile([basepath],[basename '.xml']));

totalch = par.nChannels;
sbefore = 16;%samples before/after for spike extraction
safter = 16;%... could read from SpkGroups in xml
if isfield(par,'SpkGrps')
    if isfield(par.SpkGrps,'nSamples')
        if ~isempty(par.SpkGrps(1).nSamples);
            if isfield(par.SpkGrps,'PeakSample')
                if ~isempty(par.SpkGrps(1).PeakSample);
                    sbefore = par.SpkGrps(1).PeakSample;
                    safter = par.SpkGrps(1).nSamples - par.SpkGrps(1).PeakSample;
                end
            end
        end
    end
end

if exist(rez.ops.fbinary,'file')
    datpath = rez.ops.fbinary;
else
    datpath = fullfile(basepath,[basename '.dat']);
end

%% identify timestamps that are from good clusters
spikeTimes = readNPY('spike_times.npy');
clusters = readNPY('spike_clusters.npy');
S = tdfread('cluster_group.tsv');
group = S.group;
cluster_id = S.cluster_id;
c = 0;
for idx = 1:size(group,1)
    if group(idx,1:4) == 'good'
        c  = c+1;
        GClusters(c) = idx;
    end
end
ExtClus = cluster_id(GClusters);

% Separating idx by cluster
auxiliarC = [];
for idx = 1:length(ExtClus)
    auxC = find(clusters == ExtClus(idx));
    auxiliarC = cat(1,auxiliarC,auxC);
end
[~,I] = sort(auxiliarC);
auxiliarC = auxiliarC(I);
%% getting spike information

spktimes = uint64(readNPY('spike_times.npy'));
spktimes = spktimes(auxiliarC);
clu = uint32(readNPY('spike_clusters.npy'));
clu = clu(auxiliarC);
pcFeatures = readNPY('pc_features.npy');
pcFeatures = pcFeatures(auxiliarC,:,:);
% pcFeatureInds = uint32(readNPY('pc_feature_ind.npy'))';
templates = readNPY('templates.npy');

mkdir(fullfile(savepath,'PhyClus'))
%% assigning cluster ids to shanks
shanks = readNPY('shanks.npy');
cluShank = readNPY('cluster_ids.npy');

auxC = unique(clu);
templateshankassignments = zeros(size(auxC));
for idx = 1:length(auxC)
    temp = find(cluShank == auxC(idx));
    templateshankassignments(idx) = shanks(temp);
end
grouplookup = rez.ops.kcoords;
allgroups = unique(grouplookup);
allgroups(allgroups==0) = [];

for groupidx = 1:length(allgroups)

    %if isfield(par.SpkGrps(groupidx),'Channels')
    %if ~isempty(par.SpkGrps(groupidx).Channels)
    % for each group loop through, find all templates clus
    tgroup          = allgroups(groupidx);%shank number
    ttemplateidxs   = find(templateshankassignments==tgroup);%which templates/clusters are in that shank
    %     ttemplates      = templates(:,:,ttemplateidxs);
    %     tPCFeatureInds  = pcFeatureInds(:,ttemplateidxs);

    tidx            = ismember(clu,auxC(ttemplateidxs));%find spikes indices in this shank
    tclu            = clu(tidx);%extract template/cluster assignments of spikes on this shank
    tspktimes       = spktimes(tidx);

    gidx            = find(rez.ops.kcoords == tgroup);%find all channels in this group
    channellist     = [];

    for ch = 1:length(par.spikeGroups.groups{groupidx})
        if sum(ismember(gidx,par.spikeGroups.groups{groupidx}(:)+1))
            channellist = par.spikeGroups.groups{groupidx}(:)+1;
            break
        end
    end
    if isempty(channellist)
        disp(['Cannot find spkgroup for group ' num2str(groupidx) ])
        continue
    end

    %% spike extraction from dat
    if groupidx == 1;
        dat             = memmapfile(datpath,'Format','int16');
    end
    tsampsperwave   = (sbefore+safter);
    ngroupchans     = length(channellist);
    valsperwave     = tsampsperwave * ngroupchans;
    wvforms_all     = zeros(length(tspktimes)*tsampsperwave*ngroupchans,1,'int16');
    wvranges        = zeros(length(tspktimes),ngroupchans);
    wvpowers        = zeros(1,length(tspktimes));

    for j=1:length(tspktimes)
        try
            w = dat.data((double(tspktimes(j))-sbefore).*totalch+1:(double(tspktimes(j))+safter).*totalch);
            wvforms=reshape(w,totalch,[]);
            %select needed channels
            wvforms = wvforms(channellist,:);
            %         % detrend
            %         wvforms = floor(detrend(double(wvforms)));
            % median subtract
            wvforms = wvforms - repmat(median(wvforms')',1,sbefore+safter);
            wvforms = wvforms(:);

        catch
            disp(['Error extracting spike at sample ' int2str(double(tspktimes(j))) '. Saving as zeros']);
            disp(['Time range of that spike was: ' num2str(double(tspktimes(j))-sbefore) ' to ' num2str(double(tspktimes(j))+safter) ' samples'])
            wvforms = zeros(valsperwave,1);
        end

        %some processing for fet file
        wvaswv = reshape(wvforms,tsampsperwave,ngroupchans);
        wvranges(j,:) = range(wvaswv);
        wvpowers(j) = sum(sum(wvaswv.^2));

        lastpoint = tsampsperwave*ngroupchans*(j-1);
        wvforms_all(lastpoint+1 : lastpoint+valsperwave) = wvforms;
        %     wvforms_all(j,:,:)=int16(floor(detrend(double(wvforms)')));
        if rem(j,100000) == 0
            disp([num2str(j) ' out of ' num2str(length(tspktimes)) ' done'])
        end
    end
    wvranges = wvranges';

    %% Spike features
    %     for each template, rearrange the channels to reflect the shank order
    tdx = [];
    for tn = 1:size(pcFeatures,3)
%         %         tTempPCOrder = tPCFeatureInds(:,tn);%channel sequence used for pc storage for this template
        for k = 1:length(channellist);
%             %             i = find(tTempPCOrder==channellist(k));
            if ~isempty(k)
                tdx(tn,k) = k;
            else
                tdx(tn,k) = nan;
            end
%
        end
    end

    featuresperspike = 3; % kilosort default

    % initialize fet file
    fets    = zeros(sum(tidx),size(pcFeatures,2),size(pcFeatures,3));
    pct     = pcFeatures(tidx,:,:);

    %for each cluster/template id, grab at once all spikes in that group
    %and rearrange their features to match the shank order
    allshankclu = unique(tclu);

    for tc = 1:length(allshankclu)
        tsc     = allshankclu(tc);
        i       = find(tclu==tsc);
        tforig  = pct(i,:,:);%the subset of spikes with this clu ide
        tfnew   = tforig; %will overwrite

        ii      = tdx(:,16);%handling nan cases where the template channel used was not in the shank
        gixs    = ~isnan(ii);%good vs bad channels... those shank channels that were vs were not found in template pc channels
        bixs    = isnan(ii);
        g       = ii(gixs);

        tfnew(:,:,gixs) = tforig(:,:,gixs);%replace ok elements
        tfnew(:,:,bixs) = 0;%zero out channels that are not on this shank
        try
            fets(i,:,:) = tfnew(:,:,1:size(pcFeatures,3));
        catch
            keyboard
        end
    end
    %extract for relevant spikes only...
    % and heurstically on d3 only take fets for one channel for each original channel in shank... even though kilosort pulls 12 channels of fet data regardless
    tfet1 = squeeze(fets(:,1,1:size(pct,3)));%lazy reshaping
    tfet2 = squeeze(fets(:,2,1:size(pct,3)));
    tfet3 = squeeze(fets(:,3,1:size(pct,3)));
    fets = cat(2,tfet1,tfet2,tfet3)';%     fets = h5read(tkwx,['/channel_groups/' num2str(shank) '/features_masks']);
    %     fets = double(squeeze(fets(1,:,:)));
    %mean activity per spike
    %     fetmeans = mean(fets,1);%this is pretty redundant with wvpowers
    %     %find first pcs, make means of those...
    %     featuresperspike = 4;
    %     firstpcslist = 1:featuresperspike:size(fets,1);
    %     firstpcmeans = mean(fets(firstpcslist,:),1);
    %
    %     nfets = size(fets,1)+1;
    %     fets = cat(1,fets,fetmeans,firstpcmeans,wvpowers,wvranges,double(tspktimes'));
    fets = cat(1,double(fets),double(wvpowers),double(wvranges),double(tspktimes'));
    fets = fets';
    % fets = cat(1,nfets,fets);

    %% writing to clu, res, fet, spk

    cluname = fullfile(savepath, [basename '.clu.' num2str(tgroup)]);
    resname = fullfile(savepath, [basename '.res.' num2str(tgroup)]);
    fetname = fullfile(savepath, [basename '.fet.' num2str(tgroup)]);
    spkname = fullfile(savepath, [basename '.spk.' num2str(tgroup)]);
    %fet
    SaveFetIn(fetname,fets);

    %clu
    % if ~exist(cluname,'file')
    tclu = cat(1,length(unique(tclu)),double(tclu));
    fid=fopen(cluname,'w');
    %     fprintf(fid,'%d\n',clu);
    fprintf(fid,'%.0f\n',tclu);
    fclose(fid);
    clear fid
    % end

    %res
    fid=fopen(resname,'w');
    fprintf(fid,'%.0f\n',tspktimes);
    fclose(fid);
    clear fid

    %spk
    fid=fopen(spkname,'w');
    fwrite(fid,wvforms_all,'int16');
    fclose(fid);
    clear fid

    disp(['Shank ' num2str(tgroup) ' done'])
    %end
    %end
end
clear dat
copyfile(fullfile(savepath, [basename,'.clu.*']),fullfile(savepath, 'PhyClus'))

function SaveFetIn(FileName, Fet, BufSize);

if nargin<3 | isempty(BufSize)
    BufSize = inf;
end

nFeatures = size(Fet, 2);
formatstring = '%d';
for ii=2:nFeatures
    formatstring = [formatstring,'\t%d'];
end
formatstring = [formatstring,'\n'];

outputfile = fopen(FileName,'w');
fprintf(outputfile, '%d\n', nFeatures);

if isinf(BufSize)

    temp = [round(100* Fet(:,1:end-1)) round(Fet(:,end))];
    fprintf(outputfile,formatstring,temp');
else
    nBuf = floor(size(Fet,1)/BufSize)+1;

    for i=1:nBuf
        BufInd = [(i-1)*nBuf+1:min(i*nBuf,size(Fet,1))];
        temp = [round(100* Fet(BufInd,1:end-1)) round(Fet(BufInd,end))];
        fprintf(outputfile,formatstring,temp');
    end
end
fclose(outputfile);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function p = vonMises_prob( x, m, k, use_log )
%VONMIS_PROB Calculates the probability of x coming from a Von Mises
%distribution with mean mu and concentration parameter k.
%   p = vonMises_prob( x, m, k, use_log )

if nargin < 4, use_log = 0; end
[d N] = size(x);

m = m(:);
M = m*ones(1,N);

denom = (2*pi)*besseli(0,k);
mahal = (k*cos(x-M));

if use_log
    p = mahal-log(denom);
else
    p = exp(mahal)/denom;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function y = cvx_isaffine( x, full )
error( nargchk( 1, 2, nargin ) );
y = cvx_vexity( x );
if nargin < 2,
    y = nnz( y ) == 0;
else
    y = cvx_reshape( y == 0, x.size_ );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
classdef CEC2010_F4 < PROBLEM
% <single> <real> <constrained>
% CEC'2010 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% R. Mallipeddi and P. N. Suganthan, Problem definitions and evaluation
% criteria for the CEC 2010 competition on constrained real-parameter
% optimization, Nanyang Technological University, Singapore, 2010.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;  % Optimal decision vector
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2010.mat'),'Data');
            obj.O = Data{4};
            obj.M = 1;
            if isempty(obj.D); obj.D = 10; end
            obj.D = min(obj.D,length(obj.O));
            obj.lower    = zeros(1,obj.D) - 50;
            obj.upper    = zeros(1,obj.D) + 50;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = max(Z,[],2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopCon(:,1) = abs(mean((Z.*cos(sqrt(abs(Z)))),2)) - 1e-4;
            PopCon(:,2) = abs(sum((Z(:,1:end/2-1)-Z(:,2:end/2)).^2,2)) - 1e-4;
            PopCon(:,3) = abs(sum((Z(:,end/2+1:end-1).^2-Z(:,end/2+2:end)).^2,2)) - 1e-4;
            PopCon(:,4) = abs(sum(Z,2)) - 1e-4;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function kern = polyardKernExpandParam(kern, params)

% POLYARDKERNEXPANDPARAM Create kernel structure from POLYARD kernel's parameters.
% FORMAT
% DESC returns a automatic relevance determination polynomial kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : polyardKernParamInit, polyardKernExtractParam, kernExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% KERN

kern.weightVariance = params(1);
kern.biasVariance = params(2);
kern.variance = params(3);
kern.inputScales = params(4:end);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function imchi=kkimbook2(omega,rechi,alpha)
%The program inputs are the vector of the frequency (or energy)
%components, the vector of the real part of the susceptibility
%under examination, and the value of the moment considered.
%The two vectors must have the same length
%and the frequency vector omega must be equispaced.
%If not, apply MATLAB functions such as interp.
%If rechi is the real part of a linear susceptibility,
%alpha must be 0.
%If rechi is the real part of the nth
%harmonic generation susceptibility, alpha=0,1,..2n.
%If rechi is the real part of a pump and probe
%susceptibility, alpha=0 or 1.
%This files accompanies the book
%"Kramers-Kronig Relations in Optical Materials Research"
%by Lucarini, V., Saarinen, J.J., Peiponen, K.-E., Vartiainen, E.M.
%Springer, Heidelberg, 2005
%where the theory and applications are fully developed.
%The output is the estimate of the imaginary part as obtained
%with K-K relations.
%This software is distributed under the GNU licence agreement
%by Valerio Lucarini
%email: lucarini@alum.mit.edu
%University of Camerino
%Department of Mathematics and Computer Science
%Camerino, Italy

if size(omega,1)>size(omega,2);
omega=omega';
end; if size(rechi,1)>size(rechi,2);
rechi=rechi';
end;
%Here the program rearranges the two vectors so that,
%whichever their initial shape, they become row vectors.

g=size(omega,2);
%Size of the vectors.%

imchi=zeros(size(rechi));
%The output is initialized.

a=zeros(size(rechi));
b=zeros(size(rechi));
%Two vectors for intermediate calculations are initialized

deltaomega=omega(2)-omega(1);
%Here we compute the frequency (or energy) interval

j=1;
beta1=0;
for k=2:g;
b(1)=beta1+rechi(k)*omega(k)^(2*alpha)/(omega(k)^2-omega(1)^2);
beta1=b(1);
end;
imchi(1)=-2/pi*deltaomega*b(1)*omega(1)^(1-2*alpha);
%First element of the output: the principal part integration
%is computed by excluding the first element of the input

j=g;
alpha1=0;
for k=1:g-1;
a(g)=alpha1+rechi(k)*omega(k)^(2*alpha)/(omega(k)^2-omega(g)^2);
alpha1=a(g);
end;
imchi(g)=-2/pi*deltaomega*a(g)*omega(g)^(1-2*alpha);
%Last element of the output: the principal part integration
%is computed by excluding the last element of the input.

for j=2:g-1; ;
%Loop on the inner components of the output vector.
alpha1=0;
beta1=0;
for k=1:j-1;
a(j)=alpha1+rechi(k)*omega(k)^(2*alpha)/(omega(k)^2-omega(j)^2);
alpha1=a(j);
end;
for k=j+1:g;
b(j)=beta1+rechi(k)*omega(k)^(2*alpha)/(omega(k)^2-omega(j)^2);
beta1=b(j);
end;
imchi(j)=-2/pi*deltaomega*(a(j)+b(j))*omega(j)^(1-2*alpha);
%Last element of the output: the principal part integration
%is computed by excluding the last element of the input
end;

author: Sable

repo: mcbench-benchmarks


text:
function [overall_peaks, overall_onsets] = co_ppg_peak_detector(ppg, fs, up)
%CO_PPG_PEAK_DETECTOR detects PPG pulse peaks
%
% This has been adapted by Peter Charlton from code supplied by Christina
% Orphanidou.
%
% Purpose: To detect PPG peaks.
%
% With grateful thanks to Christina Orphanidou and Alexander Darrell.
%

%% Pre-processing
% remove nans
ppg.v(isnan(ppg.v)) = mean(ppg.v(~isnan(ppg.v)));
% filter
ppgfilt.v = filtfilt(up.paramSet.CO_peak_det.PPG_smoother.coeffs, 1, ppg.v);
ppgfilt.t = ppg.t;

%% Setup
% Segmentation into 10s windows
win_length = 10;
overlap = 3;
win_starts = ppgfilt.t(1):win_length:ppgfilt.t(end);
win_ends = win_starts + win_length + overlap;
win_ends(end) = win_ends(end) - overlap;
% setup variables
[overall_peaks, overall_onsets] = deal([]);

%% Perform peak detection on each 10s window
for win_no = 1 : length(win_starts)

    rel_data = struct;

    rel_data.i = find(ppgfilt.t >= win_starts(win_no) & ppgfilt.t <= win_ends(win_no));
    rel_data.t = ppgfilt.t(rel_data.i);
    rel_data.v = ppgfilt.v(rel_data.i);

    %% Calculate thresholds
    thresh1=quantile(rel_data.v,up.paramSet.CO_peak_det.upctl);
    thresh2=quantile(rel_data.v,up.paramSet.CO_peak_det.lpctl);
    thresh3=thresh2+0.3*(thresh1-thresh2);
    thresh4=thresh2+0.7*(thresh1-thresh2);

    %% Find all peaks
    % identify peaks
    diffs_on_left_of_pt = diff(rel_data.v); diffs_on_left_of_pt = diffs_on_left_of_pt(1:(end-1)); diffs_on_left_of_pt = logical(diffs_on_left_of_pt>0);
    diffs_on_right_of_pt = diff(rel_data.v); diffs_on_right_of_pt = diffs_on_right_of_pt(2:end); level_on_right_of_pt = logical(diffs_on_right_of_pt==0); diffs_on_right_of_pt = logical(diffs_on_right_of_pt<0);
    diffs_on_second_right_of_pt = [1; diff(rel_data.v(2:end))]; diffs_on_second_right_of_pt = diffs_on_second_right_of_pt(2:end); diffs_on_second_right_of_pt = logical(diffs_on_second_right_of_pt<0);
    peaks.i = find((diffs_on_left_of_pt & diffs_on_right_of_pt) | ...
        (diffs_on_left_of_pt & level_on_right_of_pt & diffs_on_second_right_of_pt))+1;
    peaks.i = peaks.i';
    % Take maximum value in window if no peaks were found
    if isempty(peaks.i)
        peaks.i = find(max(rel_data.v) == rel_data.v);
    end
    % Extract signal values at peaks
    peaks.t = rel_data.t(peaks.i);
    peaks.v = rel_data.v(peaks.i);

    %% Classify peaks
    % Identify relevant peaks according to amplitude
    upperdiff = abs(peaks.v-thresh1);
    middlehighdiff = abs(peaks.v-thresh4);
    middlelowdiff = abs(peaks.v-thresh3);
    lowerdiff = abs(peaks.v-thresh2);
    upper_pks = find(upperdiff<middlehighdiff & upperdiff<middlelowdiff & upperdiff<lowerdiff);
    PPG_PKS.i = peaks.i(upper_pks);
    PPG_PKS.v = peaks.v(upper_pks);
    PPG_PKS.t = peaks.t(upper_pks);
    % eliminate peaks which are too close together in time
    good_els = find(diff(PPG_PKS.i)>=fs/3)+1;
    PPG_PKS.i = PPG_PKS.i(good_els);
    PPG_PKS.v = PPG_PKS.v(good_els);
    PPG_PKS.t = PPG_PKS.t(good_els);

    %% Find troughs

    PPG_TRS.i = nan(length(PPG_PKS.t)-1,1);
    for s = 1 : (length(PPG_PKS.t)-1)
        start_el = PPG_PKS.i(s);
        [~, additional_el] = min(rel_data.v(PPG_PKS.i(s):PPG_PKS.i(s+1)));
        PPG_TRS.i(s) = start_el - 1 + additional_el;
    end
    PPG_TRS.v = rel_data.v(PPG_TRS.i);
    PPG_TRS.t = rel_data.t(PPG_TRS.i);

    temp_pk_is = rel_data.i(PPG_PKS.i(:));
    temp_on_is = rel_data.i(PPG_TRS.i(:));

    overall_peaks = [overall_peaks; temp_pk_is(:)]; overall_peaks = unique(overall_peaks);
    overall_onsets = [overall_onsets; temp_on_is(:)]; overall_onsets = unique(overall_onsets);

    clear upper lower middlehigh middlelow temp_pk_is temp_on_is PPG_TRS PPG_PKS additional_el s good_els peaks lower_pks upper_pks lowerdiff middlelowdiff middlehighdiff upperdiff ind

end

author: peterhcharlton

repo: RRest


text:
function SPGR_PlotSimCurve(MTdata, MTnoise, Prot, Sim, SimCurveResults)

leg = {};
if ~isempty(MTdata)
    semilogx(Prot.Offsets, MTdata,'kx', 'MarkerSize', 8); hold on;
    leg{1} = 'Raw data';
end
if (Sim.Opt.AddNoise)
    semilogx(Prot.Offsets, MTnoise,'bo','MarkerSize',8);
end

semilogx(SimCurveResults.Offsets, SimCurveResults.curve, 'LineWidth', 2);

if (Sim.Opt.AddNoise)
   leg{2} = 'Noisy data';
end

nAngles  = size(SimCurveResults.curve,2);
nOffsets = length(Prot.Offsets)/nAngles;
for i = 1:nAngles
    if Prot.Angles(1) == Prot.Angles(2)
        leg{end+1} = sprintf('Fitted curve (angle = %0.0f)',Prot.Angles(nOffsets*i));
    else
        leg{end+1} = sprintf('Fitted curve (angle = %0.0f)',Prot.Angles(i));
    end
end

if ~moxunit_util_platform_is_octave
    hleg = legend(leg, 'FontSize', 12);
    legend('boxoff')
    set(hleg,'Location', 'best')
else
    hleg = legend(leg);
    set(hleg,'FontSize',8)
    set(hleg,'location', 'southeast')
end
xlabel('Offset (Hz)','FontWeight','bold','FontSize',10);
ylabel('|Mz|','FontWeight','bold','FontSize',10);

set(gca, 'XScale', 'log');

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%   CLASS IRNSS_SS
% =========================================================================
%
% DESCRIPTION
%   container of IRNSS Satellite System parameters
%
% REFERENCES
%   CRS parameters, according to each GNSS system CRS definition
%   (ICD document in brackets):
%
%   *_QZS --> GRS80    (IS-IRNSS 1.8E)
%   Standard: https://www.google.it/url?sa=t&rct=j&q=&esrc=s&source=web&cd=2&ved=0ahUKEwjF5LXHyOPWAhUJLZoKHenOD4MQFggvMAE&url=https%3A%2F%2Fforum.nasaspaceflight.com%2Findex.php%3Faction%3Ddlattach%3Btopic%3D36710.0%3Battach%3D634319&usg=AOvVaw0sOE3-CqbC2o1p2tjjO8IM
%
%   Other useful links
%     - http://www.navipedia.net/index.php/IRNSS_Signal_Plan
%     - Ellipsoid: http://www.unoosa.org/pdf/icg/2012/template/WGS_84.pdf
%       note that GM and OMEGAE_DOT are redefined in the standard IS-GPS200H (GPS is not using WGS_84 values)
%     - http://www.navipedia.net/index.php/Reference_Frames_in_GNSS
%     - http://gage6.upc.es/eknot/Professional_Training/PDF/Reference_Systems.pdf

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Andrea Gatti
%  Contributors:      Andrea Gatti, Giulio Tagliaferro ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

classdef IRNSS_SS < Satellite_System
    properties (Constant, Access = 'public')
        SYS_EXT_NAME = 'IRNSS';    % full name of the constellation
        SYS_NAME     = 'IRN';     % 3 characters name of the constellation, this "short name" is used as fields of the property list (struct) to identify a constellation
        SYS_C        = 'I';       % Satellite system (ss) character id

        % System frequencies as struct [MHz]
        F = struct('L5', 1176.450, ...
                   'S',  2492.028)

        % Array of supported frequencies [MHz]
        F_VEC = struct2array(IRNSS_SS.F) * 1e6;

        % Array of the corresponding wavelength - lambda => wavelengths
        L_VEC = 299792458 ./ IRNSS_SS.F_VEC;

        N_SAT = 9;                    % Maximum number of satellite in the constellation
        PRN = (1 : IRNSS_SS.N_SAT)';  % Satellites id numbers as defined in the constellation

        % CODE2DATA ftp://igs.org/pub/data/format/rinex303.pdf
        CODE_RIN3_ATTRIB  = {'XCBA F', 'XCBA F'}; % last letter of the observation code e.g. C5A - C5B - C5C - C5X
        CODE_RIN3_DEFAULT_ATTRIB  = {'C' 'C'}; % last letter of the observation code
        CODE_RIN3_2BAND  = '59';             % id for the freq as stored in F_VEC e.g. L5 -> C5A, S -> C9A
        IONO_FREE_PREF  = ['59'];  % to be evaluated which combination is really better
end

    properties (Constant, Access = 'private')
        % IRNSS (GRS80) Ellipsoid semi-major axis [m]
        ELL_A = 6378137;
        % IRNSS (GRS80) Ellipsoid flattening
        ELL_F = 1/298.257222101;
        % IRNSS (GRS80) Ellipsoid Eccentricity^2
        ELL_E2 = (1 - (1 - IRNSS_SS.ELL_F) ^ 2);
        % IRNSS (GRS80) Ellipsoid Eccentricity
        ELL_E = sqrt(IRNSS_SS.ELL_E2);
    end

    properties (Constant, Access = 'public')
        % Structure of orbital parameters (ellipsoid, GM, OMEGA_EARTH_DOT)
        ORBITAL_P = struct('GM', 3.986005e14, ...                  % Gravitational constant * (mass of Earth) [m^3/s^2]
                           'OMEGAE_DOT', 7.2921151467e-5, ...      % Angular velocity of the Earth rotation [rad/s]
                           'ELL',struct( ...                       % Ellipsoidal parameters IRNSS (GRS80)
                           'A', IRNSS_SS.ELL_A, ...             % Ellipsoid semi-major axis [m]
                           'F', IRNSS_SS.ELL_F, ...             % Ellipsoid flattening
                           'E', IRNSS_SS.ELL_E, ...             % Eccentricity
                           'E2', IRNSS_SS.ELL_E2));             % Eccentricity^2
        ORBITAL_INC = 31;    % Orbital inclination
        ORBITAL_RADIUS  = 35786000 + 6378137; % Orbital radius
    end

    methods
        function this = IRNSS_SS(offset)
            % Creator
            % SYNTAX: IRNSS_SS(<offset>)
            if (nargin == 0)
                offset = 0;
            end
            this@Satellite_System(offset);
        end

        function copy = getCopy(this)
            % Get a copy of this
            copy = IRNSS_SS(this.getOffset());
            copy.import(this);
        end
    end
end

author: goGPS-Project


text:
function [dMu,dMu_n] = log_gmm_div_mu(X,W,Mu,Sigma)
%LOG_GMM_DIV_MU Derivative of the Log liklihood Guassian Mixture Model
%
%
%   input -----------------------------------------------------------------
%
%       o X         : (D x N),      Samples
%
%       o Priors    : (1 x K),      Weights
%
%       o Mu        : (D x K),      Means
%
%       o Sigma     : (D x D x K),  Covariance
%
%
%   output ----------------------------------------------------------------
%
%       o dMu       : (D x N x K), gradients of means for each data point
%                                  and each Gaussian Component
%

K = size(Mu,2);
D = size(Mu,1);
N = size(X,2);

dMu = zeros(D,K);

if nargout > 1
   dMu_n = zeros(D,N,K);
end

for k=1:K
    if nargout > 1
         dMu_n(:,:,k) = repmat(W(:,k)',D,1) .* (inv(Sigma(:,:,k)) * (X - repmat(Mu(:,k),1,N)));
         dMu(:,k)     = mean(dMu_n(:,:,k),2);
    else
         dMu(:,k) = mean(repmat(W(:,k)',D,1) .* (inv(Sigma(:,:,k)) * (X - repmat(Mu(:,k),1,N))),2);
    end
end

end

author: epfl-lasa


text:
%%*********************************************************************
%% This is a test example in using interface
%% to solve the following FAP problem
%% max <((k-1)/2k)*L(G,W)-1/2*Diag(We),X>
%% s.t. diag(X)==1 ; X positive semidefinite
%%      X(i,j) == -1/(k-1)  for (i,j) in U_set
%%      X(i,j) >= -1/(k-1)  for (i,j) in Edge_set\U_set
%% SDPNAL+:
%% Copyright (c) 2017 by
%% Yancheng Yuan, Kim-Chuan Toh, Defeng Sun and Xinyuan Zhao
%%*********************************************************************

%% read FAP data: U_set, Edge_set
 clear all;
 fname = 'fap08';
 if exist(fname)
    fid = fopen(fname,'r');
 elseif exist([fname,'.dat']);
    fid = fopen([fname,'.dat'],'r');
 else
    error('** Problem not found. \n');
 end
 [tmpr,count] = fscanf(fid,'%c');
 datavec = sscanf(tmpr,'%f'); clear tmpr;
 n = datavec(1);
 numedges = datavec(2);
 kpara = datavec(3);
 datavec = datavec(4:length(datavec));
 len = length(datavec);
 if (len ~= 3*numedges)
    error(' fapread: numedges and data do not match.');
 end
 I = datavec(1:3:len);
 J = datavec(2:3:len);
 w = datavec(3:3:len);
 idxU = find(w==1000);
 IU = I(idxU); JU = J(idxU); wU = w(idxU);
 idxE = find(w~=1000);
 IE = I(idxE); JE = J(idxE); wE = w(idxE);
 fclose(fid);
 %% model the problem using the interface
 GE = spconvert([IE JE wE; n n 0]);
 GE = GE + GE';
 GU = spconvert([IU JU wU; n n 0]);
 GU = GU + GU';
 LG = diag(GE*ones(n,1))-GE;
 C =  -0.5*diag(GE*ones(n,1)) + ((kpara-1)/(2*kpara))*LG;
 %%*********************************************************************
 model = ccp_model('Example_FAP');
     X = var_sdp(n, n);
     model.add_variable(X);
     model.maximize(inprod(C,X));
     model.add_affine_constraint(map_diag(X) == ones(n,1));
     const = -1/(kpara-1);
     model.add_affine_constraint(X(IU,JU)==const);
     model.add_affine_constraint(X(IE,JE) >= const);
 model.solve;
 %%*********************************************************************

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function pass = test_curl( )
% Test curl

tol = 1e2*chebfunpref().cheb2Prefs.chebfun2eps;

% Check curl of a spherefun is a spherefunv.
u = spherefun(@(x,y,z) x.*y.*z);
u = curl(u);
pass(1) = isa(u,'spherefunv');

% Check curl of the zero function is zero
f = spherefun(@(x,y,z) 0*x);
pass(2) = norm(curl(f)) < tol;

%
% Check curl of a non-zero function gives the correct result
%

% Example 1
u = curl(spherefun(@(x,y,z) (x-y).*z));
% Exact curl
f = spherefun(@(x,y,z) (x-y).*y + z.^2);
g = spherefun(@(x,y,z) (y-x).*x + z.^2);
h = spherefun(@(x,y,z) -(x+y).*z);
exact = spherefunv(f,g,h);
pass(3) = norm(u-exact) < tol;

% Example 2
u = curl(spherefun(@(x,y,z) cos(4*z)));
% Exact curl
f = spherefun(@(x,y,z) -4*y.*sin(4*z));
g = spherefun(@(x,y,z) 4*x.*sin(4*z));
h = spherefun(@(x,y,z) 0*x);
exact = spherefunv(f,g,h);
pass(4) = norm(u-exact) < tol;

% Example 3
u = curl(spherefun(@(x,y,z) cos(4*x)));
% Exact curl
f = spherefun(@(x,y,z) 0*x);
g = spherefun(@(x,y,z) -4*z.*sin(4*x));
h = spherefun(@(x,y,z) 4*y.*sin(4*x));
exact = spherefunv(f,g,h);
pass(5) = norm(u-exact) < 100*tol;

% Example 4
u = curl(spherefun(@(x,y,z) cos(4*y)));
% Exact curl
f = spherefun(@(x,y,z) 4*z.*sin(4*y));
g = spherefun(@(x,y,z) 0*x);
h = spherefun(@(x,y,z) -4*x.*sin(4*y));
exact = spherefunv(f,g,h);
pass(6) = norm(u-exact) < 100*tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: September 9th, 2016
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Make_Your_Geometry_and_Input_Files()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  1024;      % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Lx = 5.0;        % Length of Eulerian Grid in x-Direction
dx = Lx/Nx;      % Grid spatial resolution

%
% Immersed Structure Geometric / Dynamic Parameters %
%
ds = 0.45*dx;               % Lagrangian Pt. Spacing (2x resolution of Eulerian grid)
struct_name = 'insect_HT'; % Name for .vertex, .spring, etc files. (must match what's in 'input2d')

% Call function to construct geometry
[xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,Nx,Lx);
yLag = yLag - 2+0.125;

% Plot Geometry to test
plot(xLag,yLag,'r-'); hold on;
plot(xLag,yLag,'*'); hold on;
xlabel('x'); ylabel('y');
axis([0 Lx 0 Lx]);

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .spring file!
%k_Spring = 2.5e4;                    % Spring stiffness (does not need to be equal for all springs)
%ds_Rest = 0.0;                       % Spring resting length (does not need to be equal for all springs)
%print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name);

% Prints .beam file!
% k_Beam = 0.5;                      % Beam Stiffness (does not need to be equal for all beams)
% C = compute_Curvatures(xLag,yLag)  % Computes curvature of initial configuration
%print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name);

% Prints .target file!
k_Target = 2e5;
print_Lagrangian_Target_Pts(xLag,k_Target,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called 'struct_name'.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag); % Total # of Lag. Pts

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints TARGET points to a file called 'struct_name'.target
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag);

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called 'struct_name'.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N );    % Print # of springs

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            else
                %Case s=N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, 1,   k_Spring, ds_Rest);
            end
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called 'struct_name'.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C(s) );
            else
                %Case s=N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C(s) );
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes "curvature" of starting configuration
%
% NOTE: not curvature in the traditional geometric sense, in the 'discrete'
% sense through cross product.
%
% NOTE: assumes a CLOSED structure
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C = compute_Curvatures(xLag,yLag)

N = length(xLag);
C = zeros( N );

%Note: needs to be done same order as you print .beam file!
for i=1:N

   % Pts Xp -> Xq -> Xr (same as beam force calc.)

   if ( (i > 1) && (i < N) )

        Xp = xLag(i-1); Xq = xLag(i); Xr = xLag(i+1);
        Yp = yLag(i-1); Yq = yLag(i); Yr = yLag(i+1);

   elseif (i==1)

        Xp = xLag(N); Xq = xLag(i); Xr = xLag(i+1);
        Yp = yLag(N); Yq = yLag(i); Yr = yLag(i+1);

   elseif (i==N)

        Xp = xLag(N-1); Xq = xLag(N); Xr = xLag(1);
        Yp = yLag(N-1); Yq = yLag(N); Yr = yLag(1);

   end

   C(i) = (Xr-Xq)*(Yq-Yp) - (Yr-Yq)*(Xq-Xp); %Cross product btwn vectors

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,Nx,Lx)

% ds: Lagrangian pt. spacing
% Nx: Eulerian grid resolution

% Gives elliptical subset for tube side
a = 0.5;    % semi-major axis (horizontal)
b = 0.45;   % semi-minor axis (vertical)
ang = pi/6; % angle start/cuttoff for ellipse
[xE,yE] = give_Me_Elliptical_Geometry(ds,a,b,ang);
xOff = 2*xE(1);

% Gives top and bottom portion of one section of tube
d = 0.2;                  % separation distance
xTube = [xE+3*xOff xE+3*xOff xE+2*xOff xE+2*xOff xE+xOff xE+xOff xE xE];          % top first, then bottom
xTube = xTube + 1;
xTube = xTube(end:-1:1);  % reverse it!
yTube = [yE+d/2 -yE-d/2 yE+d/2 -yE-d/2 yE+d/2 -yE-d/2 yE+d/2 -yE-d/2]; % top first, then bottom
yTube = yTube + Lx/2;

% Give indices for to attach valves
ind_Top = 1;                 % node index for first point on top (on RHS)
ind_Bot = length(xTube)/8+1; % node index for first point on bottom (on RHS)

xFlat = xTube(ind_Top)-Lx/20:ds:xTube(ind_Top);
yFlatTop = yTube(ind_Top)*ones(1,length(xFlat));
yFlatBot = yTube(ind_Bot)*ones(1,length(xFlat));

xLag = [xTube xFlat xFlat];
yLag = [yTube yFlatTop yFlatBot];

% Gives valve geometry
[xV,yV] = give_Valve_Geometry(d,ds);

% TESTING (i.e., PLOT IT, yo!)
%plot(xTube,yTube,'b*'); hold on
%len = length(xTube);
%plot(xTube(len/8+1),yTube(len/8+1),'r*'); hold on
%plot(xTube(1),yTube(1),'k*'); hold on;
%axis([0 5 0 5]);
%pause();

%plot(xV,yV,'g*'); hold on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: gives valve geometry inside a distance, d, between the edges of
% the tube
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xV,yV] = give_Valve_Geometry(d,ds)

% for A*x^2 geometry -> 2 eqns and 2 unknowns
%y(0) = 0 = 0
%y(0.9d/2) = A(0.9/2)^2*d^2 = 1.5*d
%
% ----> A = 1.5 / ( (0.45)^2 * d )
%
A = 1.0 / ( (0.45)^2 * d );

yV = -0:ds/2:0.9/2*d;
xV = -A*yV.^2;

%figure(2)
%plot(xV,yV,'b*'); hold on;
%pause();

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: gives elliptical geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xE,yE] = give_Me_Elliptical_Geometry(ds,a,b,ang)

dtheta = ds/max(a,b); % use minimum of axis lengths

n = 1;
xE(n) = a*cos(ang); yE(n) = b*sin(ang);
ang_n = ang + dtheta;

while ang_n < (pi/2)

    n = n+1;                % update counter
    xE(n) = a*cos(ang_n);   % next xVal along ellipse
    yE(n) = b*sin(ang_n);   % next yVal along ellipse
    ang_n = ang_n + dtheta; % update angle

end
xE(n+1) = a*cos(pi/2);
yE(n+1) = b*sin(pi/2);

xE = [xE -xE(end-1:-1:1)];
yE = [yE yE(end-1:-1:1)];

yE = yE - yE(1);           % Translate down

%plot(xE,yE,'*'); hold on;
%axis([-a a -a a]);

xLength = abs( 2*xE(1) );
fprintf('\nThe x-Length of the tube section is: %d\n',xLength);

author: nickabattista


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Universitat Politecnica de Catalunya BarcelonaTech (UPC) - Spain
%  University of California Berkeley (UCB) - USA
%
%  Jordi Pont-Tuset <jordi.pont@upc.edu>
%  Pablo Arbelaez <arbelaez@berkeley.edu>
%  June 2014
% ------------------------------------------------------------------------
% This file is part of the MCG package presented in:
%    Arbelaez P, Pont-Tuset J, Barron J, Marques F, Malik J,
%    "Multiscale Combinatorial Grouping,"
%    Computer Vision and Pattern Recognition (CVPR) 2014.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
function root_dir = root_dir()
    % Change suggested by Ahmad Humayun <ahmadh@gatech.edu>, thank you!
    root_dir = fileparts(which(mfilename));
end

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
function varargout = st_options(varargin)
% ST_OPTIONS Application M-file for st_options.fig
%    FIG = ST_OPTIONS launch st_options GUI.
%    ST_OPTIONS('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.5 24-May-2005 15:03:29

if nargin == 1  % LAUNCH GUI

	fig = openfig(mfilename,'reuse');

	% Use system color scheme for figure:
	set(fig,'Color',get(groot,'defaultUicontrolBackgroundColor'));

	% Generate a structure of handles to pass to callbacks, and store it.
	handles = guihandles(fig);
	guidata(fig, handles);

  % Set figure position
  movegui(fig, 'center');

  if varargin{1} == 1
    set(handles.lblSpacingY, 'String', 'Latitude spacing [deg]:');
    set(handles.txtSpacingY, 'String', '0.1');
    set(handles.lblSpacingX, 'String', 'Longitude spacing [deg]:');
    set(handles.txtSpacingX, 'String', '0.1');
  else
    set(handles.lblSpacingY, 'String', 'Depth spacing [km]:');
    set(handles.txtSpacingY, 'String', '1');
    set(handles.lblSpacingX, 'String', 'Horizontal spacing [km]:');
    set(handles.txtSpacingX, 'String', '1');
  end

	% Wait for callbacks to run and window to be dismissed:
	uiwait(fig);

	if nargout > 0
		varargout{1} = fig;
	end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

	try
		if (nargout)
			[varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
		else
			feval(varargin{:}); % FEVAL switchyard
		end
	catch ME
		disp(ME.message);
	end

end

% --------------------------------------------------------------------
function varargout = Gridding_Callback(h, eventdata, handles, varargin)

set(handles.radNumber, 'Value', (varargin{1} == 1));
set(handles.radRadius, 'Value', (varargin{1} == 2));
set(handles.radRectangle, 'Value', (varargin{1} == 3));

% --------------------------------------------------------------------
function varargout = btnCancel_Callback(h, eventdata, handles, varargin)

handles.answer = 0;
guidata(h, handles);
uiresume(handles.figOptions);

% --------------------------------------------------------------------
function varargout = btnOk_Callback(h, eventdata, handles, varargin)

% Validation of inputs
% --------------------

% % Hypothesis values
% if get(handles.radHNumber, 'Value') == 1
%   if ~gui_IsStringPositiveInteger(get(handles.txtHNumber, 'String'), 'Number of events (hypothesis)')
%     return;
%   end
% elseif get(handles.radHRadius, 'Value') == 1
%   if ~gui_IsStringPositiveFloat(get(handles.txtHRadius, 'String'), 'Constant radius (hypothesis)')
%     return;
%   end
% elseif get(handles.radHRectangle, 'Value') == 1
%   if ~gui_IsStringPositiveFloat(get(handles.txtHSizeRectX, 'String'), 'Rectangle size X (hypothesis)')
%     return;
%   end
%   if ~gui_IsStringPositiveFloat(get(handles.txtHSizeRectY, 'String'), 'Rectangle size Y (hypothesis)')
%     return;
%   end
% end
% if ~gui_IsStringPositiveInteger(get(handles.txtHMinimumNumber, 'String'), 'Minimum number of events (hypothesis)')
%   return;
% end
% if ~gui_IsStringPositiveFloat(get(handles.txtHMaximumRadius, 'String'), 'Maximum radius (hypothesis)')
%   return;
% end
%
% % Null hypothesis values
% if get(handles.radNNumber, 'Value') == 1
%   if ~gui_IsStringPositiveInteger(get(handles.txtNNumber, 'String'), 'Number of events (null hypothesis)')
%     return;
%   end
% elseif get(handles.radNRadius, 'Value') == 1
%   if ~gui_IsStringPositiveFloat(get(handles.txtNRadius, 'String'), 'Constant radius (null hypothesis)')
%     return;
%   end
% elseif get(handles.radNRectangle, 'Value') == 1
%   if ~gui_IsStringPositiveFloat(get(handles.txtNSizeRectX, 'String'), 'Rectangle size X (null hypothesis)')
%     return;
%   end
%   if ~gui_IsStringPositiveFloat(get(handles.txtNSizeRectY, 'String'), 'Rectangle size Y (null hypothesis)')
%     return;
%   end
% end
% if ~gui_IsStringPositiveInteger(get(handles.txtNMinimumNumber, 'String'), 'Minimum number of events (null hypothesis)')
%   return;
% end
% if ~gui_IsStringPositiveFloat(get(handles.txtNMaximumRadius, 'String'), 'Maximum radius (null hypothesis)')
%   return;
% end
%
% % Testing values
% if get(handles.radTNumber, 'Value') == 1
%   if ~gui_IsStringPositiveInteger(get(handles.txtTNumber, 'String'), 'Number of events (testing)')
%     return;
%   end
% elseif get(handles.radTRadius, 'Value') == 1
%   if ~gui_IsStringPositiveFloat(get(handles.txtTRadius, 'String'), 'Constant radius (testing)')
%     return;
%   end
% elseif get(handles.radTRectangle, 'Value') == 1
%   if ~gui_IsStringPositiveFloat(get(handles.txtTSizeRectX, 'String'), 'Rectangle size X (testing)')
%     return;
%   end
%   if ~gui_IsStringPositiveFloat(get(handles.txtTSizeRectY, 'String'), 'Rectangle size Y (testing)')
%     return;
%   end
% end
% if ~gui_IsStringPositiveInteger(get(handles.txtTMinimumNumber, 'String'), 'Minimum number of events (testing)')
%   return;
% end
% if ~gui_IsStringPositiveFloat(get(handles.txtTMaximumRadius, 'String'), 'Maximum radius (testing)')
%   return;
% end
%
% % Options
% if ~gui_IsStringPositiveFloat(get(handles.txtSpacingX, 'String'), 'Spacing X')
%   return;
% end
% if ~gui_IsStringPositiveFloat(get(handles.txtSpacingY, 'String'), 'Spacing Y')
%   return;
% end
% if ~gui_IsStringPositiveFloat(get(handles.txtSplitTime, 'String'), 'Splitting time')
%   return;
% end
% if get(handles.chkForecastPeriod, 'Value') == 1
%   if ~gui_IsStringPositiveFloat(get(handles.txtForecastPeriod, 'String'), 'Forecast period')
%     return;
%   end
% end
% if get(handles.chkLearningPeriod, 'Value') == 1
%   if ~gui_IsStringPositiveFloat(get(handles.txtLearningPeriod, 'String'), 'Learning period')
%     return;
%   end
% end
% if get(handles.chkRandomNode, 'Value') == 1
%   if ~gui_IsStringPositiveInteger(get(handles.txtNumberCalculationNode, 'String'), 'Number calculations per node')
%     return;
%   end
% end
% if get(handles.chkMinMag, 'Value') == 1
%   if ~gui_IsStringPositiveFloat(get(handles.txtMinMag, 'String'), 'Minimum magnitude')
%     return;
%   end
% end
% if ~gui_IsStringPositiveFloat(get(handles.txtMaxMag, 'String'), 'Maximum magnitude')
%   return;
% end

% OK, proceed
handles.answer = 1;
guidata(h, handles);
uiresume(handles.figOptions);

% --------------------------------------------------------------------
function varargout = chkBootstrapSecond_Callback(h, eventdata, handles, varargin)

% --------------------------------------------------------------------
function varargout = txtNumberBootstrapsSecond_Callback(h, eventdata, handles, varargin)

% --- Executes on button press in chkSaveParameter.
function chkSaveParameter_Callback(hObject, eventdata, handles)

if get(handles.chkSaveParameter, 'Value') == 1
  set(handles.btnSaveParameter, 'Enable', 'on');
else
  set(handles.btnSaveParameter, 'Enable', 'off');
end

% --- Executes on button press in btnSaveParameter.
function btnSaveParameter_Callback(hObject, eventdata, handles)

[newfile, newpath] = uiputfile('*.mat', 'Choose parameter file');
% Cancel pressed?
if isequal(newfile, 0)  ||  isequal(newpath, 0)
  return;
end
% Everything ok?
newfile = [newpath newfile];
if length(newfile) > 1
  set(handles.lblSaveParameter, 'String', newfile);
end

function lblSaveParameter_Callback(hObject, eventdata, handles)
% hObject    handle to lblSaveParameter (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of lblSaveParameter as text
%        str2double(get(hObject,'String')) returns contents of lblSaveParameter as a double

author: CelsoReyes


text:
function out = SY_SlidingWindow(y,windowStat,acrossWinStat,numSeg,incMove)
% SY_SlidingWindow  Sliding window measures of stationarity.
%
% This function is based on sliding a window along the time series, measuring
% some quantity in each window, and outputting some summary of this set of local
% estimates of that quantity.
%
% Another way of saying it: calculate 'windowStat' in each window, and computes
% 'acrossWinStat' for the set of statistics calculated in each window.
%
%---INPUTS:
%
% y, the input time series
%
% windowStat, the measure to calculate in each window:
%               (i) 'mean', mean
%               (ii) 'std', standard deviation
%               (iii) 'ent', distribution entropy
%               (iv) 'mom3', skewness
%               (v) 'mom4', kurtosis
%               (vi) 'mom5', the fifth moment of the distribution
%               (vii) 'lillie', the p-value for a Lilliefors Gaussianity test
%               (viii) 'AC1', the lag-1 autocorrelation
%               (ix) 'apen', Approximate Entropy, ApEn(1,0.2)
%               (ix) 'sampen', Sample Entropy, SampEn(2,0.1)
%
% acrossWinStat, controls how the obtained sequence of local estimates is
%                   compared (as a ratio to the full time series):
%                       (i) 'std': standard deviation
%                       (ii) 'ent' histogram entropy
%                       (iii) 'apen': Approximate Entropy, ApEn(1,0.2)
%                       (iii) 'sampen': Sample Entropy, SampEn(2,0.1)
%
% numSeg, the number of segments to divide the time series up into, thus
%       controlling the window length
%
% incMove, the increment to move the window at each iteration, as 1/fraction of the
%       window length (e.g., incMove = 2, means the window moves half the length of the
%       window at each increment)
%
% NOTE: SY_SlidingWindow(y,'mean','std',X,1) is the same as StatAvX, computed as
%                       SY_StatAv(y,'seg',X);
% cf. "Heart rate control in normal and aborted-SIDS infants", S. M. Pincus et al.
%           Am J. Physiol. Regul. Integr. Comp. Physiol. 264(3) R638 (1993)

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

doPlot = false; % plot outputs

% ------------------------------------------------------------------------------
% Check Inputs
% ------------------------------------------------------------------------------

if nargin < 2 || isempty(windowStat)
    windowStat = 'mean'; % measure within each window
end
if nargin < 3 || isempty(acrossWinStat)
    acrossWinStat = 'std'; % measure across all windows
end
if nargin < 4 || isempty(numSeg)
    numSeg = 5;
end
if nargin < 5 || isempty(incMove)
    incMove = 2;
end

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

winLength = floor(length(y)/numSeg); % length of window
if winLength==0
    warning('Time-series of length %u is too short for %u windows',length(y),numSeg);
    out = NaN;
    return
end
inc = floor(winLength/incMove); % increment to move at each step
% If increment rounded down to zero, prop it up:
if inc == 0
    inc = 1;
end

numSteps = (floor((length(y)-winLength)/inc)+1);
qs = zeros(numSteps,1);

% Convert a step index (stepInd) to a range of indices corresponding to that window:
getWindow = @(stepInd) ((stepInd-1)*inc + 1:(stepInd-1)*inc + winLength);

switch windowStat
    case 'mean' % Sliding window mean
        for i = 1:numSteps
            qs(i) = mean(y(getWindow(i)));
        end
    case 'std' % Sliding window std
        for i = 1:numSteps
            qs(i) = std(y(getWindow(i)));
        end
    case 'ent' % Sliding window distributional entropy
        for i = 1:numSteps
            qs(i) = EN_DistributionEntropy(y(getWindow(i)),'ks',[]);
        end
    case 'apen' % Sliding window ApEn
        for i = 1:numSteps
            qs(i) = EN_ApEn(y(getWindow(i)),1,0.2);
        end
    case 'sampen' % Sliding window SampEn
        for i = 1:numSteps
            sampEn_struct = EN_SampEn(y(getWindow(i)),1,0.1);
            qs(i) = sampEn_struct.sampen1;
        end
    case 'mom3' % Third moment
        for i = 1:numSteps
            qs(i) = DN_Moments(y(getWindow(i)),3);
        end
    case 'mom4' % Fourth moment
        for i = 1:numSteps
            qs(i) = DN_Moments(y(getWindow(i)),4);
        end
    case 'mom5' % Fifth moment
        for i = 1:numSteps
            qs(i) = DN_Moments(y(getWindow(i)),5);
        end
    case 'lillie' % Lilliefors test
        for i = 1:numSteps
            qs(i) = HT_DistributionTest(y(getWindow(i)),'lillie','norm');
        end
    case 'AC1' % Lag-1 autocorrelation
        for i = 1:numSteps
            qs(i) = CO_AutoCorr(y(getWindow(i)),1,'Fourier');
        end
    otherwise
        error('Unknown statistic ''%s''',windowStat)
end

% Check for all errors (e.g., short time series):
if all(isnan(qs))
    warning('These sliding window settings are not suitable for this time series');
    out = NaN;
    return
end

% ------------------------------------------------------------------------------
% Plot
% ------------------------------------------------------------------------------
if doPlot
    figure('color','w'); box('on');
    plot(round(winLength/2):inc:(numSteps-1)*inc+round(winLength/2),qs,'r');
end

% ------------------------------------------------------------------------------
% Compute the output statistic
% ------------------------------------------------------------------------------
switch acrossWinStat
    case 'std'
        % normalized by std of full time series
        out = std(qs)/std(y);
    case 'apen'
        out = EN_ApEn(qs,1,0.2); % ApEn of the sliding window measures
    case 'sampen'
        sampEn_struct = EN_SampEn(qs,2,0.15);
        out = sampEn_struct.quadSampEn1;
    case 'ent'
        % get a load of statistics from kernel-smoothed distribution (inefficient since only one is used)
        kssimpouts = DN_FitKernelSmooth(qs);
        out = kssimpouts.entropy; % distributional entropy
    otherwise
        error('Unknown statistic: ''%s''.',acrossWinStat)
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "Hill+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the spring attributes!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function springs_info = update_Springs(dt,current_time,xLag,yLag,springs_info)

%springs_info: col 1: starting spring pt (by lag. discretization)
%              col 2: ending spring pt. (by lag. discretization)
%              col 3: spring stiffness
%              col 4: spring resting lengths

t = current_time;   %CURRENT TIME IN SIM.
ds1 = 3.906250e-03; %PHASE 1 Lag. Spacing
ds2 = 3.342919e-03; %PHASE 2 Lag. Spacing

Nbody = 106;                    %Number of springs before bell connections
RL = give_Resting_Lengths();    %Resting Lengths for Phase 1 (col 1) and Phase 2 (col 2)

% CHANGE RESTING LENGTH BTWN SIDES OF JELLYFISH BELL
for i=1:length( springs_info(:,4) )

    if i <= Nbody
        springs_info(i,4) = ds1 + 5*t*(ds2-ds1);
    else
        ii = i-Nbody;
        springs_info(i,4) = RL(ii,1) + 5*t*( RL(ii,2) - RL(ii,1) );
        %springs_info(107+(i-1),4) = distVec(i)*abs( cos(4*pi*current_time) );
        %springs_info(107+(i-1),4) = distVec(i) * ( 1 - 0.75*( sin(4*pi*current_time) ) );
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: give Resting Lengths (hardcoded!)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function RL = give_Resting_Lengths()

RL = [0.150000000000000   0.075000000000000
   0.149926741000240   0.074973170433545
   0.149706855857594   0.074892624017213
   0.149340020675537   0.074758187724185
   0.148825697410294   0.074569727429823
   0.148163136564709   0.074326763924451
   0.147351380996092   0.074028781527991
   0.146388643571571   0.073675421427274
   0.145274017365354   0.073265752667258
   0.144005152820154   0.072798927606369
   0.142580953464829   0.072274183901744
   0.140998544348621   0.071690309088121
   0.139255611062067   0.071045912762170
   0.137348541189626   0.070339929225114
   0.135274492318472   0.069570469632096
   0.133030516962324   0.068735403317456
   0.130612168717920   0.067832977122889
   0.128014600466564   0.066860709378195
   0.125234180745507   0.065816586316769
   0.122265614832230   0.064697078425658
   0.119103315657513   0.063499067862592
   0.115741451424401   0.062218444447808
   0.112173999448214   0.060850772863983
   0.108392818785359   0.059391311235482
   0.104393528844001   0.057835033327000
   0.100165653754868   0.056174749467196
   0.095702659718569   0.054402747848968
   0.090998237021989   0.052510837712120
   0.086041701055792   0.050487212358191
   0.080827140474999   0.048320369497431
   0.075349280754860   0.045994841644143
   0.069598549422225   0.043493120994599
   0.063573706279126   0.040793235757662
   0.057274941348722   0.037865974506005
   0.050704050274530   0.034676918247010
   0.043870119404512   0.031179908592683
   0.036787309077097   0.027316537892343
   0.029477909226729   0.023009349350403
   0.021975429586313   0.018167707789027
   0.014318966133393   0.012694598507574
   0.006558552552061   0.006558010776142];


text:
function img = image_cap_quant(img, k, klower)
    for i = 1 : size(img, 5)
        tmp = img(:,:,:,:, i);
        sorted = sort(tmp(:), 'ascend');
        thrval = sorted(round(k * end));
        tmp(tmp > thrval) = thrval;

        thrval = sorted(round(klower * end));
        tmp(tmp < thrval) = thrval;

        img(:,:,:,:, i) = tmp;
    end
end

repo: pTVreg


text:
%IM_SCALE Fixed mapping scaling binary images to a giving fraction of pixels 'on'
%
%   B = IM_SCALE(A,P)
%   B = A*IM_SCALE([],P)
%   B = A*IM_SCALE(P)
%
% B is a zoomed in / out version of A such that about a fraction
% P of the image pixels is 'on' (1).
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES, IM_BOX, IM_CENTER

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function b = im_scale(varargin)

	argin = shiftargin(varargin,'vector');
  argin = setdefaults(argin,[],0.5);
  if mapping_task(argin,'definition')
    b = define_mapping(argin,'fixed');
    b = setname(b,'Image scale');
  else
    [a,p] = deal(argin{:});
    if isdataset(a)
      error('Command cannot be used for datasets as it may change image size')
    elseif isdatafile(a)
      isobjim(a);
      b = filtim(a,mfilename,{p});
      b = setfeatsize(b,getfeatsize(a));
    elseif isa(a,'double') || isa(a,'dip_image') % here we have a single image
      sca = sqrt(p/mean(a(:)));
      sa = size(a);
      c = imresize(double(a),round(sca*sa),'nearest');
      sc = size(c);
      d = abs(floor((sc - sa)/2));
      if sca < 1
        b = zeros(size(a));
        b(d(1)+1:d(1)+sc(1),d(2)+1:d(2)+sc(2)) = c;
      else
        b = c(d(1)+1:d(1)+sa(1),d(2)+1:d(2)+sa(2));
      end
    end
  end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
a=[2,4;3,-5;1,2;2,1];
b=[11;3;6;7];
solution=a\b

author: Eurus-Holmes


text:
function [ diag, off, seed ] = r8sr_random ( n, nz, row, col, seed )

%*****************************************************************************80
%
%% R8SR_RANDOM randomizes a R8SR matrix.
%
%  Discussion:
%
%    The R8SR storage format stores the diagonal of a sparse matrix in DIAG.
%    The off-diagonal entries of row I are stored in entries ROW(I)
%    through ROW(I+1)-1 of OFF.  COL(J) records the column index
%    of the entry in A(J).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer NZ, the number of offdiagonal nonzero elements in A.
%
%    Input, integer ROW(N+1).  The nonzero offdiagonal elements of row I of A
%    are contained in A(ROW(I)) through A(ROW(I+1)-1).
%
%    Input, integer COL(NZ), contains the column index of the element
%    in the corresponding position in A.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real DIAG(N), the diagonal elements of A.
%
%    Output, real OFF(NZ), the off-diagonal elements of A.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  for i = 1 : n
    [ diag(i), seed ] = r8_uniform_01 ( seed );
    for j = row(i) : row(i+1)-1
      [ off(j), seed ] = r8_uniform_01 ( seed );
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function cc_grids_constrained_display ( )

%*****************************************************************************80
%
%% CC_GRIDS_CONSTRAINED_DISPLAY displays grids generated by CC_GRIDS_CONSTRAINED.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 November 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_GRIDS_CONSTRAINED_DISPLAY:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Display the 2D Clenshaw-Curtis grids\n' );
  fprintf ( 1, '  generated by CC_GRIDS_CONSTRAINED.\n' );

  dim_num = 2;

  while ( 1 )
%
%  Get user input.
%
    q_max = input ( 'Enter Q_MAX or RETURN to exit;' );

    if ( isempty ( q_max ) )
      break
    end

    if ( q_max < dim_num )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  We require DIM_NUM <= Q_MAX!\n' );
      continue
    end

    alpha = input ( 'Enter [ ALPHA1, ALPHA2 ] or RETURN to exit;' );

    if ( isempty ( alpha ) )
      break
    end

    order_min = input ( 'Enter [ ORDER_MIN1, ORDER_MIN2 ] or RETURN to exit;' );

    if ( isempty ( order_min ) )
      break
    end

    order_max = input ( 'Enter [ ORDER_MAX1, ORDER_MAX2 ] or RETURN to exit;' );

    if ( isempty ( order_max ) )
      break
    end
%
%  Compute data.
%
    [ grid_num, point_num ] = cc_grids_constrained_size ( dim_num, ...
      q_max, alpha, order_min, order_max );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Number of grids is %d\n', grid_num );
    fprintf ( 1, '  Number of points is %d\n', point_num );

    [ grid_order, grid_point ] = cc_grids_constrained ( dim_num, ...
      q_max, alpha, order_min, order_max, grid_num, point_num );

    clf
%
%  We have to name the axes in order to control the grid.
%
    axes_handle = axes;
%
%  Plot the points.
%
    handle = scatter ( grid_point(1,:), grid_point(2,:), 'filled' );
%
%  Force the plotting region to be square, not rectangular.
%
    axis square
%
%  Request grid lines.
%
    grid on
%
%  Specify the location of the grid lines, and suppress labeling.
%
    set ( axes_handle, 'xtick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
    set ( axes_handle, 'xticklabel', [] );
    set ( axes_handle, 'ytick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
    set ( axes_handle, 'yticklabel', [] );
%
%  Make the plotting region slightly bigger than the data.
%
    axis ( [ -1.1, 1.1, -1.1, 1.1 ] )
%
%  Title
%
    s = sprintf ( '%d <= Q <= %d', q_min, q_max );
    title ( s );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_GRIDS_CONSTRAINED_DISPLAY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function dtiSaveFiberClusters(clusterlabels, fibergroupsvector, FiberGroup, versionNum, coordinateSpace)
% Save fibergroups (clusters of fibers)
%
%   dtiSaveFiberClusters(clusterlabels, fibergroupsvector, FiberGroup,versionNum, coordinateSpace)
%
%Save fibergroups (clusters of fibers) specified as a vector with cluster numbers -- in the
%cluster structure clusterlabels. FG is the original (clustered eventually) set of fibers.
%E.g., saveFiberClusters(clusterlabels, [4 93 87], fibergroup1, versionNum, coordinateSpace)
%
%To form fibergroupsvector use something like this:
% fibergroupsvector=[]
%for clust=1:max(clusterlabels)
% if size(find(clusterlabels==clust), 1)>10
% fibergroupsvector=[fibergroupsvector clust]
% end
% end

%ER 2007 SCSNL; edited 02/2008

%Figure out span for colors
nfibergroups=length(fibergroupsvector);
%Form vectors of colors:
vecofcolor=int16(colormap(hsv(nfibergroups))*255);
iteration=0; %iterations;

for fibergroupIndex=fibergroupsvector

iteration=iteration+1;
fg.name=['Fibergroup' num2str(fibergroupIndex)];
fg.colorRgb=vecofcolor(iteration, :);

fg.thickness=FiberGroup.thickness;
fg.visible=FiberGroup.visible;
if ~isempty(FiberGroup.seeds)
    fg.seeds=FiberGroup.seeds(clusterlabels==fibergroupIndex, :);
end
    fg.seedRadius=FiberGroup.seedRadius;
fg.seedVoxelOffsets=FiberGroup.seedVoxelOffsets;
fg.params=FiberGroup.params;

fg.fibers=FiberGroup.fibers(clusterlabels==fibergroupIndex);

filename=['Fibergroup' num2str(fibergroupIndex)];
save(filename, 'fg', 'coordinateSpace', 'versionNum');

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function r = mpower(a,b)
%MPOWER       Implements  a ^ b  for intervals
%
% either a and b are (interval) scalar or, b is scalar integer
%

% written  10/16/98     S.M. Rump
% modified 12/30/98     S.M. Rump  improved performance and even exponent
% modified 05/20/02     S.M. Rump  interval integer exponent checked
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    NaN corrected
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
%

  % for both a and b scalars, use .^ with improved diameter for even exponent
  if ( prod(size(a))==1 ) & ( prod(size(b))==1 )
    r = a .^ b ;
    return
  end

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  % Check integer exponent
  if isa(b,'intval') & all(inf(b)==sup(b))
    b = b.inf;
    if ~isreal(b)
      error('invalid call of intval mpower ^')
    end
    a = intval(a);
  end

  if isa(b,'double') & isreal(b) & prod(size(b))==1 & b==round(b)
    [m n] = size(a);
    if m~=n
      error('intval mpower of non-square matrix')
    end
    if b==0
      if issparse(a)
        r = intval(speye(size(a)));
      else
        r = intval(eye(size(a)));
      end
      index = isnan(a) & ( b==0 );
      if any(index(:))
        %VVVV  r(index) = NaN;
        s.type = '()'; s.subs = {index}; r = subsasgn(r,s,NaN);
        %AAAA  Matlab bug fix
      end
    else                        % b is integer
      b_is_negative = b<0;
      b = abs(b) - 1;           % abs(b) is at least 1
      r = a;
      while b>0
        if mod(b,2)==1
          r = r*a;
        end
        b = floor(b/2);
        if b~=0
          a = a*a;
        end
      end
      if b_is_negative
        r = inv(r);
      end
    end
  else
    error('invalid call of intval mpower ^')
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function febio(varargin)

% function febio(febFileName)
% ------------------------------------------------------------------------
% This function calls FEBio from MATLAB and runs it in the command window.
% If no input is provided this commend simply triggers FEBio to run without
% an input file.
% If the input is a valid FEBio model input file path then FEBio will run
% the analysis for that file.
% The input may also be a structure with the following default fields:
%   defaultOptionStruct.FEBioPath=getFEBioPath; %The FEBio path
%   defaultOptionStruct.run_filename=[]; %The .feb file name
%   defaultOptionStruct.run_logname=[]; %The log file name
%   defaultOptionStruct.runMode='internal'; %Running internally
%
%
% 2020/11/24 Created
% ------------------------------------------------------------------------

%% Parse input

%Create default structure
defaultOptionStruct.FEBioPath=getFEBioPath; %Get FEBio path
defaultOptionStruct.run_filename=[];
defaultOptionStruct.run_logname=[];
defaultOptionStruct.runMode='internal';

switch nargin
    case 0
        inputVar=[];
    case 1
        inputVar=varargin{1};
end

if ~isempty(inputVar)
    if ~isstruct(inputVar) %Assume a feb file name was provided
        FEBioRunStruct.run_filename=inputVar;
        [filePath,fileName,~]=fileparts(inputVar);
        logFileName=fullfile(filePath,[fileName,'.txt']);
        FEBioRunStruct.run_logname=logFileName;
    else
        FEBioRunStruct=inputVar;
    end
    %Complete struct if there are missing fields
    FEBioRunStruct=structComplete(FEBioRunStruct,defaultOptionStruct,0);
else
    FEBioRunStruct=[];
end

%% Run FEBio

if isstruct(FEBioRunStruct) && ~isfield(FEBioRunStruct,'run_string')
    switch FEBioRunStruct.runMode
        case 'external_old'
            FEBioRunStruct.run_string=['"',FEBioRunStruct.FEBioPath,'" "',FEBioRunStruct.run_filename,'" -o "',FEBioRunStruct.run_logname,'" &'];
        case 'external'
            if ismac % Code to run on Mac plaform
                %TO DO IMPROVE THIS LINE FOR MAC
                FEBioRunStruct.run_string=['nice "',FEBioRunStruct.FEBioPath,'" "',FEBioRunStruct.run_filename,'" &'];
            elseif isunix && ~ismac % Code to run on Linux plaform
                FEBioRunStruct.run_string=['nice "',FEBioRunStruct.FEBioPath,'" "',FEBioRunStruct.run_filename,'" &'];
            elseif ispc % Code to run on Windows platform
                FEBioRunStruct.run_string=['start /min "GIBBON - FEBio" "',FEBioRunStruct.FEBioPath,'" "',FEBioRunStruct.run_filename,'"'];
            else
                warning('Unknown operational system, run command might be inappropriate');
                FEBioRunStruct.run_string=['"',FEBioRunStruct.FEBioPath,'" "',FEBioRunStruct.run_filename,'" &'];
            end
        case 'internal'
            if ispc
                FEBioRunStruct.run_string=['"',FEBioRunStruct.FEBioPath,'" "',FEBioRunStruct.run_filename,'"'];
            elseif isunix
                FEBioRunStruct.run_string=['nice "',FEBioRunStruct.FEBioPath,'" "',FEBioRunStruct.run_filename,'"'];
            else
                FEBioRunStruct.run_string=['"',FEBioRunStruct.FEBioPath,'" "',FEBioRunStruct.run_filename,'"'];
            end
    end
elseif isempty(FEBioRunStruct)
    FEBioRunStruct.run_string=['"',defaultOptionStruct.FEBioPath,'"'];
end

% if isstruct(FEBioRunStruct) && ~isfield(FEBioRunStruct,'run_string')
%     switch FEBioRunStruct.runMode
%         case 'external_old'
%             FEBioRunStruct.run_string=['"',FEBioRunStruct.FEBioPath,'" -i "',FEBioRunStruct.run_filename,'" -o "',FEBioRunStruct.run_logname,'" &'];
%         case 'external'
%             if ismac % Code to run on Mac plaform
%                 %TO DO IMPROVE THIS LINE FOR MAC
%                 FEBioRunStruct.run_string=['nice "',FEBioRunStruct.FEBioPath,'"',' -i "',FEBioRunStruct.run_filename,'" -o "',FEBioRunStruct.run_logname,'" &'];
%             elseif isunix && ~ismac % Code to run on Linux plaform
%                 FEBioRunStruct.run_string=['nice "',FEBioRunStruct.FEBioPath,'"',' -i "',FEBioRunStruct.run_filename,'" -o "',FEBioRunStruct.run_logname,'" &'];
%             elseif ispc % Code to run on Windows platform
%                 FEBioRunStruct.run_string=['start /min "GIBBON - FEBio" "',FEBioRunStruct.FEBioPath,'"',' -i "',FEBioRunStruct.run_filename,'" -o "',FEBioRunStruct.run_logname,'"'];
%             else
%                 warning('Unknown operational system, run command might be inappropriate');
%                 FEBioRunStruct.run_string=['"',FEBioRunStruct.FEBioPath,'"',' -i "',FEBioRunStruct.run_filename,'" -o "',FEBioRunStruct.run_logname,'" &'];
%             end
%         case 'internal'
%             if ispc
%                 FEBioRunStruct.run_string=['"',FEBioRunStruct.FEBioPath,'" -i "',FEBioRunStruct.run_filename,'" -o "',FEBioRunStruct.run_logname,'"'];
%             elseif isunix
%                 FEBioRunStruct.run_string=['nice "',FEBioRunStruct.FEBioPath,'" -i "',FEBioRunStruct.run_filename,'" -o "',FEBioRunStruct.run_logname,'"'];
%             else
%                 FEBioRunStruct.run_string=['"',FEBioRunStruct.FEBioPath,'" -i "',FEBioRunStruct.run_filename,'" -o "',FEBioRunStruct.run_logname,'"'];
%             end
%     end
% elseif isempty(FEBioRunStruct)
%     FEBioRunStruct.run_string=['"',defaultOptionStruct.FEBioPath,'"'];
% end

% FEBioRunStruct.run_string

system(FEBioRunStruct.run_string); %START FEBio NOW!!!!!!!!

end
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%GRADIENTDERICHE  Applies Deriche filter to an image
%
%     dst = cv.GradientDeriche(op, dir)
%     dst = cv.GradientDeriche(op, dir, 'OptionName',optionValue, ...)
%
% ## Input
% * __op__ Source 8-bit or 16-bit image, 1-channel or 3-channel image.
% * __dir__ Filter direction. One of:
%   * __X__
%   * __Y__
%
% ## Output
% * __dst__ result `single` image with same number of channel than `op`.
%
% ## Options
% * __AlphaDerive__ double see paper. default 1.0
% * __AlphaMean__ double see paper. default 1.0
%
% For more details about this implementation, please see [deriche1987].
%
% ## References
% [deriche1987]:
% > Rachid Deriche. "Using Canny's criteria to derive a recursively
% > implemented optimal edge detector". International Journal of Computer
% > Vision, Volume 1 Issue 2, pages 167-187, 1987.
% > http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.476.5736&rep=rep1&type=pdf
%
% See also: imgradientxy
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%kitexture 'Texture Feature Extraction Using LAW Metrics'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros itexture.pane file
%
% Parameters:
% InputFile: i 'Input object ', required: 'input object'
% OutputFile: o 'Output object', required: 'output object'
% Toggle: L5L5 'L5L5', default: 0: 'If set, apply L5L5 texture kernel'
% Toggle: E5L5 'E5L5', default: 0: 'If set, apply E5L5 texture kernel'
% Toggle: S5L5 'S5L5', default: 0: 'If set, apply S5L5 texture kernel'
% Toggle: W5L5 'W5L5', default: 0: 'If set, apply W5L5 texture kernel'
% Toggle: R5L5 'R5L5', default: 0: 'If set, apply R5L5 texture kernel'
% Toggle: L5E5 'L5E5', default: 0: 'If set, apply L5E5 texture kernel'
% Toggle: E5E5 'E5E5', default: 0: 'If set, apply E5E5 texture kernel'
% Toggle: S5E5 'S5E5', default: 0: 'If set, apply S5E5 texture kernel'
% Toggle: W5E5 'W5E5', default: 0: 'If set, apply W5E5 texture kernel'
% Toggle: R5E5 'R5E5', default: 0: 'If set, apply R5E5 texture kernel'
% Toggle: L5S5 'L5S5', default: 0: 'If set, apply L5S5 texture kernel'
% Toggle: E5S5 'E5S5', default: 0: 'If set, apply E5S5 texture kernel'
% Toggle: S5S5 'S5S5', default: 0: 'If set, apply S5S5 texture kernel'
% Toggle: W5S5 'W5S5', default: 0: 'If set, apply W5S5 texture kernel'
% Toggle: R5S5 'R5S5', default: 0: 'If set, apply R5S5 texture kernel'
% Toggle: L5W5 'L5W5', default: 0: 'If set, apply L5W5 texture kernel'
% Toggle: E5W5 'E5W5', default: 0: 'If set, apply E5W5 texture kernel'
% Toggle: S5W5 'S5W5', default: 0: 'If set, apply S5W5 texture kernel'
% Toggle: W5W5 'W5W5', default: 0: 'If set, apply W5W5 texture kernel'
% Toggle: R5W5 'R5W5', default: 0: 'If set, apply R5W5 texture kernel'
% Toggle: L5R5 'L5R5', default: 0: 'If set, apply L5R5 texture kernel'
% Toggle: E5R5 'E5R5', default: 0: 'If set, apply E5R5 texture kernel'
% Toggle: S5R5 'S5R5', default: 0: 'If set, apply S5R5 texture kernel'
% Toggle: W5R5 'W5R5', default: 0: 'If set, apply W5R5 texture kernel'
% Toggle: R5R5 'R5R5', default: 0: 'If set, apply R5R5 texture kernel'
%
% Example: o = kitexture(i, {'i','';'o','';'L5L5',0;'E5L5',0;'S5L5',0;'W5L5',0;'R5L5',0;'L5E5',0;'E5E5',0;'S5E5',0;'W5E5',0;'R5E5',0;'L5S5',0;'E5S5',0;'S5S5',0;'W5S5',0;'R5S5',0;'L5W5',0;'E5W5',0;'S5W5',0;'W5W5',0;'R5W5',0;'L5R5',0;'E5R5',0;'S5R5',0;'W5R5',0;'R5R5',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% itexture - Texture Feature Extraction Using LAW Metrics
%
%  DESCRIPTION
% .I itexture
% will extract texture features from an image by convolving the original image
% with a set of special spatial filter kernels.  The spatial filters are
% comprised of 5 x 5 kernels
% derived from center-weighted vectors defined by the LAWS texture measure.
%
% Each texture measure calculated will be stored as width-height planes along
% the elements direction in the Output Object (o).
% The order in which the (selected) texture measures will be stored in the
% output object are:
% .DS
%         L5L5
%         L5E5
%         L5S5
%         L5W5
%         L5R5
%         E5L5
%         E5E5
%         E5S5
%         E5W5
%         E5R5
%         S5L5
%         S5E5
%         S5S5
%         S5W5
%         S5R5
%         W5L5
%         W5E5
%         W5S5
%         W5W5
%         W5R5
%         R5L5
%         R5E5
%         R5S5
%         R5W5
%         R5R5
% .DE
%
% If the input object elements dimension is greater than 1 then the results
% obtained by applying the first selected kernel to each band are stored by planes
% in the output object, followed by the results obtained by applying the
% next selected texture kernel to all bands, and so on. Thus a 3-band input
% image, operated on by the L5L5 and R5R5 kernels, will produce a 6-band
% output image with the results as follows:
% .DS
%       Band 0 - L5L5 on Band 0 of input image
%       Band 1 - L5L5 on Band 1 of input image
%       Band 2 - L5L5 on Band 2 of input image
%       Band 3 - R5R5 on Band 0 of input image
%       Band 4 - R5R5 on Band 1 of input image
%       Band 5 - R5R5 on Band 2 of input image
% .DE
%
% The five center-weighted vectors are:
% .DS
%          L5 = [ 1  4  6  4  1]
%          E5 = [-1 -2  0  2  1]
%          S5 = [-1  0  2  0 -1]
%          W5 = [-1  2  0 -2  1]
%          R5 = [ 1 -4  6 -4  1]
% .DE
%
% Each 5 x 5 kernel is derived from multiplying a particular combination of
% two of the above vectors.  This results in 25 possible 5 x 5 kernels.
% Note that 10 of the kernels are formed by taking the transpose.  For example,
% L5E5 is a 5 x 5 kernel formed by multiplying the vectors, L5 with E5.  It's
% transpose is E5L5, and is similarly formed by multiplying the vectors,
% E5 with L5.  One of these kernels is more sensitive to horizontal
% changes in texture, while the other is sensitive to vertical changes in
% texture.
%
% The kernels are applied with a centered hotspot.
%
% The sum of the elements of each kernel is zero, which results in the output
% image having a mean of zero.  Therefore the relevant texture information is
% contained in the image variance of the microtexture features.  The LAW
% kernels were designed to be sensitive to structures such as edges, ripples,
% and spots.
%
% Choosing which kernel to use will depend on the nature of the texture of
% interest, and will require trial and error to decide which produces the
% feature with the most discriminating power.  In general, if a texture kernel
% of say, L5S5 is used, one may also want to use its transpose, S5L5 for
% certain types of quasiperiodic variations commonly found in textured images.
%
% The colorspace model for the output object is always set to KNONE.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% .I itexture
% will operate on input objects of any data storage type.
%
% The output object will be of type KLONG, KDOUBLE, or
% KDCOMPLEX, determined by the data type of the input object.
% If the input object is of data type KBYTE, KUBYTE, KSHORT, or KUSHORT,
% then it is converted up to an KLONG image. If the input
% is of type KFLOAT or KDOUBLE, then the result will be KDOUBLE.
% Any complex input type is converted to type KDCOMPLEX.
%
%  REFERENCES
% K. I. Laws, "Rapid texture identification" Proc. SPIE, vol. 238, pp.376-380,
% 1980.
%
% J. Y. Hsiao, and A. A. Sawchuk, "Supervised textured image segmentation
% using feature smoothing and probabilistic relaxation techniques" IEEE Trans.
% Pattern Anal. Machine Intell., vol. 11, No. 12, pp 1279-1292, 1989.
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kitexture(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kitexture(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output';'L5L5', 0;'E5L5', 0;'S5L5', 0;'W5L5', 0;'R5L5', 0;'L5E5', 0;'E5E5', 0;'S5E5', 0;'W5E5', 0;'R5E5', 0;'L5S5', 0;'E5S5', 0;'S5S5', 0;'W5S5', 0;'R5S5', 0;'L5W5', 0;'E5W5', 0;'S5W5', 0;'W5W5', 0;'R5W5', 0;'L5R5', 0;'E5R5', 0;'S5R5', 0;'W5R5', 0;'R5R5', 0};
maxval={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
minval={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
istoggle=[0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'itexture"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [fx,dF_dX,dF_dTheta] = f_vanDerPol(Xt,Theta,ut,inF)
% evolution function for the Van Der Pol (limit cycle) oscillator.
% function [fx,dF_dX,dF_dTheta,d2F_dXdTheta] = f_vanDerPol(Xt,Theta,ut,inF)
% IN:
%   - Xt: current hidden state
%   - Theta: evolution parameter
%   - ut: current input to the system
%   - in: user-defined input structure (containing the time decimation)
% OUT:
%   - fx: the predicted hidden state
%   - dF_dX: the jacobian of the system
%   - dF_dTheta: the derivative of the evolution function w.r.t the
%   evolution parameters

deltat = inF.deltat;

x       = Xt;
mu       = Theta(1);

fx      = zeros(2,1);
fx(1)   = x(2);
fx(2)   = mu.*(1-x(1).^2).*x(2) - x(1);
fx      = deltat.*fx + x;

J       = zeros(2,2);
J(1,:)  = [0,1];
J(2,:)  = [-2.*mu.*x(1).*x(2)-1,mu.*(1-x(1).^2)];
dF_dX   = deltat.*J';

dF_dTheta           = zeros(1,2);
dF_dTheta(1,:)      = [0,(1-x(1).^2).*x(2)];
dF_dTheta           = deltat.*dF_dTheta;

author: MBB-team

repo: VBA-toolbox


text:
function c = mid(a)
%MID          Implements  mid(a)  for intervals (rounded)
%
%   c = mid(a)
%
% mid(a) and rad(a) computed such that
%    alpha  in  < mid(a) , rad(a) >  for all alpha in a
%
%For intervals at least 3 bits wide, the midpoint is always an inner point.
%

% written  10/16/98     S.M. Rump
% modified 06/22/99     S.M. Rump  for sparse matrices
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 11/20/05     S.M. Rump  fast check for rounding to nearest
% modified 09/10/07     S.M. Rump  performance, huge arrays
% modified 10/18/08     S.M. Rump  again huge arrays
% modified 07/23/09     S.M. Rump  changed formula: midpoint now in rnd to nearest
%                                    to make sure mid([1-eps,1+2eps])=1
%                                    (thanks to Gerhard Heindl for pointing to this)
% modified 10/06/09     S.M. Rump  check for rndold
%

  if a.complex                            % complex or thin interval
    c = a.mid;
  else
    e = 1e-30;
    if 1+e==1-e                         % fast check for rounding to nearest
      rndold = 0;
    else
      rndold = getround;
      setround(0)
    end
    % use a.inf + (0.5*a.sup-0.5*a.inf) for correct result in case a.sup-a.inf overflows
    [m,n] = size(a.inf);
    if m*n<2^31                           % input not huge
      c = a.inf + (0.5*a.sup-0.5*a.inf);  % make sure result correct in underflow range
      indexinf = isinf(a.inf);
      indexsup = isinf(a.sup);
      anyindexinf = any(indexinf(:));
      anyindexsup = any(indexsup(:));
      if anyindexinf                      % make sure mid([-inf,x]) is x
        c(indexinf) = a.sup(indexinf);
      end
      if anyindexsup                      % make sure mid([x,inf]) is x
          c(indexsup) = a.inf(indexsup);
      end
      if anyindexinf | anyindexsup        % some components are inf
        c(indexinf & indexsup) = 0;       % make sure mid([-inf,inf]) is 0
      end
    else                                  % take care of huge matrices
      % check some components are inf
      % careful with intervals [0,2] or [-2,0] or [-2,2]
      [Iinf,Jinf,Sinf] = find(a.inf);
      [Isup,Jsup,Ssup] = find(a.sup);
      if ( ~isempty(Iinf) ) | ( ~isempty(Isup) )
        ainfsup = sparse([Iinf;Isup],[Jinf;Jsup],[complex(Sinf,0);complex(0,Ssup)],m,n);
        [I,J,S] = find(ainfsup);
        Sold = S;
        S = real(S) + (0.5*imag(S)-0.5*real(S));
        indexinf = isinf(real(Sold));
        if any(indexinf(:))               % make sure mid([-inf,x]) is x
          S(indexinf) = imag(Sold(indexinf));
        end
        indexsup = isinf(imag(Sold));
        if any(indexsup(:))               % make sure mid([x,inf]) is x
          S(indexsup) = real(Sold(indexsup));
        end
        S(indexinf & indexsup) = 0;       % make sure mid([-inf,inf]) is 0
        c = sparse(I,J,S,m,n);
      else
        c = a.inf + 0.5*(a.sup-a.inf);    % make sure result correct in underflow range
      end
    end
    if rndold
      setround(rndold)                    % reset rounding
    end
  end


text:
function [sts, val] = subsasgn_check(item,subs,val)

% function [sts, val] = subsasgn_check(item,subs,val)
% Check whether type of val conforms to configuration tree specification.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: subsasgn_check.m 2512 2008-12-01 13:21:29Z volkmar $

rev = '$Rev: 2512 $'; %#ok

sts = true;

% check, whether arguments for 'val' are cfg_items
switch subs(1).subs
    case {'val'}
        sts = ~cfg_get_defaults('cfg_item.checkval') || subsasgn_check_valcfg(subs,val,[0 Inf]);
end;

author: spm


text:
function mean = burr_mean ( a, b, c, d )

%*****************************************************************************80
%
%% BURR_MEAN returns the mean of the Burr PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, D, the parameters of the PDF.
%    0 < B,
%    0 < C.
%
%    Output, real MEAN, the mean of the PDF.
%
  mean = a + b * gamma ( 1.0 - 1.0 / c ) * gamma ( d + 1.0 / c ) / gamma ( d );

  return
end


text:
function val = file2mat(a,varargin)
% Function for reading from file_array objects
% FORMAT val = file2mat(a,ind1,ind2,ind3,...)
% a      - file_array object
% indx   - indices for dimension x (int64)
% val    - the read values
%
% This function is normally called by file_array/subsref.
%__________________________________________________________________________
% Copyright (C) 2005-2017 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: file2mat.m 7147 2017-08-03 14:07:01Z spm $

%-This is merely the help file for the compiled routine
error('file2mat.c not compiled - see Makefile');