text:
function priors = spm_cfg_eeg_inv_priors
% Configuration file to set up priors for M/EEG source reconstruction
%__________________________________________________________________________
% Copyright (C) 2010-2017 Wellcome Trust Centre for Neuroimaging

% Gareth Barnes
% $Id: spm_cfg_eeg_inv_priors.m 7110 2017-06-15 11:45:05Z guillaume $

D = cfg_files;
D.tag = 'D';
D.name = 'M/EEG datasets';
D.filter = 'mat';
D.num = [1 Inf];
D.help = {'Select the M/EEG mat files.'};

val = cfg_entry;
val.tag = 'val';
val.name = 'Inversion index';
val.strtype = 'n';
val.help = {'Index of the cell in D.inv where the forward model can be found and the results will be stored.'};
val.val = {1};

priorname = cfg_entry;
priorname.tag = 'priorname';
priorname.name = 'Prefix for prior';
priorname.strtype = 's';
priorname.num = [1 Inf];
priorname.val = {'priorset1'};
priorname.help = {'Prefix for prior directory'};

sensorlevel = cfg_entry;
sensorlevel.tag = 'sensorlevel';
sensorlevel.name = 'RMS Sensor level noise';
sensorlevel.strtype = 'r';
sensorlevel.help = {'Assuming white sensor noise of same magnitude on all channels. Units fT for MEG. eg 10fT/rt(Hz) white noise for 100Hz bandwidth would be 100fT rms'};
sensorlevel.val = {100};

FWHMmm = cfg_entry;
FWHMmm.tag = 'FWHMmm';
FWHMmm.name = 'Patch smoothness';
FWHMmm.strtype = 'r';
FWHMmm.num = [1 1];
FWHMmm.val = {[4]};
FWHMmm.help = {'FWHM on cortex in mm. NB. Ignored for EBB and IID. Square windowed for COH '};

popular=cfg_menu;
popular.tag='popular';
popular.name='prior types';
popular.help={' Select from popular priors: IID (min norm), COH (LORETA) or EBB (beamforming)'};
popular.labels={'IID','COH','EBB','sparseEBB'};
popular.values={'IID','COH','EBB','sparseEBB'};
popular.val={'IID'};

popularpars=cfg_branch;
popularpars.tag='popularpars';
popularpars.name='Popular priors';
popularpars.help={' Select from popular priors and kernel smoothness (no smoothing used for IID)'};
popularpars.val={popular,FWHMmm};

npatches = cfg_entry;
npatches.tag = 'npatches';
npatches.name = 'Number of randomly selected patches';
npatches.strtype = 'i';
npatches.num = [1 1];
npatches.val = {[512]};
npatches.help = {'Number of randomly centred patches (priors) on each iteration'};

niter = cfg_entry;
niter.tag = 'niter';
niter.name = 'Number of iterations';
niter.strtype = 'i';
niter.num = [1 1];
niter.val = {[1]};
niter.help = {'Number of iterations'};

probfile = cfg_entry;
probfile.tag = 'probfile';
probfile.name = 'Prior to bias distribution of patches X';
probfile.strtype = 's';
%probfile.num = [1 1];
probfile.val = {''};
probfile.help = {'Name of patch probability file'};

randpatch = cfg_branch;
randpatch.tag = 'randpatch';
randpatch.name = 'Random Patches';
randpatch.help = {'Define random patches'};
randpatch.val  = {npatches,niter,FWHMmm,probfile};

fixedfile = cfg_files;
fixedfile.tag = 'fixedfile';
fixedfile.name = 'Patch definition file ';
fixedfile.filter = 'mat';
fixedfile.num = [1 1];
fixedfile.help = {'Select patch definition file (mat file with variable Ip: rows are iterations columns are patch indices '};
fixedfile.val={{''}};

fixedrows = cfg_entry;
fixedrows.tag = 'fixedrows';
fixedrows.name = 'Rows from fixed patch file';
fixedrows.strtype = 'i';
fixedrows.num = [1 2];
fixedrows.val = {[1 Inf]};
fixedrows.help = {'Select first and last row of patch sets to use from file'};

fixedpatch = cfg_branch;
fixedpatch.tag = 'fixedpatch';
fixedpatch.name = 'Fixed Patches';
fixedpatch.help = {'Define fixed patches'};
fixedpatch.val  = {fixedfile,fixedrows,FWHMmm};

sparsify = cfg_entry;
sparsify.tag = 'sparsify';
sparsify.name = 'Sparsify to N largest local maxima';
sparsify.strtype = 'i';
sparsify.num = [1 2];
sparsify.val = {[1 Inf]};
sparsify.help = {'Select first and last row of patch sets to use from file'};

space = cfg_menu;
space.tag = 'space';
space.name = 'Prior image space';
space.help = {'Space of the mask image.'};
space.labels = {'MNI', 'Native'};
space.values = {1, 0};
space.val = {1};

clean = cfg_menu;
clean.tag = 'clean';
clean.name = 'Start fresh';
clean.help = {'Start from fresh or keep existing priors'};
clean.labels = {'Fresh', 'Keep'};
clean.values = {1, 0};
clean.val = {1};

priorsmask  = cfg_files;
priorsmask.tag = 'priorsmask';
priorsmask.name = 'Priors file';
priorsmask.filter = '(.*\.gii$)|(.*\.mat$)|(.*\.nii(,\d+)?$)|(.*\.img(,\d+)?$)';
priorsmask.num = [0 1];
priorsmask.help = {'Select a mask or a mat file with priors or posterior.'};
priorsmask.val = {{''}};

funcpriors = cfg_branch;
funcpriors.tag = 'funcpriors';
funcpriors.name = 'Functionally defined priors';
funcpriors.help = {'Use posterior from previous inversion or restrict solutions to pre-specified VOIs'};
funcpriors.val  = {priorsmask, space,sparsify};

priortype = cfg_repeat;
priortype.tag = 'priortype';
priortype.name = 'Prior';
priortype.help = {'Specify the prior type to add'};
priortype.num  = [1 Inf];
priortype.values  = {randpatch,fixedpatch,funcpriors,popularpars};
priortype.val  = {randpatch};

%
% npost = cfg_entry;
% npost.tag = 'npost';
% npost.name = 'N posteriors';
% npost.strtype = 'n';
% npost.help = {'Nunber of posterior files to generate'};
% npost.val = {1};
%
% noisepost = cfg_entry;
% noisepost.tag = 'noisepost';
% noisepost.name = 'Percentage max noise';
% noisepost.strtype = 'r';
% noisepost.help = {'Add Gaussian noise with of a certain percentage of maximum power to posterior'};
% noisepost.val = {0};
%
%
% smoothpost = cfg_entry;
% smoothpost.tag = 'smoothpost';
% smoothpost.name = 'Number of smoothing iterations';
% smoothpost.strtype = 'r';
% smoothpost.help = {'Number of smoothing iterations (0 for no smoothing)'};
% smoothpost.val = {0};
%
% noisepost = cfg_entry;
% noisepost.tag = 'noisepost';
% noisepost.name = 'Percentage max noise';
% noisepost.strtype = 'r';
% noisepost.help = {'Add Gaussian noise with of a certain percentage of maximum power to posterior'};
% noisepost.val = {0};
%
% threshpost = cfg_entry;
% threshpost.tag = 'threshpost';
% threshpost.name = 'Percentage of maximum power at at which to threshold';
% threshpost.strtype = 'r';
% threshpost.help = {'Save posterior with full (0) or thresholded power distribution'};
% threshpost.val = {0};
%
% nmaxpost = cfg_entry;
% nmaxpost.tag = 'nmaxpost';
% nmaxpost.name = 'Number of maxima to keep as patch centers';
% nmaxpost.strtype = 'r';
% nmaxpost.help = {'Generate another file containing patch centres in Ip structure (0 make no file)'};
% nmaxpost.val = {0};
%
%
% postname = cfg_entry;
% postname.tag = 'postname';
% postname.name = 'Prefix for posterior';
% postname.strtype = 's';
% postname.num = [1 Inf];
% postname.val = {'postset1'};
% postname.help = {'Prefix for posterior directory'};
%

priors = cfg_exbranch;
priors.tag = 'priors';
priors.name = 'Inversion priors';
priors.val = {D, val,sensorlevel,priortype,priorname,clean};
priors.help = {'Set priors for source reconstruction'};
priors.prog = @add_priors;
priors.vout = @vout_priors;

%==========================================================================
function  out = add_priors(job)

val=job.val;
D = spm_eeg_load(job.D{val});

if ~isfield(D.inv{job.val}.inverse,'L')
    error('Need to set up spatial modes first');
else
    UL=D.inv{val}.inverse.L; %% dimension reduced lead field matrix
end

[a1,b1,c1]=fileparts(D.fname);
priordir=[D.path filesep job.priorname '_' b1];

if exist(priordir,'dir') ~= 7
    fprintf('Making directory for priors: %s\n',priordir);
    mkdir(priordir);
else
    mesh=D.inv{job.val}.mesh.tess_mni;
    [priorfiles] = spm_select('FPListRec',priordir,'.*\.mat$');

    count=size(priorfiles,1);
    fprintf('Found %d existing priors in this directory\n',count);

    if job.clean
        fprintf('Deleting existing priors\n');
        for f=1:count
            delete(deblank(priorfiles(f,:)));
        end
    else
        fprintf('Keeping existing priors\n');
    end

end

mesh=D.inv{val}.forward.mesh; %% assume this is in metres

Ns=size(mesh.vert,1);

% Set up sensor level prior
Qe{1}=eye(size(UL,1)).*job.sensorlevel;
Qp={}; %% source level covariance or svd components

% add in functional (from other modalities or experiment) hypotheses

for i=1:numel(job.priortype) %% move through different prior specifications
    allIp=[];
    Qp=[];

    % RANDOM PATCH CENTRES
    if isfield(job.priortype{i},'randpatch')

        base=job.priortype{i}.randpatch;

        Np = base.npatches;
        Npatchiter = base.niter;
        probfile=base.probfile;

        if isempty(probfile)
            prob=ones(Ns,1)./Ns;
        else

            probmesh=gifti(probfile);

            data=double(probmesh.cdata);
            M0.vertices=probmesh.vertices;
            M0.faces=probmesh.faces;

            if size(probmesh.vertices,1)==size(mesh.vert,1)
            %    fprintf('NB smoothing prior with 20 iterations\n');
                %prob=spm_mesh_smooth(M0,data',20);
                prob=data;
            else
                error('Gifti mesh size does not match');
            end
        end
        if sum(prob)>1.01
            warning('Rescaling bias file to unit sum');
        end

        prob=prob./sum(prob);
        Ip=[];

        fprintf('Using %d iterations of %d random patches\n',Npatchiter,Np);
        allIp=sparse(Npatchiter,Np);

        %rng('shuffle');

        dum=sparse(1,Ns);
        figure;
        hold on;

        for j=1:Npatchiter

            for k=1:Np
                pval=max(prob)*(k-1)./Np;
                sind=find(prob>=pval);
                randind=randperm(length(sind));
                allIp(j,k)=sind(randind(1));
                dum(allIp(j,k))=dum(allIp(j,k))+1;
            end
            %end; %% IF
            plot(1:length(prob),dum,'.');  drawnow;
        end % for j
        fprintf('Sampled %3.2f percent of cortex (%d vertices) \n',100*length(find(dum>0))./length(dum),length(find(dum>0)));

    end % random patch

    % FIXED PATCH CENTRES FROM FILE
    if isfield(job.priortype{i},'fixedpatch')

        base=job.priortype{i}.fixedpatch;
        try
            dum=load(char(base.fixedfile));
        catch
            error('failed to load fixedpatch file');
        end
        if ~isfield(dum,'Ip')
            error('Need to have patch indices in structure Ip');
        end
        allIp=dum.Ip;
        rows=base.fixedrows;
        lastrow=min(rows(2),size(allIp,1));
        rowind=[rows(1):lastrow];
        allIp=allIp(rowind,:);
    end % fixed patch

    %%%%%%%%%%%%%%% FIXED OR RANDOM PATCHES
    if isfield(job.priortype{i},'fixedpatch') || isfield(job.priortype{i},'randpatch')

        % Set up patches on cortical surface determined by indices of allIp
        % one prior files per row of Ip

        [Qp,Qe,allpriornames]=spm_eeg_invert_setuppatches(allIp,mesh,base,priordir,Qe,UL);

    end % if fixedpatch or randpatch

    % figure;
    %
    % if job.nmaxpost>0,
    %     postIpname=[postdir filesep sprintf('post_Ip.mat')];
    %     Ip=zeros(job.npost,job.nmaxpost);
    % end;
    %
    % M0=[];M0.vertices=mesh.vert; M0.faces=mesh.face;
    % for f=1:job.npost,
    %
    %     subplot(job.npost,1,f);
    %     plot(sourcevar,'k');hold on;
    %     postname=[postdir filesep sprintf('post%d.mat',f)];
    %     fprintf('adding %f percent noise to posterior\n',job.noisepost);
    %     nsourcevar=sourcevar+max(sourcevar)*job.noisepost*randn(size(sourcevar))./100;
    %     fprintf('Smoothing by %d iterations\n',job.smoothpost);
    %     nsourcevar = spm_mesh_smooth(M0, nsourcevar',job.smoothpost);
    %
    %     fprintf('Thresholding at %f percent of max power\n',job.threshpost);
    %     thresh=max(nsourcevar)*job.threshpost./100;
    %     plot([0,length(sourcevar)],[thresh thresh],':g');
    %     nsourcevar(nsourcevar<thresh)=0;
    %     plotind=find(nsourcevar);
    %     plot(plotind,nsourcevar(plotind),'ro');
    %     legend('orig','thresh','thresh noisy');
    %     nsourcevar=sparse(nsourcevar);
    %     Qp{1}=diag(nsourcevar);
    %     fprintf('SAving %s\n',postname);
    %     F=-Inf;
    %     save(postname,'Qp','Qe','UL','F','-v7.3');
    %     %%
    %     fprintf('\n Getting %d local maxima ',job.nmaxpost);
    %     Ip1 = spm_mesh_get_lm(M0,nsourcevar); %% get local maxima
    %
    %     [vals,ind]=sort(nsourcevar(Ip1),'descend');
    %     Ip(f,:)=Ip1(ind(1:job.nmaxpost));
    %
    %
    %
    % end; % for f
    % if job.nmaxpost>0,
    %     fprintf('Saving %s\n',postIpname);
    %     save(postIpname,'Ip'); % =[postdir filesep sprintf('post_Ip.mat')];
    %     %Ip=zeros(job.npost,job.nmaxpost);
    % end;

    %%%%%%%%%%% FUNCTIONALLY DEFINED PRIORS

    if isfield(job.priortype{i},'funcpriors')
        base=job.priortype{i}.funcpriors;
        P = char(base.priorsmask);
        priors{i}.priortype='funcpriors';
        priors{i}.priorfname=P;

        fprintf('Using priors from file %s\n',P);

        P=priors{i}.priorfname;
        [p,f,e] = fileparts(P);
        switch lower(e)
            case '.mat'
                a=load(P);
                if isfield(a,'Qp')

                    disp('Replacing sensor and source posterior distribution');
                    for f=1:numel(a.Qp)
                        Qp{end+1}=a.Qp{f};
                    end

                    % now need to compress the original Qe (say 275 channels) into reduced channel form
                    U=D.inv{val}.inverse.U{1}; %% get spatial modes definition
                    disp('NOT SURE IF WE NEED TO RESCALE FOR CHAN REDUCTION');
                    scaleU=size(a.Qe,1)./size(U,1); % white noise power pwer channel reduce with number of channels

                    Qe{1}=U*a.Qe*U'./scaleU;
                    fprintf('\n Reducing white noise variance from %3.2f fT^2 on each of %d channels to %3.2f fT^2 on each of %d channels\n',...
                        mean(diag(a.Qe)),size(a.Qe,1),mean(diag(Qe{1})),size(U,1));

                else
                    disp('loading just source level posterior');
                    Qp{end+1}=a.pQ;
                end

            case '.gii'
                g = gifti(P);

                for i=1:size(g.cdata,2)
                    Qp{end+1} = double(g.cdata(:,i));
                end

            case {'.img', '.nii'}
                S.D = D;
                S.fmri = P;
                S.space = job.isstandard.custom.priors.space;
                D = spm_eeg_inv_fmripriors(S);
                inverse.fmri = D.inv{D.val}.inverse.fmri;
                a=load(inverse.fmri.priors);

                Qp{end+1} = a.pQ;
            otherwise
                error('Unknown file type.');
        end % switch

        idnum=randi(1e6);

        priorfname=[priordir filesep sprintf('prior%d.mat',idnum)];
        fprintf('Saving %s\n',priorfname);
        F=[]; % no associated free energy value
        save(priorfname,'Qp','Qe','UL','F');
    end % func priors

    % POPULAR PRIORS

    if isfield(job.priortype{i},'popularpars')
        base=job.priortype{i}.popularpars;

        fprintf('Adding prior covariance matrix for %s algorithm\n',base.popular);

        if strcmp(base.popular,'sparseEBB')

            % calculate power on cortical surface
            % using beamformer assumptions
            AYYA=D.inv{val}.inverse.AYYA;

            InvCov = spm_inv(AYYA);

            allsource = sparse(Ns,1);
            Sourcepower = sparse(Ns,1);

            for bk = 1:Ns
                normpower = 1/(UL(:,bk)'*UL(:,bk));
                Sourcepower(bk) = 1/(UL(:,bk)'*InvCov*UL(:,bk));
                % normalise power by unit noise through lead fields
                allsource(bk) = Sourcepower(bk)./normpower;
            end

            % now get local maxima on mesh
            M.vertices=mesh.vert;
            M.faces=mesh.face;

            Ip = spm_mesh_get_lm(M,allsource); %% get local maxima
            figure;
            plot(allsource);
            legend('sparse EBB');
            hold on;
            maxBFpatch=100;
            fprintf('Limiting to a max of %d peaks\n',maxBFpatch);

            [vals,ind]=sort(allsource(Ip),'descend');
            Ip=Ip(ind(1:maxBFpatch));
            plot(Ip,allsource(Ip),'ro');

            [Qp,Qe,allpriornames]=spm_eeg_invert_setuppatches(Ip,mesh,base,priordir,Qe{1},UL);

        end % sparse EBB

        if strcmp(base.popular,'EBB')

            % calculate power on cortical surface
            % using beamformer assumptions
            AYYA=D.inv{val}.inverse.AYYA;

            InvCov = spm_inv(AYYA);

            allsource = sparse(Ns,1);
            Sourcepower = sparse(Ns,1);

            for bk = 1:Ns
                normpower = 1/(UL(:,bk)'*UL(:,bk));
                Sourcepower(bk) = 1/(UL(:,bk)'*InvCov*UL(:,bk));
                % normalise power by unit noise through lead fields
                allsource(bk) = Sourcepower(bk)./normpower;
            end

            figure;
            plot(allsource);
            legend('EBB');
            hold on;
            % now get local maxima on mesh
            M.vertices=mesh.vert;
            M.faces=mesh.face;
            warning('Smoothing not used');
            Qp{1}=diag(allsource);

            idnum=randi(1e6);
            priorfname=[priordir filesep sprintf('prior%d.mat',idnum)];
            F=[];
            save(priorfname,'Qp','Qe','UL','F');

        end % EBB

        if strcmp(base.popular,'COH')

            % create minimum norm prior
            %------------------------------------------------------------------
            Qp{1}  = speye(Ns,Ns);
            warning('SQUARE SMOOTHING KERNEL OF FWHM MM');

            M=[];M.faces=mesh.face;M.vertices=mesh.vert;
            dist=spm_mesh_distmtx(M,1);
            ind2=find(dist>0);
            ind=find(dist<base.FWHMmm./1000);
            ind=intersect(ind,ind2);
            Qp{1}(ind)=1;

            disp('Set up COH prior');

            idnum=randi(1e6);
            priorfname=[priordir filesep sprintf('prior%d.mat',idnum)];
            fprintf('Saving %s\n',priorfname);
            F=[]; % no associated free energy value

            save(priorfname,'Qp','Qe','UL','F');
        end % COH

        if strcmp(base.popular,'IID')

            % create minimum norm prior
            %------------------------------------------------------------------
            Qp{1}  = speye(Ns,Ns);

            disp('Set up IID prior');

            idnum=randi(1e6);

            priorfname=[priordir filesep sprintf('prior%d.mat',idnum)];
            fprintf('Saving %s\n',priorfname);
            F=[]; % no associated free energy value
            save(priorfname,'Qp','Qe','UL','F');
        end % IID

    end % popular priors

end % for i (looping over priors)

% NOW JUST WANT TO PLOT OUT THE CURRENT DENSITY PRIORS FOR A RANDOM SET OF PATCHES
% THIS SHOULD ALL GO IN A PREVIEW FUNCTION

[LCpL,Q,sumLCpL,QE,Cy,M,Cp,Cq,Lq]=spm_eeg_assemble_priors(UL,Qp,Qe);

figure;
mesh=D.inv{job.val}.mesh.tess_mni;
M0.vertices=mesh.vert;
M0.faces=mesh.face;
sourcepwr=diag(Cp);

spm_mip(sourcepwr,mesh.vert',6);
title('Prior');
colorbar;

CHECKSMOOTH=0;
if CHECKSMOOTH
    figure;
    r=randperm(length(Qp));

    sigma2=[];
    for f=1:min(length(r),4)
        if isfield(Qp{r(f)},'q')
            q=Qp{r(f)}.q;
            [val,ind]=max(q);

            order=1; %1 distance in mm, 0 distance in edges
            dist = spm_mesh_geodesic(M,ind-1,order);
            distind=intersect(find(dist<0.020),find(q>0)); % less than 2cm away

            y=spm_vec(log(q(distind)'));
            x=spm_vec(-dist(distind).^2);
            B = regress(y,[x ones(size(x,1),1)] );
            sigma2(f)=1./(2*B(1)); %% sd squared of gaussian
            plot(dist(distind).*1000,q(distind),'r.',dist(distind).*1000,max(q).*exp(-(dist(distind).^2)/(2*sigma2(f))),'go'); %% plot in mm
        end

    end

    if ~isempty(sigma2)
        title(sprintf('Prior Profile of prior amplitude across cortex'));
        meanfwhm=mean(sqrt(sigma2))*2.355*1000; %% in mm
        stdfwhm=std(sqrt(sigma2))*2.355*1000;
        legend(sprintf('Mean FWHM=%3.2f mm,std FWHM= %3.2fmm',meanfwhm,stdfwhm));
        drawnow;
    end

end %% if CHECKSMOOTH

out.D = job.D;

%==========================================================================
function dep = vout_priors(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'M/EEG dataset(s) after imaging source reconstruction';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'filter','mat'}});

% function [q,dist,useind]= gauss_patch(M,i,FWHM,q)
%
% order=1;
%
%
% sigma=FWHM./2.355;
% sigma2=sigma^2;
%
% d=spm_mesh_geodesic(M,i-1,order);
% useind=find(d<FWHM*2);
% dist=d(useind);
% q(useind)=exp(-(dist.^2)/(2*sigma2));

author: spm


text:
% Fig. 6.57a   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%
% response comparison of continuous and digital control step.

clear all;
close all;

sim('fig6_57c')
r=[1 1];   %reference input
t=[0 5];
plot(t,r,'r')
hold on
plot(ycd(:,1),ycd(:,2))
plot(ycd(:,1),ycd(:,3),'m:')
title('Figure 6.57(a) Step Responses of Digital and Continuous Controllers')
ylabel('y')
xlabel('Time (sec)')
text(1.1,1.06, '\leftarrow continuous  controller')
text(.7,1.23, '\leftarrow digital  controller')
grid
hold off

author: Sable

repo: mcbench-benchmarks


text:
function Value = bst_units_ui(Units, Value)
% BST_UNITS_IU: Scale values to international units (meters)
%
% USAGE:   Value = bst_units_ui(Units, Value)   % Scale values in input
%         Factor = bst_units_ui(Units)          % Returns scaling factor

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

% Supported units
switch strtrim(str_remove_spec_chars(lower(Units)))
    case 'mm'
        factor = 0.001;
    case 'cm'
        factor = 0.01;
    case 'm'
        factor = 1;
    otherwise
        disp(['BST> Warning: Unknown units "' Units '"']);
        factor = 1;
end
% Return scaling factor
if (nargin < 2)
    Value = factor;
% Apply scaling factor
elseif (factor ~= 1)
    Value = factor .* Value;
end

author: brainstorm-tools


text:
function var = kWaveGPUsingle(var)
%KWAVEGPUSINGLE   GPU datacasting wrapping for GPUmat library
%
% USAGE:
%       var = kWaveGPUsingle(var)
%
% ABOUT:
%       author: Bradley E. Treeby
%       date: 21st January 2010
%       last update: 21st January 2010
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% only cast non scalar variables
if ~isscalar(var)
    var = GPUsingle(var);
end

repo: Aero-matlab


text:
function [Y,FS] = spm_DEM_play_song(qU,T);
% displays the song-bird images specified by the states in qU
% FORMAT [Y,FS] = spm_DEM_play_song(qU,T);
%
% qU   - conditional moments of states (see spm_DEM)
% T    - number of seconds over which to play the sound
%
% Y    - sound image
% FS   - sampling rate (Hz)
%
% A button press on the spectrogram will play the song
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_DEM_play_song.m 6044 2014-06-14 10:22:46Z karl $

% load frequency modes
%--------------------------------------------------------------------------
try, T;           catch, T = 2;  end
try, v = qU.v{1}; catch, v = qU; end

% create sound image
%==========================================================================
[Nm m] = size(v);

% frequencies
%--------------------------------------------------------------------------
Hf  = 5000;                                % upper frequency (Hz)
Lf  = 2500;                                % lower frequency (Hz)
Nf  = 64;                                  % number of frequency bin
Hz  = linspace(Lf,Hf,64)';                 % frequencies
FS  = 2*Hz(end);                           % sampling rate (Hz)
k   = Hz/Hz(1);                            % cycles per window
n   = FS/Hz(1);                            % window length
N   = FS*T;                                % number of sonogram bins
R   = fix(N/m);                            % interpolation factor
N   = R*m;
pst = (1:N)/FS;                            % peristimulus time
sf  = 2*64^2;                              % dispersion of frequencies

% resample temporal modes
%--------------------------------------------------------------------------
for i = 1:Nm
    V(i,:) = spm_interp(v(i,:),R);
end

% create sonogram sound
%--------------------------------------------------------------------------
b     = V(1,:);                            % amplitude modulation
f     = V(2,:);                            % frequency modulation
b     = abs(b);
b     = b/max(b);
b     = tanh((b - 1/2)*6) + 1;
f     = 32*f + Lf;
S     = sparse(Nf,N);
for i = 1:N
    s      = b(i)*exp(-(Hz - f(i)).^2/sf);
    s      = sparse(s.*(s > exp(-4)));
    S(:,i) = s;
end

% inverse Fourier transform
%--------------------------------------------------------------------------
Y   = spm_iwft(S,k,n);
Y   = Y/max(Y);

% Graphics
%==========================================================================
colormap('pink')
imagesc(pst,Hz,abs(S))
axis xy
xlabel('time (sec)')
ylabel('Frequency (Hz)')

% set sound data
%--------------------------------------------------------------------------
h      = get(gca,'Children');
set(h(1),'Userdata',{Y,FS})
set(h(1),'ButtonDownFcn','spm_DEM_ButtonDownFcn')


text:
function amplitudes = LoadSpikeAmplitudes(filename,nChannels,nSamples,peak,rate)

%LoadSpikeAmplitudes - Read spike times, cluster IDs and amplitudes from disk.
%
%  USAGE
%
%    amplitudes = LoadSpikeAmplitudes(filename,nChannels,nSamples)
%
%    filename            spike waveform file name
%    nChannels           number of channels in electrode group
%    nSamples            number of samples per waveform
%    peak                peak sample (1..N)
%    rate                sampling rate
%
%  SEE
%
%    See also GetSpikeAmplitudes, LoadSpikeWaveforms.

% Copyright (C) 2004-2012 by Michaël Zugaro
%
% 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.

if nargin < 3,
  error('Incorrect number of parameters (type ''help LoadSpikeAmplitudes'' for details).');
end

[path,basename,extension] = fileparts(filename);
if isempty(path), path = '.'; end

electrodeGroupStr = extension(2:end);
electrodeGroup = str2num(electrodeGroupStr);
[~,basename] = fileparts(basename);

% Times and cluster IDs
times = LoadSpikeTimes(filename,rate);

% Load .spk file
filename = [path '/' basename '.spk.' electrodeGroupStr];
if ~exist(filename),
	error(['File ''' filename ''' not found.']);
end
amp = LoadBinary(filename,'offset',peak-1,'skip',nSamples-1,'nChannels',nChannels);

amplitudes = [times amp];

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function i = dyad(j)
% dyad -- Index entire j-th dyad of 1-d wavelet xform
%  Usage
%    ix = dyad(j);
%  Inputs
%    j     integer
%  Outputs
%    ix    list of all indices of wavelet coeffts at j-th level
%
    i = (2^(j)+1):(2^(j+1)) ;

%
% Copyright (c) 1993. David L. Donoho
%

%
%  Part of Wavelab Version 850
%  Built Tue Jan  3 13:20:40 EST 2006
%  This is Copyrighted Material
%  For Copying permissions see COPYING.m
%  Comments? e-mail wavelab@stat.stanford.edu

author: BUNPC


text:
function [ cur_group ] = kmeans_diy( input_data, k )
n = size(input_data, 1);
d = size(input_data, 2);
% initialise means from random samples
means = input_data(randsample(n, k), :);
% main loop
newmeans = zeros(k,d);
scores = zeros(n, k);
change = Inf;
while change > eps
    % assign labels to data points according to distances to the means
    for j = 1:k
        diff = bsxfun(@minus, input_data, means(j,:));
        scores(:,j) = sqrt(sum(diff.^2,2));
    end
    [~,cur_group] = min(scores, [], 2);
    % recompute means based on obtained labelings
    for j = 1:k
        newmeans(j,:) = mean(input_data(cur_group == j,:));
    end
    % stopping criteria
    change = sqrt(sum(sum((newmeans-means).^2, 1), 2))/k;
    means = newmeans;
end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function value = s_eqi ( s1, s2 )

%*****************************************************************************80
%
%% S_EQI is a case insensitive comparison of two strings for equality.
%
%  Example:
%
%    S_EQI ( 'Anjana', 'ANJANA' ) is TRUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S1, S2, the strings to compare.
%
%    Output, logical VALUE, is TRUE if the strings are equal.
%
  len1 = length ( s1 );
  len2 = length ( s2 );
  lenc = min ( len1, len2 );

  for i = 1 : lenc

    c1 = ch_cap ( s1(i) );
    c2 = ch_cap ( s2(i) );

    if ( c1 ~= c2 )
      value = 0;
      return
    end

  end

  for i = lenc + 1 : len1
    if ( s1(i) ~= ' ' )
      value = 0;
      return
    end
  end

  for i = lenc + 1 : len2
    if ( s2(i) ~= ' ' )
      value = 0;
      return
    end
  end

  value = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = sph2cart2d(theta, phi, rho)
%SPH2CART2D Convert spherical coordinates to cartesian coordinates in degrees
%
%   C = SPH2CART2(THETA, PHI, RHO)
%   C = SPH2CART2(THETA, PHI)       (assume rho = 1)
%   C = SPH2CART2(S)
%   [X, Y, Z] = SPH2CART2(THETA, PHI, RHO);
%
%   S = [phi theta rho] (spherical coordinate).
%   C = [X Y Z]  (cartesian coordinate)
%
%   The following convention is used:
%   THETA is the colatitude, in degrees, 0 for north pole, +180 degrees for
%   south pole, +90 degrees for points with z=0.
%   PHI is the azimuth, in degrees, defined as matlab cart2sph: angle from
%   Ox axis, counted counter-clockwise.
%   RHO is the distance of the point to the origin.
%   Discussion on choice for convention can be found at:
%   http://www.physics.oregonstate.edu/bridge/papers/spherical.pdf
%
%   See also:
%   angles3d, cart2sph2d, sph2cart2
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-06-29,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% Process input arguments
if nargin == 1
    phi     = theta(:, 2);
    if size(theta, 2) > 2
        rho = theta(:, 3);
    else
        rho = ones(size(phi));
    end
    theta   = theta(:, 1);

elseif nargin == 2
    rho     = ones(size(theta));

end

% conversion
rz = rho .* sind(theta);
x  = rz  .* cosd(phi);
y  = rz  .* sind(phi);
z  = rho .* cosd(theta);

% Process output arguments
if nargout == 1 || nargout == 0
    varargout{1} = [x, y, z];

else
    varargout{1} = x;
    varargout{2} = y;
    varargout{3} = z;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ n_data, n, x, fx ] = legendre_poly_values ( n_data )

%*****************************************************************************80
%
%% LEGENDRE_POLY_VALUES returns values of the Legendre polynomials.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      LegendreP [ n, x ]
%
%  Differential equation:
%
%    (1-X*X) * P(N,X)'' - 2 * X * P(N,X)' + N * (N+1) = 0
%
%  First terms:
%
%    P( 0,X) =       1
%    P( 1,X) =       1 X
%    P( 2,X) =  (    3 X^2 -       1)/2
%    P( 3,X) =  (    5 X^3 -     3 X)/2
%    P( 4,X) =  (   35 X^4 -    30 X^2 +     3)/8
%    P( 5,X) =  (   63 X^5 -    70 X^3 +    15 X)/8
%    P( 6,X) =  (  231 X^6 -   315 X^4 +   105 X^2 -     5)/16
%    P( 7,X) =  (  429 X^7 -   693 X^5 +   315 X^3 -    35 X)/16
%    P( 8,X) =  ( 6435 X^8 - 12012 X^6 +  6930 X^4 -  1260 X^2 +   35)/128
%    P( 9,X) =  (12155 X^9 - 25740 X^7 + 18018 X^5 -  4620 X^3 +  315 X)/128
%    P(10,X) =  (46189 X^10-109395 X^8 + 90090 X^6 - 30030 X^4 + 3465 X^2-63 ) /256
%
%  Recursion:
%
%    P(0,X) = 1
%    P(1,X) = X
%    P(N,X) = ( (2*N-1)*X*P(N-1,X)-(N-1)*P(N-2,X) ) / N
%
%    P'(0,X) = 0
%    P'(1,X) = 1
%    P'(N,X) = ( (2*N-1)*(P(N-1,X)+X*P'(N-1,X)-(N-1)*P'(N-2,X) ) / N
%
%  Formula:
%
%    P(N,X) = (1/2**N) * sum ( 0 <= M <= N/2 ) C(N,M) C(2N-2M,N) X^(N-2*M)
%
%  Orthogonality:
%
%    Integral ( -1 <= X <= 1 ) P(I,X) * P(J,X) dX
%      = 0 if I =/= J
%      = 2 / ( 2*I+1 ) if I = J.
%
%  Approximation:
%
%    A function F(X) defined on [-1,1] may be approximated by the series
%
%      C0*P(0,X) + C1*P(1,X) + ... + CN*P(N,X)
%
%    where
%
%      C(I) = (2*I+1)/(2) * Integral ( -1 <= X <= 1 ) F(X) P(I,X) dx.
%
%  Special values:
%
%    P(N,1) = 1.
%    P(N,-1) = (-1)**N.
%    | P(N,X) | <= 1 in [-1,1].
%
%    P(N,0,X) = P(N,X), that is, for M=0, the associated Legendre
%    function of the first kind and order N equals the Legendre polynomial
%    of the first kind and order N.
%
%    The N zeroes of P(N,X) are the abscissas used for Gauss-Legendre
%    quadrature of the integral of a function F(X) with weight function 1
%    over the interval [-1,1].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the order of the function.
%
%    Output, real X, the point where the function is evaluated.
%
%    Output, real FX, the value of the function.
%
  n_max = 22;

  fx_vec = [ ...
      0.1000000000000000E+01, ...
      0.2500000000000000E+00, ...
     -0.4062500000000000E+00, ...
     -0.3359375000000000E+00, ...
      0.1577148437500000E+00, ...
      0.3397216796875000E+00, ...
      0.2427673339843750E-01, ...
     -0.2799186706542969E+00, ...
     -0.1524540185928345E+00, ...
      0.1768244206905365E+00, ...
      0.2212002165615559E+00, ...
      0.0000000000000000E+00, ...
     -0.1475000000000000E+00, ...
     -0.2800000000000000E+00, ...
     -0.3825000000000000E+00, ...
     -0.4400000000000000E+00, ...
     -0.4375000000000000E+00, ...
     -0.3600000000000000E+00, ...
     -0.1925000000000000E+00, ...
      0.8000000000000000E-01, ...
      0.4725000000000000E+00, ...
      0.1000000000000000E+01 ];

  n_vec = [ ...
     0,  1,  2, ...
     3,  4,  5, ...
     6,  7,  8, ...
     9, 10,  3, ...
     3,  3,  3, ...
     3,  3,  3, ...
     3,  3,  3, ...
     3 ];

  x_vec = [ ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.25E+00, ...
     0.00E+00, ...
     0.10E+00, ...
     0.20E+00, ...
     0.30E+00, ...
     0.40E+00, ...
     0.50E+00, ...
     0.60E+00, ...
     0.70E+00, ...
     0.80E+00, ...
     0.90E+00, ...
     1.00E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    x = 0.0;
    fx = 0.0;
  else
    n = n_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function [biomassValues, targetLowerBounds, targetUpperBounds, plottedReactions] = multiProductionEnvelope(model, deletions, biomassRxn, geneDelFlag, nPts, plotAllFlag, plotTools)
% Calculates the byproduct secretion envelopes for
% every product (excreted metabolites with 1 or more Carbons)
%
% USAGE:
%
%    [biomassValues, targetValues] = multiProductionEnvelope(model, deletions, biomassRxn, geneDelFlag, nPts, plotTools)
%
% INPUT:
%    model:            COBRA model structure
%
% OPTIONAL INPUT:
%    deletions:        List of reaction or gene deletions (empty if wild type)
%                      (Default = {})
%    biomassRxn:       Biomass `rxn` name (Default = whatever is defined in model)
%    geneDelFlag:      Perform gene and not reaction deletions (Default = false)
%    nPts:             Number of points in the plot (Default = 20)
%    plotTools:        boolean (default = false) - add tools for editing the figure and its properties
%    plotAllFlag:      plot all envelopes, even ones that are not growth
%                      coupled (Default = false)
%    plotTools:        boolean (default = false) - add tools for editing the figure and its properties
%
% OUTPUT:
%    biomassValues:         Biomass values for plotting
%    targetUpperBounds:     Target upper bounds for plotting (
%                           biomassvalues x reactions)
%    targetLowerBounds:     Target lower bounds for plotting (
%                           biomassvalues x reactions)
%    plottedReactions:      Reactions that led to relevant side product
%
% .. Author: - Jeff Orth 8/16/07

if (nargin < 2)
    deletions = {};
end
if (nargin < 3)
    % Biomass flux
    biomassRxn = model.rxns(model.c==1);
end
if (nargin < 4)
    % Gene or rxn deletions
    geneDelFlag = false;
end
if (nargin < 5)
    nPts = 20;
end
if ~exist('plotTools','var')
    plotTools = false;
end

if ~exist('plotAllFlag','var')
    plotAllFlag = false;
end

% Create model with deletions
if (length(deletions) > 0)
    if (geneDelFlag)
        model = deleteModelGenes(model,deletions);
    else
        model = changeRxnBounds(model,deletions,zeros(size(deletions)),'b');
    end
end

%get all C exchange reactions
excRxns = model.rxns(findExcRxns(model,false,false));
CRxns = findCarbonRxns(model,1);
CExcRxns = intersect(excRxns,CRxns);
substrateIDs = find(model.lb(findRxnIDs(model,CExcRxns))<0);
%remove the substrate reactions
for i = 1:length(substrateIDs)
    j = substrateIDs(i);
    if j == 1
        CExcRxns = CExcRxns(2:length(CExcRxns));
    elseif j == length(CExcRxns)
        CExcRxns = CExcRxns(1:length(CExcRxns)-1);
    else
        CExcRxns = cat(1,CExcRxns(1:j-i),CExcRxns(j-i+2:length(CExcRxns)));
    end
end

% Run FBA to get upper bound for biomass
model = changeObjective(model,biomassRxn,1);
solMax = optimizeCbModel(model,'max');
solMin = optimizeCbModel(model,'min');

% Create biomass range vector
biomassValues = linspace(solMin.f,solMax.f,nPts);

plottedReactions = {};
targetUpperBounds = zeros(0,0);
targetLowerBounds = zeros(0,0);
% Max/min for target production
for i = 1:length(CExcRxns)
    model = changeObjective(model,CExcRxns(i),1);
    model2 = changeRxnBounds(model,biomassRxn,max(biomassValues),'b');
    fbasol2 = optimizeCbModel(model2,'max');
    maxRate = fbasol2.f; %find max production at max growth rate
    if (plotAllFlag)||(maxRate > getCobraSolverParams('LP','feasTol')) %Change to plot everything that is above the detection limit
        plottedReactions = [plottedReactions,CExcRxns(i)];
        targetUpperBounds(end+1,1) = 0;
        targetLowerBounds(end+1,1) = 0;
        for j = 1:length(biomassValues)
            model = changeRxnBounds(model,biomassRxn,biomassValues(j),'b');
            sol = optimizeCbModel(model,'max');
            if (sol.stat == 1)
                targetUpperBounds(end,j) = sol.f;
            elseif (sol.stat == 2)
                targetUpperBounds(end,j) = Inf;
            else
                targetUpperBounds(end,j) = NaN;
            end
            sol = optimizeCbModel(model,'min');
            if (sol.stat == 1)
                targetLowerBounds(end,j) = sol.f;
            elseif (sol.stat == 2)
                  targetUpperBounds(end,j) = Inf;
            else
                targetLowerBounds(end,j) = NaN;
            end
        end
    end
end

% Plot results
colors = {'b','g','r','c','m','y','k'};
for i = 1:length(plottedReactions)
    plot([biomassValues fliplr(biomassValues)],[targetUpperBounds(i,:) fliplr(targetLowerBounds(i,:))],colors{mod(i-1,length(colors))+1},'LineWidth',2)
    axis tight;
    hold on;
end
hold off;
legend(plottedReactions);
legend off;
ylabel('Production Rate (mmol/gDW h)');
xlabel('Growth Rate (1/h)');
if plotTools
    plottools, plotbrowser('on'), figurepalette('hide'), propertyeditor('off');
end
biomassValues = biomassValues';
targetLowerBounds = targetLowerBounds';
targetUpperBounds = targetUpperBounds';

author: opencobra

repo: cobratoolbox


text:
% Minimize sidelobe level of an array with arbitrary 2-D geometry
% "Convex optimization examples" lecture notes (EE364) by S. Boyd
% "Antenna array pattern synthesis via convex optimization"
% by H. Lebret and S. Boyd
% (figures are generated)
%
% Designs an antenna array such that:
% - it minimizes sidelobe level outside the beamwidth of the pattern
% - it has a unit sensitivity at some target direction
% - it has nulls (zero sensitivity) at specified direction(s) (optional)
%
% This is a convex problem (after sampling it can be formulated as an SOCP).
%
%   minimize   max |y(theta)|     for theta outside the beam
%       s.t.   y(theta_tar) = 1
%              y(theta_null) = 0  (optional)
%
% where y is the antenna array gain pattern (complex function) and
% variables are w (antenna array weights or shading coefficients).
% Gain pattern is a linear function of w: y(theta) = w'*a(theta)
% for some a(theta) describing antenna array configuration and specs.
%
% Written for CVX by Almir Mutapcic 02/02/06

% select array geometry
ARRAY_GEOMETRY = '2D_RANDOM';
% ARRAY_GEOMETRY = '1D_UNIFORM_LINE';
% ARRAY_GEOMETRY = '2D_UNIFORM_LATTICE';

% select if the optimal array pattern should enforce nulls or not
HAS_NULLS = 0; % HAS_NULLS = 1;

%********************************************************************
% problem specs
%********************************************************************
lambda = 1;           % wavelength
theta_tar = 60;       % target direction (should be an integer -- discretization)
half_beamwidth = 10;  % half beamwidth around the target direction

% angles where we want nulls (optional)
if HAS_NULLS
  theta_nulls = [95 110 120 140 225];
end

%********************************************************************
% random array of n antenna elements
%********************************************************************
if strcmp( ARRAY_GEOMETRY, '2D_RANDOM' )
  % set random seed to repeat experiments
  rand('state',0);

  % (uniformly distributed on [0,L]-by-[0,L] square)
  n = 40;
  L = 5;
  loc = L*rand(n,2);
  angleRange = 360;

%********************************************************************
% uniform 1D array with n elements with inter-element spacing d
%********************************************************************
elseif strcmp( ARRAY_GEOMETRY, '1D_UNIFORM_LINE' )
  % (unifrom array on a line)
  n = 30;
  d = 0.45*lambda;
  loc = [d*[0:n-1]' zeros(n,1)];
  angleRange = 180;

%********************************************************************
% uniform 2D array with m-by-m element with d spacing
%********************************************************************
elseif strcmp( ARRAY_GEOMETRY, '2D_UNIFORM_LATTICE' )
  m = 6; n = m^2;
  d = 0.45*lambda;

  loc = zeros(n,2);
  for x = 0:m-1
    for y = 0:m-1
      loc(m*y+x+1,:) = [x y];
    end
  end
  loc = loc*d;
  angleRange = 360;

else
  error('Undefined array geometry')
end

%********************************************************************
% construct optimization data
%********************************************************************
% build matrix A that relates w and y(theta), ie, y = A*w
theta = [1:angleRange]';
A = kron(cos(pi*theta/180), loc(:,1)') + kron(sin(pi*theta/180), loc(:,2)');
A = exp(2*pi*i/lambda*A);

% target constraint matrix
[diff_closest, ind_closest] = min( abs(theta - theta_tar) );
Atar = A(ind_closest,:);

% nulls constraint matrix
if HAS_NULLS
  Anull = []; ind_nulls = [];
  for k = 1:length(theta_nulls)
    [diff_closest, ind_closest] = min( abs(theta - theta_nulls(k)) );
    Anull = [Anull; A(ind_closest,:)];
    ind_nulls = [ind_nulls ind_closest];
  end
end

% stopband constraint matrix
ind = find(theta <= (theta_tar-half_beamwidth) | ...
           theta >= (theta_tar+half_beamwidth) );
if HAS_NULLS, ind = setdiff(ind,ind_nulls); end;
As = A(ind,:);

%********************************************************************
% optimization problem
%********************************************************************
cvx_begin
  variable w(n) complex
  minimize( max( abs(As*w) ) )
  subject to
    Atar*w == 1;   % target constraint
    if HAS_NULLS   % nulls constraints
      Anull*w == 0;
    end
cvx_end

% check if problem was successfully solved
disp(['Problem is ' cvx_status])
if ~strfind(cvx_status,'Solved')
  return
end

min_sidelobe_level = 20*log10( max(abs(As*w)) );
fprintf(1,'The minimum sidelobe level is %3.2f dB.\n\n',...
          min_sidelobe_level );

%********************************************************************
% plots
%********************************************************************
figure(1), clf
plot(loc(:,1),loc(:,2),'o')
title('Antenna locations')

% plot array pattern
if angleRange == 180,
    theta = [1:360]';
    A = [ A; -A ];
end
y = A*w;
figure(2), clf
ymin = floor(0.1*min_sidelobe_level)*10-10; ymax = 0;
plot([1:360], 20*log10(abs(y)), ...
     [theta_tar theta_tar],[ymin ymax],'r--',...
     [theta_tar+half_beamwidth theta_tar+half_beamwidth],[ymin ymax],'g--',...
     [theta_tar-half_beamwidth theta_tar-half_beamwidth],[ymin ymax],'g--');
if HAS_NULLS % add lines that represent null positions
  hold on;
  for k = 1:length(theta_nulls)
    plot([theta_nulls(k) theta_nulls(k)],[ymin ymax],'m--');
  end
  hold off;
end
xlabel('look angle'), ylabel('mag y(theta) in dB');
axis([0 360 ymin ymax]);

% polar plot
figure(3), clf
zerodB = -ymin;
dBY = 20*log10(abs(y)) + zerodB;
ind = find( dBY <= 0 ); dBY(ind) = 0;
plot(dBY.*cos(pi*theta/180), dBY.*sin(pi*theta/180), '-');
axis([-zerodB zerodB -zerodB zerodB]), axis('off'), axis('square')
hold on
plot(zerodB*cos(pi*theta/180),zerodB*sin(pi*theta/180),'k:') % 0 dB
plot( (min_sidelobe_level + zerodB)*cos(pi*theta/180), ...
      (min_sidelobe_level + zerodB)*sin(pi*theta/180),'k:')  % min level
text(-zerodB,0,'0 dB')
tt = text(-(min_sidelobe_level + zerodB),0,sprintf('%0.1f dB',min_sidelobe_level));
set(tt,'HorizontalAlignment','right');
theta_1 = theta_tar+half_beamwidth;
theta_2 = theta_tar-half_beamwidth;
plot([0 55*cos(theta_tar*pi/180)], [0 55*sin(theta_tar*pi/180)], 'k:')
plot([0 55*cos(theta_1*pi/180)], [0 55*sin(theta_1*pi/180)], 'k:')
plot([0 55*cos(theta_2*pi/180)], [0 55*sin(theta_2*pi/180)], 'k:')
if HAS_NULLS % add lines that represent null positions
  for k = 1:length(theta_nulls)
    plot([0 55*cos(theta_nulls(k)*pi/180)], ...
         [0 55*sin(theta_nulls(k)*pi/180)], 'k:')
  end
end
hold off

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [ mD ] = CalcDistanceMatrixARows( mX )
% ----------------------------------------------------------------------------------------------- %
% [ mD ] = CalcDistanceMatrixARows( mX )
%   Calculates the distance matrix for the input data. The distance matrix
%   is a symmetric matrix where 'mD(ii, jj) = dist(mX(ii, :), mX(jj, :));'.
%   This function uses the squared Euclidean Distance for the distance
%   metric.
% Input:
%   - mX            -   Data Matrix.
%                       Each data sample is a row of the matrix.
%                       Structure: Matrix (numVars x varDim).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% Output:
%   - mD            -   Distance Matrix.
%                       A symmetric matrix where 'mD(ii, jj) = dist(mX(ii,
%                       :), mX(jj, :));'.
%                       Structure: Matrix (numVars x numVars).
%                       Type: 'Single' / 'Double'.
%                       Range: [0, inf).
% References
%   1.  A
% Remarks:
%   1.  B
% TODO:
%   1.  C
% Release Notes:
%   -   1.0.000     01/01/2021  Royi Avital	RoyiAvital@yahoo.com
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

vSsqX   = sum(mX .^ 2, 2);
mD      = vSsqX.'+ vSsqX - (2 * (mX * mX.'));

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function pdf = circdensity(bc, bp, sig,varargin)
% Wrapped gaussian fit for circular, axial data in the range 0:pi
% Input:
%    bc  -  list of angles (= bincenters)
%    bp  -  list of values (= Bin population)
%    sig -  sigma of Gaussian used to smooth the distribution
%
% Options:
%    sum -  do not normalize to area = 1 (default) but sum of radii =1
%
% Output:
%    pdf - circular density function

% Todo: make it work for other functions as well
bc=reshape(bc,[],1);
bp=reshape(bp,[],1);
% check input
if bc(1)==bc(end)
    bc(end)=[];
    bp(end)=[];
    circflag = 1;
end

% make Gaussian kernel
G = (exp( -((bc - pi).^2)./ (2*(sig)^2) ));
% convolve kernel with data
pdf = ifft(fft(bp).*fft(G));
% normalize values
if check_option(varargin,'sum')
% to area. Is that useful?
pdf = pdf./sum(pdf);
else
pdf = pdf./sqrt(polyarea(cos(bc).*pdf,sin(bc).*pdf));
end

if circflag==1
 pdf(end+1) = pdf(1);
end

end

author: mtex-toolbox

repo: mtex

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


text:
function c = tanh(a)
% TANH for adiff objects.

c = adiff( tanh(a.x), rowmult(1./cosh(a.x).^2, a.dx), a.root);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function r8lib_test ( )

%*****************************************************************************80
%
%% R8LIB_TEST tests the R8LIB library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2015
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8LIB_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the R8LIB library.\n' );

  i4int_to_r8int_test ( );

  perm0_check_test ( );
  perm0_uniform_test ( );

  perm1_check_test ( );
  perm1_uniform_test ( );

  r8_abs_test ( );
  r8_acos_test ( );
  r8_acosh_test ( );
  r8_asinh_test ( );
  r8_atan_test ( );
  r8_atanh_test ( );
  r8_big_test ( );
  r8_cas_test ( );
  r8_ceiling_test ( );
  r8_choose_test ( );
  r8_cosd_test ( );
  r8_cotd_test ( );
  r8_cscd_test ( );
  r8_cube_root_test ( );
  r8_diff_test ( );
  r8_digit_test ( );
  r8_e_test ( );
  r8_epsilon_test ( );
  r8_epsilon_compute_test ( );
  r8_factorial_test ( );
  r8_factorial2_test ( );
  r8_fall_test ( );
  r8_fractional_test ( );
  r8_gamma_test ( );
  r8_gamma_log_test ( );
  r8_huge_test ( );
  r8_log_2_test ( );
  r8_log_b_test ( );
  r8_mant_test ( );
  r8_max_test ( );
  r8_min_test ( );
  r8_mod_test ( );
  r8_modp_test ( );
  r8_mop_test ( );
  r8_nint_test ( );
  r8_normal_01_test ( );
  r8_pi_test ( );
  r8_power_test ( );
  r8_power_fast_test ( );
  r8_rise_test ( );
  r8_round2_test ( );
  r8_roundb_test ( );
  r8_roundx_test ( );
  r8_secd_test ( );
  r8_sign_test ( );
  r8_sign3_test ( );
  r8_sind_test ( );
  r8_swap_test ( );
  r8_swap3_test ( );
  r8_tand_test ( );
  r8_to_r8_discrete_test ( );
  r8_to_i4_test ( );
  r8_uniform_ab_test ( );
  r8_uniform_01_test ( );
  r8_walsh_1d_test ( );
  r8_wrap_test ( );

  r82col_print_part_test ( );

  r82poly2_type_test ( );

  r82row_order_type_test ( );
  r82row_part_quick_a_test ( );
  r82row_print_part_test ( );
  r82row_sort_heap_index_a_test ( );
  r82row_sort_quick_a_test ( );

  r83col_print_part_test ( );

  r83row_print_part_test ( );

  r8block_expand_linear_test ( );
  r8block_print_test ( );

  r8col_find_test ( );
  r8col_insert_test ( );
  r8col_max_test ( );
  r8col_mean_test ( );
  r8col_min_test ( );
  r8col_permute_test ( );
  r8col_sort_heap_a_test ( );
  r8col_sort_heap_index_a_test ( );
  r8col_sort_quick_a_test ( );
  r8col_sorted_tol_unique_test ( );
  r8col_sorted_tol_unique_count_test ( );
  r8col_sorted_tol_undex_test ( );
  r8col_sortr_a_test ( );
  r8col_sum_test ( );
  r8col_swap_test ( );
  r8col_to_r8vec_test ( );
  r8col_tol_undex_test ( );
  r8col_undex_test ( );
  r8col_unique_count_test ( );
  r8col_variance_test ( );

  r8int_to_i4int_test ( );

  r8r8vec_index_insert_unique_test ( );

  r8r8r8vec_index_insert_unique_test ( );

  r8mat_cholesky_inverse_test ( );
  r8mat_cholesky_solve_test ( );
  r8mat_cholesky_solve_upper_test ( );
  r8mat_det_2d_test ( );
  r8mat_det_3d_test ( );
  r8mat_det_4d_test ( );
  r8mat_det_5d_test ( );
  r8mat_expand_linear_test ( );
  r8mat_expand_linear2_test ( );
  r8mat_fs_test ( );
  r8mat_fss_test ( );
  r8mat_givens_post_test ( );
  r8mat_givens_pre_test ( );
  r8mat_hess_test ( );
  r8mat_house_axh_test ( );
  r8mat_house_form_test ( );
  r8mat_house_post_test ( );
  r8mat_house_pre_test ( );
  r8mat_indicator_test ( );
  r8mat_inverse_2d_test ( );
  r8mat_inverse_3d_test ( );
  r8mat_inverse_4d_test ( );
  r8mat_jac_test ( );
  r8mat_kronecker_test ( );
  r8mat_l_inverse_test ( );
  r8mat_l_print_test ( );
  r8mat_l1_inverse_test ( );
  r8mat_lu_test ( );
  r8mat_max_test ( );
  r8mat_max_index_test ( );
  r8mat_min_test ( );
  r8mat_min_index_test ( );
  r8mat_maxcol_minrow_test ( );
  r8mat_maxrow_mincol_test ( );
  r8mat_mincol_maxrow_test ( );
  r8mat_minrow_maxcol_test ( );
  r8mat_mm_test ( );
  r8mat_mv_test ( );
  r8mat_nint_test ( );
  r8mat_nonzeros_test ( );
  r8mat_norm_fro_test ( );
  r8mat_norm_l1_test ( );
  r8mat_nullspace_test ( );
  r8mat_nullspace_size_test ( );
  r8mat_orth_uniform_test ( );
  r8mat_plot_test ( );
  r8mat_power_method_test ( );
  r8mat_print_test ( );
  r8mat_print_some_test ( );
  r8mat_ref_test ( );
  r8mat_rref_test ( );
  r8mat_solve_test ( );
  r8mat_solve_2d_test ( );
  r8mat_solve_3d_test ( );
  r8mat_solve2_test ( );
  r8mat_sub_test ( );
  r8mat_symm_jacobi_test ( );
  r8mat_to_r8plu_test ( );
  r8mat_trace_test ( );
  r8mat_transpose_print_test ( );
  r8mat_transpose_print_test ( );
  r8mat_u_inverse_test ( );
  r8mat_u1_inverse_test ( );
  r8mat_uniform_ab_test ( );

  r8plu_det_test ( );
  r8plu_inverse_test ( );
  r8plu_mul_test ( );
  r8plu_sol_test ( );
  r8plu_to_r8mat_test ( );

  r8poly_degree_test ( );
  r8poly_deriv_test ( );
  r8poly_lagrange_coef_test ( );
  r8poly_lagrange_0_test ( );
  r8poly_lagrange_1_test ( );
  r8poly_lagrange_2_test ( );
  r8poly_lagrange_factor_test ( );
  r8poly_lagrange_val_test ( );
  r8poly_print_test ( );
  r8poly_value_horner_test ( );
  r8poly_values_horner_test ( );

  r8poly2_ex_test ( );
  r8poly2_ex2_test ( );
  r8poly2_val_test ( );
  r8poly2_val2_test ( );
  r8poly2_root_test ( );

  r8poly3_root_test ( );

  r8poly4_root_test ( );

  r8row_max_test ( );
  r8row_mean_test ( );
  r8row_min_test ( );
  r8row_sum_test ( );
  r8row_swap_test ( );
  r8row_to_r8vec_test ( );
  r8row_variance_test ( );

  r8slmat_print_test ( );

  r8vec_amax_test ( );
  r8vec_amin_test ( );
  r8vec_bracket_test ( );
  r8vec_bracket2_test ( );
  r8vec_bracket3_test ( );
  r8vec_bracket5_test ( );
  r8vec_chebyspace_test ( );
  r8vec_concatenate_test ( );
  r8vec_convolution_test ( );
  r8vec_convolution_circ_test ( );
  r8vec_dif_test ( );
  r8vec_direct_product_test ( );
  r8vec_direct_product2_test ( );
  r8vec_even_test ( );
  r8vec_even2_test ( );
  r8vec_even3_test ( );
  r8vec_expand_linear_test ( );
  r8vec_frac_test ( );
  r8vec_histogram_test ( );
  r8vec_house_column_test ( );
  r8vec_index_delete_all_test ( );
  r8vec_index_delete_dupes_test ( );
  r8vec_index_delete_one_test ( );
  r8vec_index_insert_test ( );
  r8vec_index_insert_unique_test ( );
  r8vec_index_order_test ( );
  r8vec_index_search_test ( );
  r8vec_index_sorted_range_test ( );
  r8vec_indexed_heap_d_test ( );
  r8vec_indexed_heap_d_extract_test ( );
  r8vec_indexed_heap_d_insert_test ( );
  r8vec_indexed_heap_d_max_test ( );
  r8vec_indicator0_test ( );
  r8vec_indicator1_test ( );
  r8vec_legendre_test ( );
  r8vec_linspace_test ( );
  r8vec_max_test ( );
  r8vec_max_index_test ( );
  r8vec_mean_test ( );
  r8vec_median_test ( );
  r8vec_midspace_test ( );
  r8vec_min_test ( );
  r8vec_min_index_test ( );
  r8vec_nint_test ( );
  r8vec_norm_l0_test ( );
  r8vec_norm_l1_test ( );
  r8vec_norm_l2_test ( );
  r8vec_norm_li_test ( );
  r8vec_normal_01_test ( );
  r8vec_normalize_l1_test ( );
  r8vec_order_type_test ( );
  r8vec_permute_test ( );
  r8vec_permute_uniform_test ( );
  r8vec_polarize_test ( );
  r8vec_print_test ( );
  r8vec_rotate_test ( );
  r8vec_reverse_test ( );
  r8vec_search_binary_a_test ( );
  r8vec_sort_bubble_a_test ( );
  r8vec_sort_heap_a_test ( );
  r8vec_sort_heap_d_test ( );
  r8vec_sort_heap_index_a_test ( );
  r8vec_sort_heap_index_d_test ( );
  r8vec_sort_heap_mask_a_test ( );
  r8vec_sort_insert_a_test ( );
  r8vec_sort_insert_index_a_test ( );
  r8vec_sort_quick_a_test ( );
  r8vec_sorted_merge_a_test ( );
  r8vec_sorted_nearest_test ( );
  r8vec_sorted_range_test ( );
  r8vec_sorted_split_test ( );
  r8vec_sorted_undex_test ( );
  r8vec_sorted_unique_test ( );
  r8vec_sorted_unique_count_test ( );
  r8vec_sorted_unique_hist_test ( );
  r8vec_split_test ( );
  r8vec_transpose_print_test ( );
  r8vec_undex_test ( );
  r8vec_uniform_01_test ( );
  r8vec_uniform_ab_test ( );
  r8vec_variance_test ( );

  r8vec2_sort_a_test ( );
  r8vec2_sort_d_test ( );
  r8vec2_sort_heap_index_a_test ( );
  r8vec2_sorted_unique_test ( );
  r8vec2_sorted_unique_index_test ( );
  r8vec2_sum_max_index_test ( );

  roots_to_r8poly_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8LIB_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
%%% Projection onto Feasible motion space

% Course: Robotic Manipulation and Mobility
% Advisor: Dr. V. Krovi
%
% Homework Number: MIDTERM
%
% Names: Sourish Chakravarty
% 	Hrishi Lalit Shah

function [dX]= ROBO_midterm_f1(t,X)

global l1 lc1 m1 j1 tau1 l2 lc2 m2 j2 tau2 g

th1=X(1);
x2= X(2);
y2= X(3);
th2=X(4);
th1d= X(5);
th2d= X(6);
%%%%%%%%%%%%%%%%%%%%%%% CREATING IMPORTANT MATRICES IN THE GOVERNING EQUATION
s1=sin(th1);
s2=sin(th2);
c1=cos(th1);
c2=cos(th2);

%%%% Element - 1
M1= [j1+m1*(lc1^2)];
E1= [1,-1];
G1= [m1*lc1*c1*g];
% B11 = [-l1*s1, l1*c1];
% B12 = [-m1*lc1*c1];

%%%% Element - 2
M2= [m2, 0, 0;
     0, m2, 0;
     0, 0, j2];
E2= [ 0; 0; 1];
G2= [0; m2*g; 0];
% B21= [-1, 0;
%      0, -1;
%      -lc2*s2, lc2*c2];
% B22= [0, -m2, 0]';

%%%% Constraints
C= [x2-l1*c1-lc2*c2;
    y2-l1*s1-lc2*s2];
A= [l1*s1, 1, 0, lc2*s2;
    -l1*c1, 0, 1, -lc2*c2];% Jacobian of constraint matrix
S = [1, 0;
    -l1*s1, -lc2*s2;
    l1*c1, lc2*c2;
    0, 1]; % Null space of A or Feasible Motion
Sd = [0, 0;
    -l1*c1*th1d, -lc2*c2*th2d;
    -l1*s1*th1d, -lc2*s2*th2d;
    0, 0]; % Derivative of S matrix

%%%%%%%%%%%%%%%%%%%%%%% PROJECTION ONTO FEASIBLE MOTION SPACE
M=[M1, zeros(1,3);
    zeros(3,1), M2];
E=[E1;
    zeros(3,1),E2];
G= [G1;G2];
T =[tau1; tau2];
V=zeros(4,1);

t1= inv(S'*M*S);
t2= S'*(M*Sd*[th1d,th2d]'+V+G-E*T);

dX(1:4,1)=S*[th1d,th2d]';
dX(5:6,1)= -t1*t2;
% Cerr1=([Cerr1;abs(C')]);
% Tstore1=[Tstore1,t];
return


text:
function [area smoothedArea] = roiSurfaceArea(vw, roi, msh)
% Get the surface area of an ROI on a mesh
%
% [area smoothedArea] = roiSurfaceArea(vw, roi, msh)
%

% Argument checks
if notDefined('vw'), vw = getSelectedGray; end
switch lower(viewGet(vw, 'viewType'))
    case {'gray' 'volume'}
        % ok
    otherwise
        error('Must be in Gray/Volume view to get ROI surface area');
end
if notDefined('msh'), msh = viewGet(vw, 'curmesh'); end
if notDefined('roi'), roi = viewGet(vw, 'selectedROI'); end

% Set current mesh in view struct
if isnumeric(msh),
    vw = viewSet(vw, 'curmeshn', msh);
    msh = viewGet(vw, 'curmesh');
end

verts = viewGet(vw, 'roivertinds', roi);

if isempty(verts),
    area = 0;
    smoothedArea = 0;
    return;
end

[areaList, smoothAreaList] = mrmComputeMeshArea(msh, verts);
area            = sum(areaList);
smoothedArea    = sum(smoothAreaList);

fprintf('ROI surface area: %0.1f mm^2 (%0.1f mm^2 on smoothed mesh)\n', area, smoothedArea);

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_suite = test_drawLine
% Tests the function 'drawLine'
%   output = test_drawLine(input)
%
%   Example
%   testClipLine
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-04-22,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.
% Licensed under the terms of the LGPL, see the file "license.txt"

test_suite = functiontests(localfunctions);

function testHoriz(testCase) %#ok<*DEFNU>
% test clipping of horizontal lines

box = [0 100 0 100];

hf = figure; clf;
axis(box);

% inside, to the right
line = [30 40 10 0];
edge = [0 40 100 40];

hl = drawLine(line);
testCase.assertEqual(edge([1 3]), get(hl, 'xdata'), 'AbsTol', .01);
testCase.assertEqual(edge([2 4]), get(hl, 'ydata'), 'AbsTol', .01);

% inside, to the left
line = [30 40 -10 0];
edge = [100 40 0 40];
hl = drawLine(line);
testCase.assertEqual(edge([1 3]), get(hl, 'xdata'), 'AbsTol', .01);
testCase.assertEqual(edge([2 4]), get(hl, 'ydata'), 'AbsTol', .01);

% outside
line = [30 140 10 0];
hl = drawLine(line);
testCase.assertEqual(-1, hl);

close(hf);

function testVert(testCase)
% test clipping of vertical lines

box = [0 100 0 100];

hf = figure; clf;
axis(box);

% inside, upward
line = [30 40 0 10];
edge = [30 0 30 100];
hl = drawLine(line);
testCase.assertEqual(edge([1 3]), get(hl, 'xdata'), 'AbsTol', .01);
testCase.assertEqual(edge([2 4]), get(hl, 'ydata'), 'AbsTol', .01);

% inside, downward
line = [30 40 0 -10];
edge = [30 100 30 0];
hl = drawLine(line);
testCase.assertEqual(edge([1 3]), get(hl, 'xdata'), 'AbsTol', .01);
testCase.assertEqual(edge([2 4]), get(hl, 'ydata'), 'AbsTol', .01);

% outside
line = [140 30 0 10];
hl = drawLine(line);
testCase.assertEqual(-1, hl);

close(hf);

function testDiagUp(testCase)
% test clipping of upward diagonal lines

box = [0 100 0 100];

hf = figure; clf;
axis(box); hold on;

% inside, top right corner
line = [80 30 10 10];
edge = [50 0 100 50];
hl = drawLine(line);
testCase.assertEqual(edge([1 3]), get(hl, 'xdata'), 'AbsTol', .01);
testCase.assertEqual(edge([2 4]), get(hl, 'ydata'), 'AbsTol', .01);

% inside, down right corner
line = [20 70 10 10];
edge = [0 50 50 100];
hl = drawLine(line);
testCase.assertEqual(edge([1 3]), get(hl, 'xdata'), 'AbsTol', .01);
testCase.assertEqual(edge([2 4]), get(hl, 'ydata'), 'AbsTol', .01);

% outside
line = [140 -30 10 10];
hl = drawLine(line);
testCase.assertEqual(-1, hl);

line = [-40 130 10 10];
hl = drawLine(line);
testCase.assertEqual(-1, hl);

close(hf);

function testMultiLines(testCase)

box = [0 100 0 100];

hf = figure; clf;
axis(box);

% inside, top right corner
line = [...
    80 30 10 10; ...
    20 70 10 10; ...
    140 -30 10 10; ...
    -40 130 10 10];
edge = [...
    50 0 100 50; ...
    0 50 50 100];
hl = drawLine(line);
testCase.assertEqual(4, length(hl));
testCase.assertEqual(edge(1, [1 3]), get(hl(1), 'xdata'), 'AbsTol', .01);
testCase.assertEqual(edge(1, [2 4]), get(hl(1), 'ydata'), 'AbsTol', .01);
testCase.assertEqual(edge(2, [1 3]), get(hl(2), 'xdata'), 'AbsTol', .01);
testCase.assertEqual(edge(2, [2 4]), get(hl(2), 'ydata'), 'AbsTol', .01);
testCase.assertEqual(-1, hl(3));
testCase.assertEqual(-1, hl(4));

close(hf);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [llh] = tapas_dlinear_hier_llh(data, theta, ptheta)
%% Likelihood of the nodes of a model with diagonal precision matrix.
% Note that this is the likelihood of several parameters theta 1 to n with mean
% mu and precision matrix diag(pe).
%

%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

y = data.y;
u = data.u;

theta = theta.y;

[np, nc] = size(y);

llh = zeros(np, nc);

ln2pi = log(2 * pi);

for j = 1:nc
    % Mean
    mu = theta{j}.mu;
    % Precision
    pe = theta{j}.pe;
    lpe = log(pe);

   for i = 1:np
        r = y{i, j} - mu;
        llh(i, j) = sum(- 0.5 * ln2pi + 0.5 * lpe - 0.5 * pe .* r .* r);
    end
end

end

author: translationalneuromodeling

repo: tapas


text:
classdef GenericSumOfSinesSteeringModel < AbstractSteeringModel
    %GenericSumOfSinesSteeringModel Summary of this class goes here
    %   Detailed explanation goes here

    properties
        gammaAngleModel(1,1) SumOfSinesModel = SumOfSinesModel(0);
        betaAngleModel(1,1) SumOfSinesModel = SumOfSinesModel(0);
        alphaAngleModel(1,1) SumOfSinesModel = SumOfSinesModel(0);

        gammaContinuity(1,1) logical = false;
        betaContinuity(1,1) logical = false;
        alphaContinuity(1,1) logical = false;

        refFrame AbstractReferenceFrame
        controlFrame AbstractControlFrame
    end

    methods
        function dcm = getBody2InertialDcmAtTime(obj, ut, rVect, vVect, bodyInfo)
            gammaAng = obj.gammaAngleModel.getValueAtTime(ut);
            betaAng = obj.betaAngleModel.getValueAtTime(ut);
            alphaAng = obj.alphaAngleModel.getValueAtTime(ut);

%             elemSet = CartesianElementSet(ut, rVect(:), vVect(:), bodyInfo.getBodyCenteredInertialFrame());
%             if(not(isempty(obj.refFrame.getOriginBody())))
%                 elemSet = elemSet.convertToFrame(obj.refFrame);
%             end

            dcm = real(obj.controlFrame.computeDcmToInertialFrame(ut, rVect, vVect, bodyInfo, gammaAng, betaAng, alphaAng, obj.refFrame));
        end

        function [angleModel, continuity] = getAngleNModel(obj, n)
            angleModel = SumOfSinesModel.empty(1,0);

            switch n
                case 1
                    angleModel = obj.gammaAngleModel;
                    continuity = obj.gammaContinuity;
                case 2
                    angleModel = obj.betaAngleModel;
                    continuity = obj.betaContinuity;
                case 3
                    angleModel = obj.alphaAngleModel;
                    continuity = obj.alphaContinuity;
            end
        end

        function t0 = getT0(obj)
            t0 = obj.alphaAngleModel.getT0();
        end

        function setT0(obj, newT0)
            obj.gammaAngleModel.setT0(newT0);
            obj.betaAngleModel.setT0(newT0);
            obj.alphaAngleModel.setT0(newT0);
        end

        function setTimeOffsets(obj, timeOffset)
            obj.gammaAngleModel.setTimeOffset(timeOffset);
            obj.betaAngleModel.setTimeOffset(timeOffset);
            obj.alphaAngleModel.setTimeOffset(timeOffset);
        end

        function [angle1Cont, angle2Cont, angle3Cont] = getContinuityTerms(obj)
            angle1Cont = obj.gammaContinuity;
            angle2Cont = obj.betaContinuity;
            angle3Cont = obj.alphaContinuity;
        end

        function setContinuityTerms(obj, angle1Cont, angle2Cont, angle3Cont)
            obj.gammaContinuity = angle1Cont;
            obj.betaContinuity = angle2Cont;
            obj.alphaContinuity = angle3Cont;
        end

        function setConstsFromDcmAndContinuitySettings(obj, dcm, ut, rVect, vVect, bodyInfo)
            if(obj.gammaContinuity || obj.betaContinuity || obj.alphaContinuity)
                elemSet = CartesianElementSet(ut, rVect(:), vVect(:), bodyInfo.getBodyCenteredInertialFrame());

%                 if(not(isempty(obj.refFrame.getOriginBody())))
%                     elemSet = elemSet.convertToFrame(obj.refFrame);
%                 end

                [gammaAngle, betaAngle, alphaAngle] = obj.controlFrame.getAnglesFromInertialBodyAxes(LaunchVehicleAttitudeState(dcm), elemSet.time, elemSet.rVect(:), elemSet.vVect(:), elemSet.frame.getOriginBody(), obj.refFrame);

                if(obj.gammaContinuity)
                    obj.gammaAngleModel.setConstValueForContinuity(gammaAngle);
                end

                if(obj.betaContinuity)
                    obj.betaAngleModel.setConstValueForContinuity(betaAngle);
                end

                if(obj.alphaContinuity)
                    obj.alphaAngleModel.setConstValueForContinuity(alphaAngle);
                end
            end
        end

        function setInitialAttitudeFromState(obj, stateLogEntry, tOffsetDelta)
            t0 = stateLogEntry.time;
            obj.setT0(t0);

            obj.gammaAngleModel.setTimeOffset(obj.gammaAngleModel.getTimeOffset() + tOffsetDelta);
            obj.betaAngleModel.setTimeOffset(obj.betaAngleModel.getTimeOffset() + tOffsetDelta);
            obj.alphaAngleModel.setTimeOffset(obj.alphaAngleModel.getTimeOffset() + tOffsetDelta);
        end

        function [angle1Name, angle2Name, angle3Name] = getAngleNames(obj)
            angleNames = obj.controlFrame.getControlFrameEnum().angleNames;
            angle1Name = angleNames{1};
            angle2Name = angleNames{2};
            angle3Name = angleNames{3};
        end

        function newSteeringModel = deepCopy(obj)
            newSteeringModel = GenericSumOfSinesSteeringModel(obj.gammaAngleModel.deepCopy(), obj.betaAngleModel.deepCopy(), obj.alphaAngleModel.deepCopy());

            newSteeringModel.gammaContinuity = obj.gammaContinuity;
            newSteeringModel.betaContinuity = obj.betaContinuity;
            newSteeringModel.alphaContinuity = obj.alphaContinuity;
            newSteeringModel.refFrame = obj.refFrame;
            newSteeringModel.controlFrame = obj.controlFrame;
        end

        function optVar = getNewOptVar(obj)
            optVar = SetGenericSumOfSinesSteeringModelActionOptimVar(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function tf = usesRefFrame(~)
            tf = true;
        end

        function refFrame = getRefFrame(obj)
            refFrame = obj.refFrame;
        end

        function setRefFrame(obj, refFrame)
            obj.refFrame = refFrame;
        end

        function tf = usesControlFrame(~)
            tf = true;
        end

        function cFrame = getControlFrame(obj)
            cFrame = obj.controlFrame;
        end

        function setControlFrame(obj, cFrame)
            obj.controlFrame = cFrame;
        end

        function tf = requiresReInitAfterSoIChange(~)
            tf = false;
        end

        function enum = getSteeringModelTypeEnum(~)
            enum = SteerModelTypeEnum.SumOfSinesAngles;
        end
    end

    methods(Access=private)
        function obj = GenericSumOfSinesSteeringModel(gammaAngleModel, betaAngleModel, alphaAngleModel)
            obj.gammaAngleModel = gammaAngleModel;
            obj.betaAngleModel = betaAngleModel;
            obj.alphaAngleModel = alphaAngleModel;

            obj.controlFrame = InertialControlFrame();
        end
    end

    methods(Static)
        function model = getDefaultSteeringModel()
            gammaAngleModel = SumOfSinesModel(0);
            betaAngleModel = SumOfSinesModel(0);
            alphaAngleModel = SumOfSinesModel(0);

            model = GenericSumOfSinesSteeringModel(gammaAngleModel, betaAngleModel, alphaAngleModel);
        end

        function typeStr = getTypeNameStr()
            typeStr = SteeringModelEnum.GenericSumOfSines.nameStr;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [theta] = normalEqn(X, y)
%NORMALEQN Computes the closed-form solution to linear regression
%   NORMALEQN(X,y) computes the closed-form solution to linear
%   regression using the normal equations.

theta = pinv(X' * X) * X' * y;

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the code to compute the closed form solution
%               to linear regression and put the result in theta.
%

% ---------------------- Sample Solution ----------------------

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

% ============================================================

end

author: zsiciarz

repo: ml-coursera


text:
function [ftH2O] = Pa2ftH2O(Pa)
% Convert pressure from pascals to feet of 4-degree-C water.
% Chad Greene 2012
ftH2O = Pa*0.000334553;

author: Sable


text:
% Test file for chebtech/diff.m

function pass = test_diff(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebtech.techPref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
x = 2 * rand(100, 1) - 1;

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    %%
    % Spot-check derivatives for a couple of functions.

    f = testclass.make(@(x) exp(x) - x, [], pref);
    df = diff(f);
    df_exact = @(x) exp(x) - 1;
    err = norm(df_exact(x) - feval(df, x), inf);
    tol = 100*vscale(df)*eps;
    pass(n, 1) = err < tol;

    f = testclass.make(@(x) atan(x), [], pref);
    df = diff(f);
    df_exact = @(x) 1./(1 + x.^2);
    err = norm(df_exact(x) - feval(df, x), inf);
    tol = 500*vscale(df)*eps;
    pass(n, 2) = err < 10*tol;

    f = testclass.make(@(x) sin(x), [], pref);
    df = diff(f);
    df_exact = @(x) cos(x);
    err = norm(df_exact(x) - feval(df, x), inf);
    tol = 100*vscale(df)*eps;
    pass(n, 3) = err < tol;

    z = exp(2*pi*1i/3);
    f = testclass.make(@(t) airy(z*t), [], pref);
    df = diff(f);
    df_exact = @(t) z*airy(1, z*t);
    err = norm(df_exact(x) - feval(df, x), inf);
    tol = 1e3*vscale(df)*eps;
    pass(n, 4) = err < tol;

    %%
    % Verify that calling diff() gives the same answer as direct construction.

    f = testclass.make(@(x) 0.5*x - 0.0625*sin(8*x), [], pref);
    df = testclass.make(@(x) sin(4*x).^2, [], pref);
    err = diff(f) - df;
    pass(n, 5) = (norm(err.coeffs, inf) < 100*vscale(df)*eps);

    %%
    % Verify basic differentiation rules.

    f = testclass.make(@(x) x.*sin(x.^2) - 1, [], pref);
    df = diff(f);
    g = testclass.make(@(x) exp(-x.^2), [], pref);
    dg = diff(g);
    tol_f = 10*vscale(df)*eps;
    tol_g = 10*vscale(dg)*eps;

    errfn = diff(f + g) - (df + dg);
    err = feval(errfn, x);
    pass(n, 6) = (norm(err, inf) < max(tol_f, tol_g));

    errfn = diff(f.*g) - (f.*dg + g.*df);
    err = feval(errfn, x);
    pass(n, 7) = (norm(err, inf) < 10*length(f)*max(tol_f, tol_g));

    const = testclass.make(@(x) ones(size(x)), [], pref);
    dconst = diff(const);
    err = feval(dconst, x);
    pass(n, 8) = (norm(err, inf) == 0);

    %%
    % Check higher-order derivatives.  (NB:  We relax the tolerance by n + 1
    % factors of 10, where n is the number of derivatives taken.)

    f = testclass.make(@(x) x.*atan(x) - x - 0.5*log(1 + x.^2), [], pref);
    df2 = diff(f, 2);
    df2_exact = @(x) 1./(1 + x.^2);
    err = df2_exact(x) - feval(df2, x);
    pass(n, 9) = (norm(err, inf) < 1e6*vscale(df2)*eps);

    f = testclass.make(@(x) sin(x), [], pref);
    df4 = diff(f, 4);
    df4_exact = @(x) sin(x);
    err = norm(df4_exact(x) - feval(df4, x), inf);
    tol = 1e7*vscale(df4)*eps;
    pass(n, 10) = err < tol;

    f = testclass.make(@(x) x.^5 + 3*x.^3 - 2*x.^2 + 4, [], pref);
    df6 = diff(f, 6);
    df6_exact = @(x) zeros(size(x));
    err = df6_exact(x) - feval(df6, x);
    tol = 900*vscale(df4)*eps;
    pass(n, 11) = (norm(err, inf) == 0);

    %%
    % Check operation for array-valued chebtech objects.

    f = testclass.make(@(x) [sin(x) x.^2 exp(1i*x)], [], pref);
    df = diff(f);
    df_exact = @(x) [cos(x) 2*x 1i*exp(1i*x)];
    err = feval(df, x) - df_exact(x);
    pass(n, 12) = (norm(err(:), inf) < 1e2*max(vscale(df)*eps));

    % DIM option.
    dim2df = diff(f, 1, 2);
    g = @(x) [(x.^2 - sin(x)) (exp(1i*x) - x.^2)];
    err = feval(dim2df, x) - g(x);
    pass(n, 13) = isequal(size(vscale(dim2df)), [1 2]) && ...
        (norm(err(:), inf) < 10*max(vscale(dim2df)*eps));

    dim2df2 = diff(f, 2, 2);
    g = @(x) exp(1i*x) - 2*x.^2 + sin(x);
    err = feval(dim2df2, x) - g(x);
    pass(n, 14) = isequal(size(vscale(dim2df2)), [1 1]) && ...
        (norm(err(:), inf) < 10*max(vscale(dim2df2)*eps));

    % DIM option should return an empty chebtech for non-array-valued input.
    f = testclass.make(@(x) x.^3);
    dim2df = diff(f, 1, 2);
    pass(n, 15) = (isempty(dim2df.coeffs));

    % Check for #1641.
    f = testclass.make(@(x) [1 + x + x.^2, 1 - x + 2*x.^2]);
    df = diff(f);
    err = norm(df.coeffs - [1 -1 ; 2 4], 'fro');
    tol = 10*eps;
    pass(n, 16) = err < tol;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function g=TaylorLinearTapering(nBar,sidelobedB,xPoints,a)
%%TAYLORLINEARTAPERING The Taylor tapering is a set of amplitude weights
%          for a continuous linear (narrowband) aperture antenna that will
%          reduce the height of nBar of the close-in sidelobes at a cost of
%          widening the beam. Such a tapering can be discretized and
%          applied to the elements in a linear phased array (An array of
%          antenna elements can be viewed as a discrete approximation to a
%          continuous aperture). This function will provide the tapering
%          values at a set of discrete points given by xPoints (the origin
%          is taken to be the center of the aperture). The radius of the
%          aperture can either be provided or is taken as value of the
%          farthest point provided.
%
%INPUTS: nBar The number of terms to use in the expansion. High values can
%             produce undesirable results.
%  sidelobedB The number of decibels of the ratio of the close-in sidelobe
%             voltages to the main lobe voltage. This must be a negative
%             number. A typical value is -30.
%    xyPoints An NX1 or 1XN set of N points at which the tapering values
%             should be evaluated. The center of the aperture is taken to
%             be the origin.
%           a The radius of the aperture. Tapering weights for points in
%             xPoints outside of the aperture are taken to be zero. If
%             this parameter is omitted or an empty matrix is passed, then
%             the radius is taken to be the distance of the farthest point
%             from the origin in xPoints.
%
%OUTPUTS: g The NX1 set of discretized Taylor tapering values evaluated at
%           the points given in xPoints. All Taylor tapering values are
%           positive and real. The coefficients are not normalized.
%
%This function implements the algorithm in [1].
%
%EXAMPLE 1:
%Here, we evaluate the tapering values for 30dB down on a large array of
%points to see what the tapering looks like.
% nBar=4;
% sidelobedB=-30;
% numPoints=100;
% %Generate points symmetric about the origin. This makes the aperture size
% %1.
% k=0:(numPoints-1);
% xPoints=(k-(1/2)*numPoints+1/2)/numPoints;
% g=TaylorLinearTapering(nBar,sidelobedB,xPoints);
%
% figure(1)
% clf
% plot(xPoints,g)
% h1=xlabel('x');
% h2=ylabel('y');
% title('Taylor Tapering Weight')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%EXAMPLE 2:
%Here, we consider the array response when using tapering values for 30dB
%sidelobes on a linear array with lambda/2 spacing between elements.
% nBar=4;
% sidelobedB=-30;
% numPoints=30;
% %Generate points symmetric about the origin.
% k=0:(numPoints-1);
% xPoints=(k-(1/2)*numPoints+1/2)/2;
% g=TaylorLinearTapering(nBar,sidelobedB,xPoints);
%
% %Tapering matrix
% T=diag(g);
%
% %Now, display the response with the tapering
% [Rsp,U]=standardUVBeamPattern(T,xPoints,'NormPowGain');
%
% figure(2)
% clf
% plot(U,10*log10(Rsp),'-b','LineWidth',2);
% axis([-1, 1, -40 1])
% h1=xlabel('u');
% h2=ylabel('Amplitude Response (decibels)');
% title('Taylor Weighted Array Response')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] T. T. Taylor, "Design of line-source antennas for narrow beamwidth and
%    low side lobes," IRE Transactions on Antennas and Propagation, vol. 3,
%    no. 1, pp. 16-28, Jan. 1955.
%
%August 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%After Equation 35 in [1], it is noted that cosh(pi*A) is the side-lobe
%ratio. Thus, we get the following expression for A in terms of the
%sidelobe level in decibels:
eta=10^(sidelobedB/(-20));
A=acosh(eta)/pi;

%The square of Sigma in Equation 41 in [1].
sigma2=nBar^2/(A^2+(nBar-1/2)^2);

%The suggestion for C after Equation 48 in [1].
C=cosh(pi*A);

%Fm holds values of F in Equation 48 for z being an integer>=1.
Fm=zeros(nBar-1,1);
for m=1:(nBar-1)
%Equation 48 contains singularities when z=m is an integer. To evaluate it
%when z is an integer, we first simplify Q in Equation 47 to eliminate the
%singularity.

    %Equation 47 in [1], taking the limit of z to an integer to eliminate
    %the singularity.
    np=[1:(m-1),(m+1):(nBar-1)];
    Qm=(-1)^(m+1)/(2*prod(1-(m./np).^2));

    n=1:(nBar-1);
    Fm(m)=C*Qm*prod(1-m^2./(sigma2*(A^2+(n-1/2).^2)));
end
F0=C;%The case where z=0.

numPoints=length(xPoints);

if(nargin<4||isempty(a))
    %The maximum distance from the origin to a point is taken to be the
    %radius of the aperture.
    a=sqrt(max(sum(xPoints.^2,1)));
end

%The loop below implements Equation 63 in [1].
g=zeros(numPoints,1);
g(:)=F0/(2*pi);
for curPoint=1:numPoints
    rho=norm(xPoints(curPoint));

    if(rho<=a)
        %The normalized radius at this point.
        p=pi*rho/a;

        m=(1:(nBar-1))';
        g(curPoint)=g(curPoint)+(1/pi)*sum(Fm.*cos(m*p));
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [ n_data_new, x, y, fxy ] = beta_values ( n_data )

%*****************************************************************************80
%
%% BETA_VALUES returns some values of the Beta function.
%
%  Formula:
%
%    BETA(X,Y) = ( GAMMA(X) * GAMMA(Y) ) / GAMMA(X+Y)
%
%  Restrictions:
%
%    Both X and Y must be greater than 0.
%
%  Properties:
%
%    BETA(X,Y) = BETA(Y,X).
%    BETA(X,Y) = Integral ( 0 <= T <= 1 ) T**(X-1) (1-T)**(Y-1) dT.
%    BETA(X,Y) = GAMMA(X) * GAMMA(Y) / GAMMA(X+Y)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N_DATA, indicates the index of the previous test data
%    returned, or is 0 if this is the first call.  For repeated calls,
%    set the input value of N_DATA to the output value of N_DATA_NEW
%    from the previous call.
%
%    Output, integer N_DATA_NEW, the index of the test data.
%
%    Output, real X, Y, the arguments of the function.
%
%    Output, real FXY, the value of the function.
%
  n_max = 17;
  b_vec = [ ...
    5.000000E+00, 2.500000E+00, 1.666667E+00, 1.250000E+00, ...
    5.000000E+00, 2.500000E+00, 1.000000E+00, 1.666667E-01, ...
    0.333333E-01, 7.142857E-03, 1.587302E-03, 0.238095E-01, ...
    5.952381E-03, 1.984127E-03, 7.936508E-04, 3.607504E-04, ...
    8.325008E-05 ];
  x_vec = [ ...
    0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, ...
    1.0E+00, 1.0E+00, 1.0E+00, 2.0E+00, ...
    3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00, ...
    6.0E+00, 6.0E+00, 6.0E+00, 6.0E+00, ...
    7.0E+00 ];
  y_vec = [ ...
    1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, ...
    0.2E+00, 0.4E+00, 1.0E+00, 2.0E+00, ...
    3.0E+00, 4.0E+00, 5.0E+00, 2.0E+00, ...
    3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00, ...
    7.0E+00 ];

  n_data_new = n_data;

  if ( n_data_new < 0 )
    n_data_new = 0;
  end

  n_data_new = n_data_new + 1;

  if ( n_max < n_data_new )
    n_data_new = 0;
    x = 0.0E+00;
    y = 0.0E+00;
    fxy = 0.0E+00;
  else
    x = x_vec(n_data_new);
    y = y_vec(n_data_new);
    fxy = b_vec(n_data_new);
  end

  return
end


text:
function r = reverse(x)
% reverse -- Reverse order of elements in 1-d signal
%  Usage
%    r = reverse(x)
%  Inputs
%    x     1-d signal
%  Outputs
%    r     1-d time-reversed signal
%
%  See Also
%    flipud, fliplr
%
   r = x(length(x):-1:1);

%
% Copyright (c) 1993. David L. Donoho
%

%
%  Part of Wavelab Version 850
%  Built Tue Jan  3 13:20:40 EST 2006
%  This is Copyrighted Material
%  For Copying permissions see COPYING.m
%  Comments? e-mail wavelab@stat.stanford.edu

author: thomaskuestner


text:
classdef ScalarNodalVariablePrinter < ResultsPrinter

    properties (Access = protected)
        simulationStr
        hasGaussData = false;
    end

    properties (Access = private)
        fieldName
    end

    methods (Access = public)

        function obj = ScalarNodalVariablePrinter(d)
            obj.init(d);
            obj.fieldName = d.fieldName;
            obj.simulationStr = d.fieldName;
        end

        function printResults(obj,iter,fileID)
            f = obj.fields;
            dS = obj.createScalarDataBase(iter,fileID,f,obj.fieldName,'OnNodes');
            ScalarNodalPrinter(dS);
        end

    end

    methods (Access = protected)

        function storeFieldsToPrint(obj,d)
            obj.fields = d.fields;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function i4vec_part_test ( )

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

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_PART_TEST\n' );
  fprintf ( 1, '  I4VEC_PART partitions an integer.\n' );

  nval = 17;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NVAL = %d\n', nval );

  a = i4vec_part ( n, nval );

  i4vec_print ( n, a, '  Partitioned:' )

  nval = -49;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NVAL = %d\n', nval );

  a = i4vec_part ( n, nval );

  i4vec_print ( n, a, '  Partitioned:' );

  return
end

author: johannesgerer


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcPolygonSeg(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

handles.V.Color_map='Gray';
set(handles.Color_map_popmenu,'Value',1);
handles=MU_update_image(handles.Matrix_display_axes,{handles.TMatrix,handles.Mask},handles,0);
MU_enable('off',[],handles);
Seg_h=impoly;
wait(Seg_h);
if ~isvalid(Seg_h) % detect when main display is deleted
    return;
end
handles.V.Segs{end+1,1}='impoly';
handles.V.Segs{end,2}=handles.V.Slice;
handles.V.Segs{end,3}=1;
handles.V.Segs{end,4}=getPosition(Seg_h);
MU_enable('on',{'Color_map_popmenu'},handles);
BW=createMask(Seg_h);
Temp=handles.Mask;
Temp2=Temp(:,:,handles.V.Slice);
Temp2(BW~=0)=1;
Temp(:,:,handles.V.Slice)=Temp2;
handles.Mask=Temp;
handles=MU_update_image(handles.Matrix_display_axes,{handles.TMatrix,handles.Mask},handles,0);
handles.V.ROI=struct(...
                     'ROI_flag', 8,...
                     'ROI_mov',[],...  % ROI movement track
                     'ROI_Stat_h', [],...
                     'ROI_h', Seg_h ...
                     );
guidata(handles.MU_matrix_display, handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [rec,prec,ap] = VOCevaldet(VOCopts,id,cls,draw)
%% VOCevaldet
% inputs: VOCopts = imdb.details.VOCopts

% load test set
[gtids,t]=textread(sprintf(VOCopts.imgsetpath,VOCopts.testset),'%s %d');

% load ground truth objects
tic;
npos=0;
gt(length(gtids))=struct('BB',[],'diff',[],'det',[]);
for i=1:length(gtids)
    % display progress
    if toc>1
        fprintf('%s: pr: load: %d/%d\n',cls,i,length(gtids));
        drawnow;
        tic;
    end

    % read annotation
    rec=PASreadrecord(sprintf(VOCopts.annopath,gtids{i}));

    % extract objects of class
    clsinds=strmatch(cls,lower( {rec.objects(:).class} ),'exact');
    gt(i).BB=cat(1,rec.objects(clsinds).bbox)';
    gt(i).diff=[rec.objects(clsinds).difficult];
    gt(i).det=false(length(clsinds),1);
    npos=npos+sum(~gt(i).diff);
end

% load results
[ids,confidence,b1,b2,b3,b4]=textread(sprintf(VOCopts.detrespath,id,cls),'%s %f %f %f %f %f');
BB=[b1 b2 b3 b4]';

% sort detections by decreasing confidence
[sc,si]=sort(-confidence);
ids=ids(si);
BB=BB(:,si);

% assign detections to ground truth objects
nd=length(confidence);
tp=zeros(nd,1);
fp=zeros(nd,1);
tic;
for d=1:nd
    % display progress
    if toc>1
        fprintf('%s: pr: compute: %d/%d\n',cls,d,nd);
        drawnow;
        tic;
    end

    % find ground truth image
    i=strmatch(ids{d},gtids,'exact');
    if isempty(i)
        error('unrecognized image "%s"',ids{d});
    elseif length(i)>1
        error('multiple image "%s"',ids{d});
    end

    % assign detection to ground truth object if any
    bb=BB(:,d);
    ovmax=-inf;
    for j=1:size(gt(i).BB,2)
        bbgt=gt(i).BB(:,j);
        bi=[max(bb(1),bbgt(1)) ; max(bb(2),bbgt(2)) ; min(bb(3),bbgt(3)) ; min(bb(4),bbgt(4))];
        iw=bi(3)-bi(1)+1;
        ih=bi(4)-bi(2)+1;
        if iw>0 & ih>0
            % compute overlap as area of intersection / area of union
            ua=(bb(3)-bb(1)+1)*(bb(4)-bb(2)+1)+...
               (bbgt(3)-bbgt(1)+1)*(bbgt(4)-bbgt(2)+1)-...
               iw*ih;
            ov=iw*ih/ua;
            if ov>ovmax
                ovmax=ov;
                jmax=j;
            end
        end
    end
    % assign detection as true positive/don't care/false positive
    if ovmax>=VOCopts.minoverlap
        if ~gt(i).diff(jmax)
            if ~gt(i).det(jmax)
                tp(d)=1;            % true positive
		gt(i).det(jmax)=true;
            else
                fp(d)=1;            % false positive (multiple detection)
            end
        end
    else
        fp(d)=1;                    % false positive
    end
end

% compute precision/recall
fp=cumsum(fp);
tp=cumsum(tp);
rec=tp/npos;
prec=tp./(fp+tp);

% compute average precision

ap=0;
for t=0:0.1:1
    p=max(prec(rec>=t));
    if isempty(p)
        p=0;
    end
    ap=ap+p/11;
end

if draw
    % plot precision/recall
    plot(rec,prec,'-');
    grid;
    xlabel 'recall'
    ylabel 'precision'
    title(sprintf('class: %s, subset: %s, AP = %.3f',cls,VOCopts.testset,ap));
end

author: CrazyStoneonRoad

repo: TGRS-HRRSD-Dataset


text:
function [configure, obj] = Estimate_Weight(configure, Seqs)

% initialization
configure.weight = rand(length(configure.id),1);
tau = configure.tau;
obj = zeros(configure.epoch * length(Seqs), 1);

tic
for n = 1:configure.epoch
    ind = randperm(length(Seqs));
    lr = configure.lr * (0.9)^(n-1);
    for m = 1:length(Seqs)
        X = [Seqs(ind(m)).Time; Seqs(ind(m)).Mark];
        [Prob, Delta] = ChosenProbability(X, configure);

        grad = 0;
        for t1 = 1:size(Prob,1)-tau
            for t2 = 1:size(Prob,2)-tau
                if t1~=t2
                    obj((n-1)*length(Seqs)+m) = obj((n-1)*length(Seqs)+m) +...
                       Prob(t1, t2)*Prob(t1+tau, t2+tau);
                    grad = grad + ...
                        2*Prob(t1, t2)*Prob(t1+tau, t2+tau)*...
                        (Delta(:,t1,t2)+Delta(:,t1+tau,t2+tau));
                end
            end
        end
        configure.weight = configure.weight - lr * grad;

        fprintf('epoch=%d, #seq=%d/%d, obj=%f, ||grad||=%.4f, time=%.2fsec\n',...
            n, m, length(Seqs), obj((n-1)*length(Seqs)+m), norm(grad), toc);
    end

end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function y = logsumexp( varargin )

%LOGSUMEXP    log(sum(exp(x))).
%   LOGSUMEXP(X) = LOG_SUM_EXP(X) = LOG(SUM(EXP(X)). We have replaced this
%   function with LOG_SUM_EXP to better match our function naming
%   conventions. Please start using it instead.

warning( 'CVX:Renamed', [ ...
    'The function "logsumexp" has been renamed "log_sum_exp". Please start\n', ...
    'using the new name. The old name will be removed in a future release.' ], 1 );

y = log_sum_exp( varargin{:} );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr


text:
% Demo for Edge Boxes (please see readme.txt first).

%% load pre-trained edge detection model and set opts (see edgesDemo.m)
model=load('models/forest/modelBsds'); model=model.model;
model.opts.multiscale=0; model.opts.sharpen=2; model.opts.nThreads=4;

%% set up opts for edgeBoxes (see edgeBoxes.m)
opts = edgeBoxes;
opts.alpha = .65;     % step size of sliding window search
opts.beta  = .75;     % nms threshold for object proposals
opts.minScore = .01;  % min score of boxes to detect
opts.maxBoxes = 1e4;  % max number of boxes to detect

%% detect Edge Box bounding box proposals (see edgeBoxes.m)
I = imread('peppers.png');
tic, bbs=edgeBoxes(I,model,opts); toc

%% show evaluation results (using pre-defined or interactive boxes)
gt=[122 248 92 65; 193 82 71 53; 410 237 101 81; 204 160 114 95; ...
  9 185 86 90; 389 93 120 117; 253 103 107 57; 81 140 91 63];
if(0), gt='Please select an object box.'; disp(gt); figure(1); imshow(I);
  title(gt); [~,gt]=imRectRot('rotate',0); gt=gt.getPos(); end
gt(:,5)=0; [gtRes,dtRes]=bbGt('evalRes',gt,double(bbs),.7);
figure(1); bbGt('showRes',I,gtRes,dtRes(dtRes(:,6)==1,:));
title('green=matched gt  red=missed gt  dashed-green=matched detect');

%% run and evaluate on entire dataset (see boxesData.m and boxesEval.m)
if(~exist('boxes/VOCdevkit/','dir')), return; end
split='val'; data=boxesData('split',split);
nm='EdgeBoxes70'; opts.name=['boxes/' nm '-' split '.mat'];
edgeBoxes(data.imgs,model,opts); opts.name=[];
boxesEval('data',data,'names',nm,'thrs',.7,'show',2);
boxesEval('data',data,'names',nm,'thrs',.5:.05:1,'cnts',1000,'show',3);

author: pdollar

repo: edges

save_path: github-repos/MATLAB/pdollar-edges