text:
function nmfmse( V, rdim, fname, showflag )
%

% Check that we have non-negative data
if min(V(:))<0, error('Negative values in data!'); end

% Globally rescale data to avoid potential overflow/underflow
V = V/max(V(:));

% Dimensions
vdim = size(V,1);
samples = size(V,2);

% Create initial matrices
W = abs(randn(vdim,rdim));
H = abs(randn(rdim,samples));

% Initialize displays
if showflag,
   figure(1); clf; % this will show the energies and sparsenesses
   figure(2); clf; % this will show the objective function
   drawnow;
end

% Calculate initial objective
objhistory = 0.5*sum(sum((V-W*H).^2));

timestarted = clock;

% Start iteration
iter = 0;
while 1,

    % Show progress
    fprintf('[%d]: %.5f \n',iter,objhistory(end));

    % Save every once in a while
    if rem(iter,5)==0,
	elapsed = etime(clock,timestarted);
	fprintf('Saving...');
	save(fname,'W','H','iter','objhistory','elapsed');
	fprintf('Done!\n');
    end

    % Show stats
    if showflag & (rem(iter,5)==0),
	figure(1);
	cursW = (sqrt(vdim)-(sum(W)./sqrt(sum(W.^2))))/(sqrt(vdim)-1);
	cursH = (sqrt(samples)-(sum(H')./sqrt(sum(H'.^2))))/(sqrt(samples)-1);
	subplot(3,1,1); bar(sqrt(sum(W.^2)));
	subplot(3,1,2); bar(cursW);
	subplot(3,1,3); bar(cursH);
	if iter>1,
	    figure(2);
	    plot(objhistory(2:end));
	end
	drawnow;
    end

    % Update iteration count
    iter = iter+1;

    % Save old values
    Wold = W;
    Hold = H;

    % Compute new W and H (Lee and Seung; NIPS*2000)
    H = H.*(W'*V)./(W'*W*H + 1e-9);
    W = W.*(V*H')./(W*H*H' + 1e-9);

    % Renormalize so rows of H have constant energy
    norms = sqrt(sum(H'.^2));
    H = H./(norms'*ones(1,samples));
    W = W.*(ones(vdim,1)*norms);

    % Calculate objective
    newobj = 0.5*sum(sum((V-W*H).^2));
    objhistory = [objhistory newobj];

end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function u=randDirVec(numDim,N)
%%RANDDIRVEC Generate a uniformly distributed random direction vector
%            (unit vector) in an arbitrary number of dimensions.
%
%INPUTS: numDim The number of dimensions of the random unit vector. This
%               is >=1.
%             N The number of random direction vectors to generate. If
%               this parameter is omitted, then N=1 is used.
%
%OUTPUTS: u A numDimXN matrix of N numDim-dimensional unit vectors that
%           are uniformly distributed on the unit-sphere (or hypersphere).
%
%One cannot simply generate a uniformly-distributed direction vector by
%setting each dimensions to a Uniform(-1,1) random variable and then
%normalizing the result. Rather, following the algorithm in Chapter 3.4.1,
%E6 of [1], the elements of the vector must be generated as normal 0-1
%random variables and then the vector normalized.
%
%REFERENCES:
%[1] D. Knuth, The Art of Computer Programming: Seminumerical Algorithms,
%    3rd ed. Reading, MA: Addison-Wesley, 1998, vol. 2.
%
%September 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2)
    N=1;
end

u=zeros(numDim,N);

for curN=1:N
    %Knuth, Chapter 3.4.1, E6 (Random point on an n-dimensional sphere with
    %radius one.
    u(:,curN)=randn(numDim,1);
    u(:,curN)=u(:,curN)/norm(u(:,curN));

    %Deal with the instance that norm(u)=0, which will occur with an extremely
    %small probability due to the limited precision of the random number
    %generation.
    if(~isfinite(u(:,curN)))
       u(:,curN)=zeros(numDim,1);
       u(1,curN)=1;
    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 r8vec_legendre_test ( )

%*****************************************************************************80
%
%% R8VEC_LEGENDRE_TEST tests R8VEC_LEGENDRE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_LEGENDRE_TEST\n' );
  fprintf ( 1, '  R8VEC_LEGENDRE computes N Legendre points in [R1,R2].\n' );

  r1 = -1.0;
  r2 = +1.0;
  n = 5;

  r = r8vec_legendre ( n, r1, r2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  N = %d, R1 = %g, R2 = %g\n', n, r1, r2 );

  r8vec_print ( n, r, '  Legendre points:' );

  r1 =   0.0;
  r2 = +10.0;
  n = 7;

  r = r8vec_legendre ( n, r1, r2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  N = %d, R1 = %g, R2 = %g\n', n, r1, r2 );

  r8vec_print ( n, r, '  Legendre points:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function quadrule_test13 ( )

%*****************************************************************************80
%
%% TEST13 tests LAGUERRE_EK_COMPUTE and LAGUERRE_SUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2011
%
%  Author:
%
%    John Burkardt
%
  order_max = 20;

  nfunc = func_set ( 'COUNT', 'DUMMY' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST13\n' );
  fprintf ( 1, '  LAGUERRE_EK_COMPUTE computes a Gauss-Laguerre rule;\n' );
  fprintf ( 1, '  LAGUERRE_SUM carries it out.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Quadrature order will vary.\n' );
  fprintf ( 1, '  Integrand will vary.\n' );
  fprintf ( 1, '  The weight function is EXP ( - X ).\n' );
  fprintf ( 1, '\n' );

  a = 1.0;

  fprintf ( 1, '  The integration interval is [ %f, +oo ).\n', a );
  fprintf ( 1, '\n' );

  for ilo = 1 : 5 : nfunc

    ihi = min ( ilo + 4, nfunc );

    fprintf ( 1, '\n' );
    fprintf ( 1, '    ' );
    for i = ilo : ihi
      fprintf ( '%14s', fname(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );

    for norder = 1 : order_max

      fprintf ( 1, '  %2d', norder );

      for i = ilo : ihi

        func_set ( 'SET', i );

        [ xtab, weight ] = laguerre_ek_compute ( norder );

        result(i) = laguerre_sum ( @func, a, norder, xtab, weight );

        fprintf ( 1, '  %12f', result(i) );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function radDepthV = getRadiologicalDepth(xV, yV, zV, gantryAngle, isocenter, isodistance, structNum, scanSet, planC)
%"getRadiologicalDepth"
%   Returns the radiological depth at the points in xV, yV, zV, as
%   determined from scanSet with a point source at sourceV.  structNum
%   tells the interpolation routine to only include the densities of
%   points within a certain structure when calculating the radiological
%   depth.
%
%JRA 06/16/05
%
%Usage:
%  function radDepthV = getRadiologicalDepth(xV, yV, zV, gA, iC, iD, structNum, scanSet, planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC')
    global planC
end
indexS = planC{end};

siz = size(xV);

%Get mask of requested structure.
rS = getRasterSegments(structNum, planC);
[structMask, slicesV] = rasterToMask(rS, scanSet, planC);
% rS = getRasterSegments(structNum, planC);
% [structMask, slicesV] = rasterToMask(rS, scanSet, planC);
% [surfPoints] = getSurfacePoints(structMask);

%Get coordinates of scan voxels inside the structure.
[xS, yS, zS] = getScanXYZVals(planC{indexS.scan}(scanSet));
[xM, yM, zM] = meshgrid(xS, yS, zS(slicesV));
xM = xM(structMask);
yM = yM(structMask);
zM = zM(structMask);
% [xS, yS, zS] = getScanXYZVals(planC{indexS.scan}(scanSet));
% zS = zS(slicesV);
% xM = xS(surfPoints(:,1));
% yM = yS(surfPoints(:,2));
% zM = zS(surfPoints(:,3));
% clear surfPoints

%Convert points to collimator coordinates.
coll3V = scan2Collimator([xM(:) yM(:) zM(:)], gantryAngle, 0, 0, isocenter, isodistance);
clear xM yM zM;

%Get distance from source to all points.
distsquared = sepsq(coll3V', [0 0 0]');

%Project vector to points to the perpendicular plane at isocenter.
coll3V = coll3V./repmat(coll3V(:,3), [1 3]) * -isodistance;

%Find the extents of X,Y,D for the region covering the structure (ie, all
%nonzero scan numbers since we will set scan densities outside of structure
%to be 0.
minExt = min(coll3V);
maxExt = max(coll3V);
clear coll3V;

minProjX = minExt(1);
maxProjX = maxExt(1);

minProjY = minExt(2);
maxProjY = maxExt(2);

maxDsq = max(distsquared);
minDsq = min(distsquared);
clear distsquared;

%Set points in scanArray outside of structure to 0.
sA = planC{indexS.scan}(scanSet).scanArray;
for i=1:size(sA, 3);
    [ismem, loc] = ismember(i, slicesV);
    if ismem
        slc = sA(:,:,i);
        slc(~structMask(:,:,loc)) = 0;
        sA(:,:,i) = slc;
    else
        sA(:,:,i) = zeros(size(sA(:,:,i)));
    end
end
clear structMask;

%Convert requested points to collimator coordinates.
coll3V = scan2Collimator([xV(:) yV(:) zV(:)], gantryAngle, 0, 0, isocenter, isodistance);
clear xM yM zM;

%Get distance from source to all requested points.
distsquared = sepsq(coll3V', [0 0 0]');

%Project vector to points to the perpendicular plane at isocenter.
coll3V = coll3V./repmat(coll3V(:,3), [1 3]) * -isodistance;

clear xV yV zV

%Find out of bounds points, for which calculation is not necessary,
%radDepth is zero.
OOBPoints = coll3V(:,1) < minProjX | coll3V(:,1) > maxProjX | coll3V(:,2) < minProjY | coll3V(:,2) > maxProjY;

%Clip the r values so they take the density of the nearest calculated
%point (points less than minDsq would take zero, and greater than maxDsq
%would take the last radDepthV value along that ray.
r = clip(distsquared, minDsq, maxDsq, 'limits');

%Set the resolution of the radiological depth calculation (at isocenter).
dX = .5;
dY = .5;
dR = .25;

%Mesh the radiological depth calculation grid in collimator coordinates.
projXV = minProjX:dX:maxProjX;
projYV = minProjY:dY:maxProjY;
distV = sqrt(minDsq):dR:sqrt(maxDsq);
[projXVM, projYVM, distVM] = meshgrid(projXV, projYV, distV);

projDistV = sqrt(sepsq([projXVM(:) projYVM(:) -isodistance*ones(size(projXVM(:)))]', [0 0 0]'));
projXVM = projXVM(:) .* (distVM(:) ./ projDistV);
projYVM = projYVM(:) .* (distVM(:) ./ projDistV);
projZVM = -isodistance*ones(size(distVM(:))) .* (distVM(:) ./ projDistV);

%Convert the depth calculation grid to rectangular coordinates.
newXYZ = collimator2Scan([projXVM(:) projYVM(:) projZVM(:)], gantryAngle, 0, 0, isocenter, isodistance);
% clear projXVM projYVM distVM

%Get the scan values at each point in the radiological depth grid.
[xVS, yVS, zVS] = getScanXYZVals(planC{indexS.scan}(scanSet));

%Interpolate to values in xV/yV/zV from the scanarray, 0 if out of bounds.
scanVals = finterp3(newXYZ(:,1), newXYZ(:,2), newXYZ(:,3), sA, [xVS(1) xVS(2)-xVS(1) xVS(end)], [yVS(1) yVS(2)-yVS(1) yVS(end)], zVS, 0);

scanVals = reshape(scanVals, length(projXV), length(projYV), length(distV));

%Scale values by the resolution of R to get the linear radiological depth.
%Divide by 1024 to get the radiological depth in units of water.
scanVals = scanVals .* (dR / 1024);

%Take the cumsum to get the actual radiological depth at each point in the
%depth calculation grid.
cumDens = cumsum(scanVals, 3);
clear scanVals

cumDens = reshape(cumDens, size(distVM));

%Interpolate to values in xV/yV/zV from the scanarray, 0 if out of bounds.
radDepthV = finterp3(coll3V(:,1), coll3V(:,2), sqrt(distsquared), cumDens, [projXV(1) projXV(2)-projXV(1) projXV(end)], [projYV(1) projYV(2)-projYV(1) projYV(end)], distV, 0);

radDepthV(OOBPoints) = 0;
radDepthV = reshape(radDepthV, siz);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function st = type( x )

s = size( x );
len = prod( s );
if len == 1,
    nzs = len;
    st = 'scalar';
else
    nzs = nnz( any( cvx_basis( x ), 1 ) );
    nd = length( s );
    st = sprintf( '%dx', s );
    st = st( 1 : end - 1 );
    if nd > 2,
        st = [ st, ' array' ];
    elseif any( s == len ),
        st = [ st, ' vector' ];
    else
        st = [ st, ' matrix' ];
    end
end
if nzs < len,
    if nzs > 1,
        st = sprintf( '%s, %d nonzeros', st, nzs );
    elseif nzs == 1,
        st = sprintf( '%s, 1 nonzero', st );
    end
end

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

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function AIFPScaled = refRegionScale(AIFP,inputS,shiftS,paramFile,headerS,muscleTimeCourse3M,muscleMaskM)
% Function to pull muscle time course from ROI mask & convert to concentration
% Loads average AIF for the study (or individual), runs Tofts fit over range of AIF amplitudes
% to find AIF amplitude that gives Tofts fit ve = 0.1 for Muscle time course data.
% Returns this AIF time course.
%
% Kristen Zakian
% AI 7/06/16
% AI 7/21/16  Updated to write ROI signal shift to parameter file

% Get parameters
T10Muscle_ms = inputS.T10ref;
TRms = headerS.TRms;
tDel = headerS.tDel ;
FA = headerS.FA;
frames = headerS.frames;
r1 = 3.9*(headerS.B0 < 2.0)*strcmp(inputS.cAgent,'gd')+ 3.3*(headerS.B0 > 2.0)*strcmp(inputS.cAgent,'gd') + ...
    8.1*(headerS.B0 < 2.0)*strcmp(inputS.cAgent,'mh') + 6.3*(headerS.B0 > 2.0)*strcmp(inputS.cAgent,'mh');

% Generate vector of sampling times  (1 x frames size)
ind = 0:(frames-1);
tsec = double(ind)*tDel;
tmin = tsec./60.;                      %ktrans, etc. are in /min.  Parker coefs are in min or /min.  Tofts fit will use minutes.

% Generate an average time course for the muscle ROI from the slice
avgROISigV = zeros(1,20);
for t = 1:frames
    muscleMaskSum = sum(muscleTimeCourse3M(:,:,t));
    avgROISigV(t) = sum(muscleMaskSum(:))/nnz(muscleMaskM);         % avg of masked elements
end

%  Display for deciding muscle ROI point shift. Must be shifted to zero for Tofts fit

%Plot the Muscle ROI time course
figTitle = 'Start of muscle time course';
figure ('Name', figTitle);
plot (tmin,avgROISigV,'-d');
pause(2);

%Pass AIFAmp from parameter file if available
%(AIFAmp is the value that gives ve < 0.1)
useValue = 'No';
shiftC = fieldnames(shiftS);
shiftIn = any(strcmp(shiftC,'muscleShift'));
if shiftIn
    if ~isempty(shiftS.muscleShift)
        aShiftIn = shiftS.muscleShift;
        useValue = questdlg(sprintf('Use value from parameter file:\n muscleShift = %g ?',aShiftIn),...
            'AIF scaling','Yes','No','Yes');
    end
end

if strcmp(useValue,'No')
    aShift = input('\nInput the muscle point shift (e.g. 2, to shift left 2 points):   ');
    %Store shift to parameter file
    fid = fopen(paramFile);
    inputC = textscan(fid, '%s %s %n','endofline','\r\n');
    fmt = '\r\n%s\t%s\t%d';
    fclose(fid);
    if shiftIn
        shiftIdx = strcmp (inputC{1},'muscleShift');
        inputC{3}(shiftIdx) = aShift;
        fid = fopen(paramFile,'w+');
        for lineNum = 1:size(inputC{1},1)
            col1 = inputC{1}(lineNum);
            col2 = inputC{2}(lineNum);
            col3 = inputC{3}(lineNum);
            fprintf(fid,fmt,col1{1},col2{1},col3);
        end
        fclose(fid);
    else
        fid = fopen(paramFile,'a');
        newEntryC = {'muscleShift','na',aShift};
        fprintf(fid,fmt,newEntryC{:});
        fclose(fid);
    end
else
    aShift = aShiftIn;
end
close(gcf);
muscleShiftV = circshift(avgROISigV', -aShift);
rawMuscleSigV = muscleShiftV';
shiftRawMuscleSigV = rawMuscleSigV(frames-aShift-8:frames-aShift);
meanRawMuscle = mean(shiftRawMuscleSigV(:));
rawMuscleSigV(frames-aShift+1:frames)= meanRawMuscle;
baseptsMuscle = aShift-1; % set number of baseline points = number of shifted points -1

% Display unshifted and shifted aif curve
figtitle = ['Muscle time course unshifted and shifted. Shift = ',num2str(aShift)];
figure ('Name', figtitle);
plot (tmin,avgROISigV,tmin,rawMuscleSigV,'--');
pause(4);
close(gcf);

% Convert ROI average time course to Concentration
% (need muscle pre-contrast T1 and average baseline value S0)
TR = TRms/1000.;                    %(TR in DICOM header is in ms. Change to s)
T10 = T10Muscle_ms/1000;
base = avgROISigV(1:baseptsMuscle); % baseline signal value obtained from time course BEFORE SHIFT
s0 = mean(base);
pCalcConcMuscle = [r1, TR, FA, T10];
muscleConcTCourseV = calcConc(rawMuscleSigV,pCalcConcMuscle,s0);
figtitle = ('Muscle concentration time course ');
figure ('Name', figtitle);
plot (tmin,muscleConcTCourseV);
pause(1);
close(gcf);

%  Fit muscle time course to Tofts model until ve = 0.1.
%Start with average AIF.
%Vary the amplitude of the average AIF until you get a ve of 0.1.
%That will then be the AIF that you use.

inputsC = fieldnames(inputS);                         %Pass AIFAmp from parameter file if available
if any(strcmp(inputsC,'AIFAmp'))                      %(AIFAmp is the value that gives ve < 0.1)
    if ~isempty(inputS.AIFAmp)
        AIFPScaled = refRegionToftsAIFScale(tmin,AIFP,muscleConcTCourseV,paramFile,inputS.AIFAmp);
    end
else
    AIFPScaled = refRegionToftsAIFScale(tmin,AIFP,muscleConcTCourseV,paramFile);
end
%Display figure
figtitle = 'AIF plasma and AIF scaled ';
figure ('Name', figtitle);
plot (tmin,AIFP,tmin,AIFPScaled,'--');
pause(1);
close(gcf);

end


text:
function [splitPsi, logQ] = sampleSplitConfig( Psi, data, anchorObjIDs, featIDs, algParams, TargetPsi )
% Sample a proposal config of (F,z) that has one extra state than input
%INPUT:
%  Psi : input model config (contains F,z)
%  data : data object
%  anchorObjIDs : length 2 vector giving ids of anchor objects for proposal
%  featIDs : integer id of feature in Psi to split into two
%  TargetPsi : target model config.
%    when provided, we don't actually sample, but instead calculate
%    probability "logQ" of sampling the target given the input Psi config
%OUTPUT:
%  splitPsi : output configuration (proposal)
%    has all its observed data sufficent statistics pre-calculated
%    for easy use in later calculations (accept ratios, etc.
%  logQ     : struct of transition probabilities (in log space)

% ------------------------------------------------------- INTERP INPUT
anchorObjIDs = anchorObjIDs(:)';
ii = anchorObjIDs(1);
jj = anchorObjIDs(2);

kold = featIDs(1);
otherObjIDs = setdiff( find( Psi.F( :, kold ) > 0   )', anchorObjIDs );
activeObjIDs = [anchorObjIDs otherObjIDs];

if exist( 'TargetPsi', 'var' ) && ~isempty( 'TargetPsi' )
   kA = TargetPsi.activeFeatIDs(1);
   kB = TargetPsi.activeFeatIDs(2);
else
    TargetPsi = [];
    kA = kold;
    kB = size(Psi.F,2) + 1;
end

% ------------------------------------------------------- INIT PROPOSAL
propF = Psi.F == 1;
propF(:, kold) = 0;
propF( ii, kA) = 1;
propF( jj, kB) = 1;
propFeatIDs = [kA kB];

propStateSeq = Psi.stateSeq;
for aa = 1:data.N
    propStateSeq( aa ).z( Psi.stateSeq(aa).z == kold ) = 0;
end
propStateSeq( ii ).z(  Psi.stateSeq( ii ).z == kold )  = kA;
propStateSeq( jj ).z(  Psi.stateSeq( jj ).z == kold )  = kB;

% Create deterministic trans params Eta
propTransM = Psi.TransM.getAllEta_PriorMean( activeObjIDs, propF, propFeatIDs );

% Create deterministic Theta for ALL features: others + [kA kB]
propThetaM = Psi.ThetaM.getAllTheta_PosteriorMean( data, propStateSeq, propFeatIDs, featIDs);
assert( propThetaM.K == length(propThetaM.Xstats), 'Badness');

% Init suff stats for F
featCounts = sum( propF(:,propFeatIDs) , 1 );
nObj = size(propF, 1) - length(otherObjIDs);
logQ.F = 0;
logQ.z = 0;

proposalsON = false(1, size(propF,2) );
proposalsON(propFeatIDs) = 1;
for aa = [ otherObjIDs( randperm(length(otherObjIDs)) ) ii jj]

    if aa == ii || aa == jj
        % Make sure to update F suff. stats for anchors
        %  since we've already "seen" them before, unlike other sequences
        featCounts = featCounts - propF(aa, propFeatIDs);
        nObj = nObj - 1;
        propThetaM = propThetaM.decXStats( aa, data, propStateSeq, propFeatIDs );
    end

    % ====================================================  Sample F(aa,:)
    % Calc Soft Evidence for the current sequence, including both new feats
    allActiveFeatIDs = find( propF(aa,:) | proposalsON );
    seqSoftEv = propThetaM.calcLogSoftEv( aa, data, allActiveFeatIDs );

    % Propose new feature assignments for the cur. seq.
    priorPrFeats = featCounts ./ ( nObj + Psi.bpM.c );
    [newFaa, logQ_Faa] = sampleTwoSharedFeatsForSeq_Gibbs( aa, propFeatIDs, seqSoftEv, propTransM.seq(aa), priorPrFeats, anchorObjIDs, TargetPsi);
    propF(aa, propFeatIDs) = newFaa==1;
    logQ.F = logQ.F + logQ_Faa;

    % ===============================================  Sample stateSeq(aa)
    seqTransM = propTransM.getSeqEtaWithFeats( aa, propF(aa,:) );
    propTransM = propTransM.setEta( aa, propF(aa,:), seqTransM.eta );
    [propStateSeq(aa).z, logQ_zaa] = sampleSingleStateSeq_WithSoftEv( aa, seqTransM, seqSoftEv(propF(aa,:), :), TargetPsi);
    logQ.z = logQ.z + logQ_zaa;

    % Update suff. stats for F
    featCounts = featCounts + propF(aa, propFeatIDs);
    nObj = nObj + 1;

    % Update Emission Params for new proposed features [kA, kB]
    %   using the newly allocated data from the current sequence
    propThetaM = propThetaM.incXStats( aa, data, propStateSeq, propFeatIDs );
    for kkID = 1:length(propFeatIDs)
        kk = propFeatIDs(kkID);
        PN = propThetaM.getPosteriorParams( propThetaM.Xstats(kk) );
        propThetaM.theta(kk) = propThetaM.getTheta_Mean( PN );
    end
end

if isempty( TargetPsi )
    propThetaM = propThetaM.updateAllXSuffStats( horzcat(propStateSeq(:).z), data );
end

splitPsi.F = propF;
splitPsi.stateSeq = propStateSeq;
splitPsi.ThetaM = propThetaM;
splitPsi.TransM = propTransM;
splitPsi.bpM = Psi.bpM;
splitPsi.activeFeatIDs = propFeatIDs;
splitPsi.anchorIDs  = anchorObjIDs;

author: michaelchughes

repo: NPBayesHMM


text:
function nword = file_word_count ( file_name )

%*****************************************************************************80
%
%% FILE_WORD_COUNT counts the number of words in a file.
%
%  Discussion:
%
%    The file is assumed to be a simple text file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the name of the file.
%
%    Output, integer NWORD, the number of words found in the file.
%
  nword = 0;
%
%  Open the file.
%
  fid = fopen ( file_name );

  if ( fid = -1 )
    nword = -1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_WORD_COUNT - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file:\n' );
    fprintf ( 1, '"%s".\n', file_name );
    return
  end
%
%  Read the lines.
%
  while ( 1 )

    line = fgets ( fid );

    if ( line == -1 )
      break;
    end

    nplus = s_word_count ( line );

    nword = nword + nplus;

  end

  fclose ( fid )

  return
end


text:
function problem = EmptyProblem()
%EmptyProblem - Create an empty problem
% Create an empty problem with all the required fields

% Basic elements (these are mandatory)
problem.type= 'none'; % What kind of problem is this?
problem.N = 0; % How many measurements this problem has?
problem.outliers = []; % Among the measurements, what are the indices of the outliers?
problem.priors = [];  % Amont the measurements, what are the indices of the element that will never be outliers?
problem.dof = 0; % Degrees of freedom of the noise process

% Problem specific data
% Add here any problem specific field ...

problem.time = now;
% Unique identifier
problem.uuid = java.util.UUID.randomUUID.toString;
end

author: MIT-SPARK

repo: GNC-and-ADAPT


text:
function addAnnotationTracks(infile)

Q = load(infile);
Q.trx(2) = Q.trx(1);
Q.trx(3) = Q.trx(1);
Q.trx(1).x(:) = Q.trx(1).arena.food(1);
Q.trx(1).y(:) = Q.trx(1).arena.food(2);
Q.trx(2).x(:) = Q.trx(1).arena.mouth(1);
Q.trx(2).y(:) = Q.trx(1).arena.mouth(2);
Q.trx(3).x(:) = Q.trx(1).arena.perch(1);
Q.trx(3).y(:) = Q.trx(1).arena.perch(2);

save(infile,'-struct','Q');

author: kristinbranson

repo: JAABA


text:
% system process model in local s0 frame with dcm formulation, Cs0e is the
% transformation from the sensor local s0 frame to e-frame
% written by Huai
% rqs02e, r s0 in e, q s0 2e, rvqs0, rs in s0, vs in s0, qs0 2s
% acc, m/s^2, acc by imu in s frame, gyro, rad/s, angular rate by imu
% w.r.t i frame coordinated in s frame, dt, time interval for covariance update
% modelNo, imu model for bias and scale factor errors
function [STM Qd]=sys_local_dcm_v001(rqs02e, rvqs0, acc, gyro, dt, imutype, modelNo)
% the covariance corresponds to states,
% rs in s0, v s in s0, q s02s, gravity in s0, ba, bg, sa, sg, qs2c, Ts in c
% for modelNo=3 and imutype=5, MEMS
% acc bias drift, random walk
% gyro bias drift, random walk
% acc scale factor, random walk
% gyro scale factor, random walk
% for modelNo=1 the acc and gyro turn on bias are removed

%system disturbance coefs
navStates=12;
N=zeros(navStates,6);
N(7:9,4:6)=eye(3); %attitude
Cs02s=quat2dcm_v000(rvqs0(7:10));
N(4:6,1:3)=Cs02s'; %velocity

%system matrix
A=zeros(navStates);
A(1:3, 4:6)=eye(3); % rs in s0
%Velocity
WIE_E=7292115e-11;
wie2s0=quatrot_v000(rqs02e(4:7),[0;0;WIE_E], 1);
% A(4:6,1:3)=-skew(wie2s0)^2; % for low grade IMU, w_ie^e is buried
% $\delta\bar{g}^{s0}(\mathbf{x}_s^{s0})$ is noisy, but buried by low grade
% IMU noise
% A(4:6,4:6)=-2*skew(wie2s0);
A(4:6,7:9)=-Cs02s'*skew(acc);
% A(4:6,10:12)=eye(3); % for low grade IMU, w_ie^e is buried

%Attitude
A(7:9,7:9)=skew(-gyro);
Anav=A;
Nnav=N;
% X(k+1) = ffun[X(k),U(k),V(k)]
% X(k+1) = Ak*X(k)+Gk*Vk

%%%%Imu error model parameters
[Aimu_d, Qimu_d, Cimu, Rimu]=imu_err_model_v001(acc, gyro, dt, imutype, modelNo);
%%%%Combine and discretize nav and imu models
% this discretization can also be accomplished by Loan's matrix exponential
% method, see sys_metric_phipsi_v000.m

Anav_d=eye(navStates)+dt*Anav;  %Use 1st order taylor series to discretize Anav
Qnav=Nnav*Rimu*Nnav';
Qnav_d=dt/2*(Anav_d*Qnav+Qnav*Anav_d');      %Use trapezoidal rule to discretize Rimu

STM=zeros(navStates+size(Aimu_d,1));
STM(1:navStates,1:navStates)=Anav_d;
STM(1:navStates,1+navStates:end)=Nnav*Cimu*dt;
STM(1+navStates:end,1+navStates:end)=Aimu_d;

Qd=zeros(navStates+size(Aimu_d,1));
Qd(1:navStates,1:navStates)=Qnav_d;
Qd(1+navStates:end,1+navStates:end)=Qimu_d;
Qd(1:navStates,1+navStates:end)=Nnav*Cimu*Qimu_d*dt/2;
Qd(1+navStates:end,1:navStates)=Qd(1:navStates,1+navStates:end)';
end

repo: ekfmonoslam


text:
function [maxel,IJ]= max2(M,userows,usecols)
% finds the location of the single overall maximum element in a 2-d array
% usage: [maxel,IJ] = max2(M)
% usage: [maxel,IJ] = max2(M,userows,usecols)
%
% The location in a 2-d array of the overall
% maximum element (or the first incidence of
% several, if the maximum is not unique), where
% you may restrict the search to a set of
% specified rows and/or columns.
%
% Note that max2 does NOT convert the matrix to
% linear indexing, so that really huge arrays
% can be worked with.
%
% arguments: (input)
%  M - an (nxm) 2-dimensional numeric array (or
%      vector) that max is able to operate on. M
%      may contain inf or -inf elements.
%
%  userows - (OPTIONAL) a list of the rows to be
%      searched for the maximum. The search will
%      be restricted to this set of rows. If empty.
%      there will be no row restriction.
%
%      userows must be a list of integers
%
%  usecols - (OPTIONAL) a list of the columns to be
%      searched for the maximum. The search will
%      be restricted to this set of columnss. If
%      empty. there will be no column restriction.
%
% arguments: (output)
%  maxel - overall maximum element found. If the
%      maximum was ot unique, then this is the
%      first element identified. Ties will be
%      resolved in a way consistent with find.
%
%  IJ - a (1x2) row vector, comtaining respectively
%      the row and column indices of the maximum as
%      found.
%
% Example:
%  M = magic(4)
% ans =
%    16     2     3    13
%     5    11    10     8
%     9     7     6    12
%     4    14    15     1
%
% % the overall maximum
%  [maxel,IJ] = max2(M)
% maxel =
%      16
% IJ =
%      1     1
%
%
% % a restricted maximum
%  [maxel,IJ] = max2(M,[1 2 3],[2 3])
% maxel =
%      11
% IJ =
%     2     2
%
%
% See also: max2, max, min, find
%
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 2/16/09

% check the arguments
if (nargin<1) || (nargin>3)
  error('max2 may have 1, 2, or 3 arguments only')
end

if length(size(M)) > 2
  error('M must be a 2-d array or a vector')
end
[n,m] = size(M);

% default for userows?
if (nargin<2) || isempty(userows)
  userows = 1:n;
else
  userows = unique(userows);
  if ~isnumeric(userows) || any(diff(userows)==0) || ...
      any(userows<1) || any(userows>n) || any(userows~=round(userows))
    error('userows must be a valid set of indices into the rows of M')
  end
end

% default for usecols?
if (nargin<3) || isempty(usecols)
  usecols = 1:m;
else
  usecols = unique(usecols);
  if ~isnumeric(usecols) || any(diff(usecols)==0) || ...
      any(usecols<1) || any(usecols>m) || any(usecols~=round(usecols))
    error('usecols must be a valid set of indices into the columns of M')
  end
end

% restrict the search
Muse = M(userows,usecols);

% The maximum down the rows
[maxrows,rowind] = max(Muse,[],1);

% find the best of these maxima
% across the columns
[maxel,colind] = max(maxrows,[],2);
rowind = rowind(colind);

% package the row and column indices
% together, in terms of the original
% matrix in case there was a restiction.
IJ = [userows(rowind),usecols(colind)];

author: Sable

repo: mcbench-benchmarks


text:
%+========================================================================+
%|                                                                        |
%|            This script uses the GYPSILAB toolbox for Matlab            |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             francois.alouges@polytechnique.edu                         |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : nrtDomND.m                                    |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Multi-dimensional quadrature                  |
%|  `---'  |                                                              |
%+========================================================================+

% Cleaning
clear all
close all
clc

% Gypsilab path
run('../../addpathGypsilab.m')

% Tetra mesh
Nvtx = 1e2;
L    = [5 4 3];
mesh = mshCube(Nvtx,L);

% Boundary mesh
bound = mesh.bnd;
ctr   = bound.ctr;
bound.col(ctr(:,3)==0.5*L(3)) = 1;

% Upperface
upperface = bound.sub(bound.col==1);
rectangle = upperface.bnd;

% Domain
omega = dom(mesh,5);
figure
plot(mesh)
hold on
plot(omega)
alpha(0.2)
axis equal
view(30,30)

% Boundary domain
sigma = dom(bound,3);
figure
plot(bound)
hold on
plot(sigma,'.y')
plotNrm(sigma,'y')
alpha(0.5)
axis equal
view(30,30)

% Rectangle
gamma = dom(rectangle,3);
figure
plot(rectangle)
hold on
plot(gamma)
alpha(0.5)
axis equal
view(30,30)

% Numerical integration on the cube
fct = @(X) ones(size(X,1),1);
V   = integral(omega,fct);
norm(V-prod(L),'inf')

% Numerical integration on the bundary of the cube
fct = @(X) ones(size(X,1),1);
S   = integral(sigma,fct);
norm(S-2*(L(1)*L(2)+L(2)*L(3)+L(1)*L(3)),'inf')

% Numerical integration on the rectangle
fct = @(X) ones(size(X,1),1);
P   = integral(gamma,fct);
norm(P-2*(L(1)+L(2)),'inf')

disp('~~> Michto gypsilab !')

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%% QUAD_SCRIPT is a script to call QUAD_FUN.
%
%  Discussion:
%
%    The BATCH command runs scripts, not functions.  So we have to write
%    this short script if we want to work with BATCH!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2010
%
%  Author:
%
%    John Burkardt
%
  n = 100000;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUAD_SCRIPT\n' );
  fprintf ( 1, '  Estimate integral using %d points\n', n );

  value = quad_fun ( n );


text:
function f = testfcn1(x)

if strcmpi(x,'init')
    f.options.PopInitRange = [0, 0; 4, 4] ;
    f.options.Vectorized = 'on' ;
    f.options.HybridFcn = {} ;
    f.options.Generations = 500 ;
    f.options.ConstrBoundary = 'penalize' ;
    f.LB = [0,0] ; f.UB = [] ;
    f.Aeq = [] ; f.beq = [] ;
    f.Aineq = [-1 0] ; f.bineq = [-0.01] ;
%     f.Aineq = [] ; f.bineq = [] ;
    f.nonlcon = [] ;
else
    a = 1 ; b = 0.01 ;
    f = (x(:,1).^b + x(:,2).^(1-b)).^a ;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function Serr=specerr(S,J,err,trialave,numsp)
% Function to compute lower and upper confidence intervals on the spectrum
% Usage: Serr=specerr(S,J,err,trialave,numsp)
% Outputs: Serr (Serr(1,...) - lower confidence level, Serr(2,...) upper confidence level)
%
% Inputs:
% S - spectrum
% J - tapered fourier transforms
% err - [errtype p] (errtype=1 - asymptotic estimates; errchk=2 - Jackknife estimates;
%                   p - p value for error estimates)
% trialave - 0: no averaging over trials/channels
%            1 : perform trial averaging
% numsp    - number of spikes in each channel. specify only when finite
%            size correction required (and of course, only for point
%            process data)
%
% Outputs:
% Serr - error estimates. Only for err(1)>=1. If err=[1 p] or [2 p] Serr(...,1) and Serr(...,2)
% contain the lower and upper error bars with the specified method.
if nargin < 4; error('Need at least 4 input arguments'); end;
if err(1)==0; error('Need err=[1 p] or [2 p] for error bar calculation. Make sure you are not asking for the output of Serr'); end;
[nf,K,C]=size(J);
errchk=err(1);
p=err(2);
pp=1-p/2;
qq=1-pp;

if trialave
   dim=K*C;
   C=1;
   dof=2*dim;
   if nargin==5; dof = fix(1/(1/dof + 1/(2*sum(numsp)))); end
   J=reshape(J,nf,dim);
else
   dim=K;
   dof=2*dim*ones(1,C);
   for ch=1:C;
     if nargin==5; dof(ch) = fix(1/(1/dof + 1/(2*numsp(ch)))); end
   end;
end;
Serr=zeros(2,nf,C);
if errchk==1;
   Qp=chi2inv(pp,dof);
   Qq=chi2inv(qq,dof);
   Serr(1,:,:)=dof(ones(nf,1),:).*S./Qp(ones(nf,1),:);
   Serr(2,:,:)=dof(ones(nf,1),:).*S./Qq(ones(nf,1),:);
elseif errchk==2;
   tcrit=tinv(pp,dim-1);
   for k=1:dim;
       indices=setdiff(1:dim,k);
       Jjk=J(:,indices,:); % 1-drop projection
       eJjk=squeeze(sum(Jjk.*conj(Jjk),2));
       Sjk(k,:,:)=eJjk/(dim-1); % 1-drop spectrum
   end;
   sigma=sqrt(dim-1)*squeeze(std(log(Sjk),1,1)); if C==1; sigma=sigma'; end;
   conf=repmat(tcrit,nf,C).*sigma;
   conf=squeeze(conf);
   Serr(1,:,:)=S.*exp(-conf); Serr(2,:,:)=S.*exp(conf);
end;
Serr=squeeze(Serr);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function varargout = ddgaussian(varargin)
% VL_DDGAUSSIAN Second derivative of the Gaussian density function
%   Y=VL_DDGAUSSIAN(X) computes the second derivative of the standard
%   Gaussian density.
%
%   To obtain the second derivative of the Gaussian density of
%   standard deviation S, do
%
%     Y = 1/S^3 * VL_DDGAUSSIAN(X/S) .
%
%   See also: VL_GAUSSIAN(), VL_DGAUSSIAN(), VL_HELP().
[varargout{1:nargout}] = vl_ddgaussian(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function S = HALMSinit(w0,mu,fb,ff,delay)

% HALMSinit     Initialize Parameter Structure of the LMS Algorithm
%               for Hearing Aids. Refer to Fig. 1.9
%
% Arguments:
% w0            Coefficients of FIR filter at start (@n=1)
% mu            Step size for LMS algorithm
% ff            Feedforward path
% fb            Feedback path
% delay         Delay to adaptive filter
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

% Assign structure fields

S.coeffs     = w0(:);       % Weight (column) vector of FIR filter
S.step       = mu;          % Step size of the LMS algorithm
S.iter       = 0;           % Iteration count
S.ff         = ff;          % Feedforward path
S.fb         = fb;          % Feedback path
S.delay      = delay;       % Delay
S.AdaptStart = length(w0);  % Running effect of adaptive filter

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function quadrule_test12 ( )

%*****************************************************************************80
%
%% TEST12 tests FEJER2_INTEGRATE_FAST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2007
%
%  Author:
%
%    John Burkardt
%
  order_max = 16;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST12\n' );
  fprintf ( 1, '  FEJER2_INTEGRATE_FAST applies a Fejer Type 2 rule;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The integration interval is [-1,1].\n' );
  fprintf ( 1, '  Quadrature order will vary.\n' );
  fprintf ( 1, '  Integrand will vary.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   Order          F1              F2              F3\n' );
  fprintf ( 1, '\n' );

  for order = 1 : order_max

    result(1) = fejer2_integrate_fast ( @f1, order );
    result(2) = fejer2_integrate_fast ( @f2, order );
    result(3) = fejer2_integrate_fast ( @f3, order );

    fprintf ( 1, '  %6d  %14f  %14f  %14f\n', order, result(1:3) );

  end

  fprintf ( 1, '\n' );

  result(1) = 46.0 * sinh ( 1.0 ) / 25.0 - 2.0 * sin ( 1.0 );
  result(2) = 1.5822329637296729331;
  result(3) = ( sqrt ( 2.0 ) + 3.0 * sqrt ( 6.0 ) ) / 6.0;

  fprintf ( 1, '  Exact   %14f  %14f  %14f\n', result(1:3) );

  return
end


text:
function bvanofit(mycat, seg1, seg2)
    % bvanofit  Calculates Freq-Mag functions (b-value) for two time-segments
    %   finds best fit to the foreground for a modified background
    %   assuming a change in time of the following types:
    %   Mnew = Mold + d     , i.e. Simple magnitude shift
    %   Mnew = c*Mold + d   , i.e. Mag stretch plus shift
    %   Nnew = fac*Nold     , i.e. Rate change (N = number of events)
    %
    % bvanofit(catalog)
    %                                      R. Zuniga IGF-UNAM/GI-UAF  6/94
    % turned into function by Celso G Reyes 2017

    % used to work from newt2 and set newcat

    report_this_filefun();
    ZG=ZmapGlobal.Data;
    ms3 = 5;
    t1p=seg1(1);
    t2p=seg1(2);
    t3p=seg2(1);
    t4p=seg2(2);
    % This is the info window text
    %
    ttlStr='Comparing Seismicity rates ';
    hlpStr1map= ...
        ['                                                '
        ' To be Implemented                              '
        '                                                '];
    % Find out if figure already exists
    %
    bvfig=findobj('Type','Figure','-and','Name','Compare two rates');

    % Set up the Seismicity Map window Enviroment
    %
    if isempty(bvfig)
        bvfig= figure_w_normalized_uicontrolunits( ...
            'Name','Compare two rates',...
            'NumberTitle','off', ...
            'backingstore','on',...
            'Visible','on');%, ...
        %'Position',position_in_current_monitor(ZG.map_len(1), ZG.map_len(2)+200));

        %{
        uicontrol('Units','normal',...
            'Position',[.0 .93 .08 .06],'String','Print ',...
            'callback',@callbackfun_001)

        uicontrol('Units','normal',...
            'Position',[.0 .75 .08 .06],'String','Close ',...
            'callback',@callbackfun_002)

        uicontrol('Units','normal',...
            'Position',[.0 .85 .08 .06],'String','Info ',...
            'callback',@callbackfun_003)
        %}
        axis off

    end % if figure exits

    figure(bvfig);
    clf
    %delete(findobj(bvfig,'Type','axes'));
    format short;

    if isempty(mycat),
        mycat = ZG.primeCatalog;  % points to same thing
    end
    [minmag2, maxmag] = bounds(mycat.Magnitude);
    n = mycat.Count;
    tdiff = round(days(mycat.DateSpan()));

    % number of mag units
    %nmagu = (maxmag*10)+1;

    %backg_beN = [ ];
    %backg_abN = [ ];
    td12 = t2p - t1p;
    td34 = t4p - t3p;

    magsteps_desc = (maxmag:-0.1:minmag2);
    xt3edges = fliplr([magsteps_desc+.05 magsteps_desc(end)-.05]);

    %bg_seg = segmentAnalytics(mycat,'background',t1p,t2p); %instead of first segment
    %fg_seg = segmentAnalytics(mycat,'foreground',t3p,t4p); %instead of second segment

    %% first segment
    l = mycat.Date > t1p & mycat.Date < t2p ;
    backg =  mycat.subset(l);
    [bval,~] = histcounts(backg.Magnitude,xt3edges);%hist
    xt2=magsteps_desc; % based off of magnitude centers, not edges
    bval = bval / days(td12);                      % normalization
    assert(size(bval,1)<=1,'expecting a long row');
    bvalsum = cumsum(bval);                        % N for M <=
    bvalsum_bkw = cumsum(fliplr(bval));    % N for M >= (counted backwards) (was valsum3)
    %[cumux, xt] = hist(mycat.Date(l),t1p:days(ZG.bin_dur):t2p);
    [cumux, xt] = histcounts(mycat.Date(l),'BinWidth',ZG.bin_dur);
    mean1 = mean(cumux);
    var1 = cov(cumux);

    %% second segment
    l = mycat.Date > t3p & mycat.Date < t4p ;
    foreg = mycat.subset(l);
    [bval2,xt3a] = histcounts(foreg.Magnitude,xt3edges); %was histogram
    bval2 = bval2/days(td34);
    bvalsum2 = cumsum(bval2);
    bvalsum2_bkw=cumsum(fliplr(bval2)); % was bvalsum4
    % [cumux2, xt] = hist(mycat.Date(l),t3p:days(ZG.bin_dur):t4p);
    [cumux2, xt] = histcounts(mycat.Date(l),'BinWidth',ZG.bin_dur);
    mean2 = mean(cumux2);
    var2 = cov(cumux2);
    zscore = (mean1 - mean2)/(sqrt(var1/length(cumux)+var2/length(cumux2)));

    %change in percent
    R1 = backg.Count/days(t2p-t1p);
    R2 = foreg.Count/days(t4p-t3p);
    change = -((R1-R2)/R1)*100;

    %{
    backg_be = log10(bvalsum);
    backg_ab = log10(bvalsum_bkw);
    foreg_be = log10(bvalsum2);
    foreg_ab = log10(bvalsum2_bkw);
    %}

    % plot b-value plot
    %
    orient tall
    set(gcf,'PaperPosition',[2 1 5.5 7.5])
    p1 = subplot(5,1,[1 2]);
    %rect = [0.20,  0.7, 0.70, 0.25];           % plot Freq-Mag curves
    %axes('position',rect)
    figure(bvfig);

    d2char=@(d)[char(d,'yyyy-') num2str(days(d-dateshift(d,'start','year')),'%06.2f') ];

    segA_linespec = 'xb-';
    segA_dispname = [ d2char(t3p) ' - '  d2char(t4p)];

    segB_linespec = 'om-.';
    segB_dispname = [d2char(t1p) ' - ' d2char(t2p)]

    semilogy(p1,magsteps_desc,bvalsum2_bkw,segA_linespec,'MarkerSize',ms3,'DisplayName',segA_dispname);
    set(gca,'NextPlot','add')
    semilogy(p1,magsteps_desc,bvalsum_bkw,segB_linespec,'MarkerSize',ms3,'DisplayName',segB_dispname);
    legend show

    te1 = max([bvalsum  bvalsum2 bvalsum2_bkw bvalsum_bkw]);
    te1 = te1 - 0.2*te1;

    ylabel(p1,'Cum. rate/year','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold')
    str = { 'Event rate comparison for 2 time periods',['Change in %: ' num2str(change,6) ]};

    title(p1,str,'FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold')

    %  find b-values;
    set(p1,'box','on',...
        'SortMethod','childorder','TickDir','out','FontWeight',...
        'bold','FontSize',ZmapGlobal.Data.fontsz.s,'Linewidth',1.0)

    % Plot histogram
    %rect = [0.20,  0.40 0.70, 0.25];
    %axes('position',rect)
    p2=subplot(5,1,[3 4]);
    plot(p2,xt2,bval2,segA_linespec,'MarkerSize',ms3,'LineWidth',1.0)
    set(gca,'NextPlot','add')
    plot(p2,xt2,bval,segB_linespec,'MarkerSize',ms3,'LineWidth',1.0)
    disp([' Summation: ' num2str(sum(bval-bval2))])
    v = axis;
    xlabel(p2,'Magnitude ','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold')
    ylabel(p2,'rate/year','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','bold')
    set(p2,'box','on',...
        'SortMethod','childorder','TickDir','out','FontWeight',...
        'bold','FontSize',ZmapGlobal.Data.fontsz.s,'Linewidth',1.0)

    uic = uicontrol('Units','normal','Position',[.35 .05 .30 .06],'String','Magnitude Signature? ', 'callback',@callbackfun_004,'enable','off');

    watchoff

    % Plot he b-value comparison
    ZG.hold_state=false;
    bdiff(backg)
    ZG.hold_state=true;
    bdiff(foreg)

    function callbackfun_004(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        delete(uic);
        synsig3;
    end

    function seg = segmentAnalytics(mycat, name, t1, t2)
        % alternate way of dong this for the functions
        seg=struct('name',name,'start',t1,'end',t2);
        seg.duration= seg.end - seg.start;
        index = mycat.Date > seg.start & mycat.Date < seg.end;
        seg.catalog = mycat.subset(index);
        seg.bvals = histcounts(seg.catalog.Magnitude, xt3edges);
        seg.scaledbvals = seg.bvals / days(seg.duration);
        seg.bvalcum = cumsum(seg.scaledbvals);
        seg.bvalcum_bkw = cumsum(fliplr(seg.scaledbvals));
        [seg.datecounts, seg.date_edges] = histcounts(seg.catalog.Date,'BinWidth',ZG.bin_dur);
        seg.meanEventsTime = mean(seg.datecounts);
        seg.covEventsTime = cov(seg.datecounts);
    end

end

author: CelsoReyes


text:
% Assigns the given value to the named property of an object or structure.
% This function can deal with nested properties.
%
% Input arguments:
% obj:
%    the structure, handle or value object the value should be assigned to
% name:
%    a property name with dot (.) separating property names at
%    different hierarchy levels
% value:
%    the value to assign to the property at the deepest hierarchy
%    level
%
% Output arguments:
% obj:
%    the updated object or structure, optional for handle objects
%
% Example:
%    obj = struct('surface', struct('nested', 10));
%    obj = nestedassign(obj, 'surface.nested', 23);
%    disp(obj.surface.nested);  % prints 23
%
% See also: nestedfetch

% Copyright 2010 Levente Hunyadi
function obj = nestedassign(obj, name, value)
    if ~iscell(name)
        nameparts = strsplit(name, '.');
    else
        nameparts = name;
    end
    obj = nestedassign_recurse(obj, nameparts, value);
end

function obj = nestedassign_recurse(obj, name, value)
% Assigns the given value to the named property of an object.
%
% Input arguments:
% obj:
%    the handle or value object the value should be assigned to
% name:
%    a cell array of the composite property name
% value:
%    the value to assign to the property at the deepest hierarchy
%    level
    if numel(name) > 1
        obj.(name{1}) = nestedassign_recurse(obj.(name{1}), name(2:end), value);
    else
        obj.(name{1}) = value;
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
%   TTeMPS Toolbox.
%   Michael Steinlechner, 2013-2016
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

function [eta, B1,B3] = precond_laplace_overlapJacobi( L, xi, xL, xR, G, B1, B3 )
% L is a cell of operators

r = xi.rank;
n = xi.size;
d = xi.order;

% If B1 and B3 are not given as arguments, we need to precalculate them
if nargin < 7
%     % if applying L is expensive (not just tridiag), one can store all
%     % applications with xL and compute the ones for xR with G.
%     % You need to first store LUl
%     LUl = cell(d,1);
%     for idx = 1:d
%         LUl{idx} = tensorprod_ttemps( xL.U{idx}, L{idx}, 2 );
%     end
%     % and then change to LUr in the loop for B3 below
%     %         if idx+1==d
%     %              LUr = tensorprod_ttemps( LUl{idx+1}, G{idx}, 1, true);
%     %         else
%     %             LUr = tensorprod_ttemps( tensorprod_ttemps( LUl{idx+1}, G{idx+1}', 3), G{idx}, 1, true);
%     %         end

    B1 = cell(d,1);
    B1{1} = 0;
    for idx = 2:d
        LUl = tensorprod_ttemps( xL.U{idx-1}, L{idx-1}, 2 );
        if idx>2
            TT = tensorprod_ttemps( xL.U{idx-1}, B1{idx-1}, 1 );
        else
            TT = 0;
        end
        B1{idx} = unfold(xL.U{idx-1},'left')'*unfold(TT + LUl,'left');
    end

    B3 = cell(d,1);
    for idx = d-1:-1:1
        LUr = tensorprod_ttemps( xR.U{idx+1}, L{idx+1}, 2 );
        if idx<d-1
            TT = tensorprod_ttemps( xR.U{idx+1}, B3{idx+1}, 3 );
        else
            TT = 0;
        end
        B3{idx} = unfold(xR.U{idx+1},'right')*unfold(TT + LUr,'right')';
    end
    B3{d} = 0;
end

eta = xi;
xi = tangent_to_TTeMPS( xi );

% % 1. STEP: Project right hand side
% below is hard-coded version of
% for ii=1:d
%     eta_partial_ii = TTeMPS_partial_project_overlap( xL, xR, xi, ii);
%     Y{ii} = eta_partial_ii.dU{ii};
% end

% TODO, it seems that the left and right cell arrays consist of a lot of
% identities and zeros.
Y = cell(1,d);
% precompute inner products
left = innerprod( xL, xi, 'LR', d-1, true );
right = innerprod( xR, xi, 'RL', 2, true );

% contract to first core
Y{1} = tensorprod_ttemps( xi.U{1}, right{2}, 3 );
% contract to first core
for idx = 2:d-1
    res = tensorprod_ttemps( xi.U{idx}, left{idx-1}, 1 );
    Y{idx} = tensorprod_ttemps( res, right{idx+1}, 3 );
end
% contract to last core
Y{d} = tensorprod_ttemps( xi.U{d}, left{d-1}, 1 );

% 2. STEP: Solve ALS systems:

% B1 and B3 were precalculated before
for idx = 1:d
    rl = r(idx);
    rr = r(idx+1);

    B2 = L{idx};

    % Solve via the diagonalization trick
    [V1,E1] = eig(B1{idx}); [V3,E3] = eig(B3{idx});
    V = kron(V3,V1);
    EE = diag(E1)*ones(1,rr) + ones(rl,1)*diag(E3)'; E = EE(:);

    rhs = matricize( Y{idx}, 2 ) * V;
    Z = zeros(size(rhs));
    for i=1:length(E)
        Z(:,i) = (B2 + E(i)*speye(n(idx))) \ rhs(:,i);
    end
    eta.dU{idx} = tensorize( Z*V', 2, [rl, n(idx), rr] );
end

eta = TTeMPS_tangent_orth( xL, xR, eta );   % todo? Can we improve efficiency since eta is not a generic TTeMPS but shares the same x.U as xL and xR

end

author: NicolasBoumal

repo: manopt


text:
function [xl,yl,zl]=axesLimits(V)
%% calculate axes limits for plotting
% V can be a nX3 matrix of 3D vertices or a cell array which contains nX3
% matrices representing different time frames.
%%
if ~iscell(V)
    Vtemp=V;
    V=cell(1);
    V{1}=Vtemp;
end

% xl=[0 0]; yl=[0 0]; zl=[0 0];
xl=[min(V{1}(:,1)) max(V{1}(:,1))]; yl=[min(V{1}(:,2)) max(V{1}(:,2))]; zl=[min(V{1}(:,3)) max(V{1}(:,3))];
for it=1:numel(V)
    xl(1)=min([min(V{it}(:,1)) xl(1)]);
    xl(2)=max([max(V{it}(:,1)) xl(2)]);
    yl(1)=min([min(V{it}(:,2)) yl(1)]);
    yl(2)=max([max(V{it}(:,2)) yl(2)]);
    zl(1)=min([min(V{it}(:,3)) zl(1)]);
    zl(2)=max([max(V{it}(:,3)) zl(2)]);
end

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function r = VBA_spm_gamrnd(a,b,varargin)
% Random arrays from gamma distribution - a compiled routine
% FORMAT r = spm_gamrnd(a,b,m,n,...)
%
% a        - shape parameter
% b        - scale parameter
% m,n,...  - dimensions of the output array [optional]
%
% r        - array of random numbers chosen from the gamma distribution
%__________________________________________________________________________
%
% Reference
%
% George Marsaglia and Wai Wan Tsang, "A Simple Method for Generating Gamma
% Variables": ACM Transactions on Mathematical Software, Vol. 26, No. 3,
% September 2000, Pages 363-372
% http://portal.acm.org/citation.cfm?id=358414
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_gamrnd.m 3251 2009-07-06 17:29:44Z guillaume $

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

author: MBB-team

repo: VBA-toolbox


text:
function Gxy = ffmGreenKernel(X,Y,green,k)
%+========================================================================+
%|                                                                        |
%|         OPENFFM - LIBRARY FOR FAST AND FREE MEMORY CONVOLUTION         |
%|           openFfm is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2019.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : ffmGreenKernel.m                              |
%|    #    |   VERSION    : 0.6                                           |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2019                                    |
%| ( === ) |   SYNOPSIS   : Green kernel computation using string         |
%|  `---'  |                definition                                    |
%+========================================================================+

% Security
if (size(X,2) ~= 3) || (size(Y,2) ~= 3)
    error('ffmGreenKernel.m : unavailable case')
end
if isempty(k)
    k = 0;
end

% Distances between particles
Rxy = sqrt( ...
    (X(:,1) - Y(:,1)).^2 + ...
    (X(:,2) - Y(:,2)).^2 + ...
    (X(:,3) - Y(:,3)).^2 );

% For empty wave-number
if isempty(k)
    k = 0;
end

% Green kernel definition
if strcmp(green,'[1/r]')
    Gxy = 1./Rxy;

elseif strcmp(green,'[exp(ikr)/r]')
    Gxy = exp(1i*k*Rxy)./Rxy;

elseif strcmp(green(1:end-1),'gradx[1/r]')
    j = str2double(green(end));
    Gxy = - (X(:,j)-Y(:,j)) ./ (Rxy.^3);

elseif strcmp(green(1:end-1),'grady[1/r]')
    j = str2double(green(end));
    Gxy = (X(:,j)-Y(:,j)) ./ (Rxy.^3);

elseif strcmp(green(1:end-1),'gradx[exp(ikr)/r]')
    j = str2double(green(end));
    Gxy = (1i*k - 1./Rxy) .* exp(1i*k.*Rxy) .* ...
        (X(:,j)-Y(:,j)) ./ (Rxy.^2);

elseif strcmp(green(1:end-1),'grady[exp(ikr)/r]')
    j = str2double(green(end));
    Gxy = - (1i*k - 1./Rxy) .* exp(1i*k.*Rxy) .* ...
        (X(:,j)-Y(:,j)) ./ (Rxy.^2);

elseif strcmp(green(1:end-2),'[ij/r+rirj/r^3]')
    i = str2double(green(end-1));
    j = str2double(green(end));
    Gxy = (i==j)./Rxy + (X(:,i)-Y(:,i)).*(X(:,j)-Y(:,j))./(Rxy.^3);

elseif strcmp(green(1:end-3),'[rirjrk/r^5]')
    i = str2double(green(end-2));
    j = str2double(green(end-1));
    k = str2double(green(end));
    Gxy = (X(:,i)-Y(:,i)).*(X(:,j)-Y(:,j)).*(X(:,k)-Y(:,k))./(Rxy.^5);

else
    error('Error in ffmGreenKernel.m : unknown green kernel')
end

% Singularity
if strcmp(green,'[exp(ikr)/r]')
    Gxy(Rxy<1e-8) = 0 + 1i*k;
else
    Gxy(Rxy<1e-8) = 0;
end
end


text:
function computeAllModelChoice_batchVASARI_LGG(pathExperiments,maxOrder,nBoot,imbalance,nBatch,matlabPATH,seed)
% -------------------------------------------------------------------------
% function computeAllModelChoice_batchVASARI_LGG(pathExperiments,maxOrder,nBoot,imbalance,nBatch,matlabPATH,seed)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes feature set selection for all feature set types
% and for all experiments with different degrees of freedom. See ref. [1]
% for more details.
% -------------------------------------------------------------------------
% REFERENCES:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathExperiments: Full path to where all experiments need to be
%                     performed.
%                     --> Ex: /myProject/WORKSPACE/VASARI
% 2. maxOrder: Integer specifying the maximal model order to construct.
%              --> Ex: 10
% 3. nBoot: Number of bootstrap samples to use.
%           --> Ex: 100
% 4. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (formal reference to come).
%               --> Ex: 'IALR'
% 5. nBatch: Number of parallel batch.
%            --> Ex: 8
% 6. matlabPATH: Full path to the MATLAB executable on the system.
%                --> 'matlab' if a symbolic link to the matlab executable
%                     was previously created.
% 7. seed: Numerical number to use as seed for bootstrapping experiment
%          --> Ex: 54288
% -------------------------------------------------------------------------
% OUTPUTS: Mutivariable models are saved in a folder named 'MODELS' in the
% corresponding folder of 'pathExperiments'.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;

% INITIALIZATON
time = 60; % Number of seconds to wait before checking if parallel computations are done
cd(pathExperiments), load('outcomes')
pathSet = fullfile(pwd,'FSET'); mkdir('MODELS'), cd('MODELS'), pathModels = pwd;
mkdir('batchLog_Models'), cd('batchLog_Models'), pathBatch = pwd;
setNames = {'VASARI'};
[param] = batchExperiments(setNames,outcomes,nBatch); nBatch = length(param);

% PRODUCE BATCH COMPUTATIONS
save('workspace','pathSet','pathModels','outcomes','param','maxOrder','nBoot','imbalance','seed'), pause(5);
for i = 1:nBatch
    nameScript = ['batch',num2str(i),'_script.m'];
    fid = fopen(nameScript,'w');
    fprintf(fid,'load(''workspace'')\n');
    for j = 1:numel(param{i})
        fprintf(fid,['computeAllModelChoice_VASARI_LGG(pathSet,pathModels,outcomes,param{',num2str(i),'}{',num2str(j),'},maxOrder,nBoot,imbalance,seed)\n']);
    end
    fprintf(fid,['system(''touch batch',num2str(i),'_end'');\n']);
    fprintf(fid,'clear all');
    fclose(fid);
    system([matlabPATH,' -nojvm -nodisplay -nodesktop -nosplash < ',nameScript,' >& ',nameScript(1:end-1),'log &']);
end

% WAITING LOOP
waitBatch(pathBatch,time,nBatch)
delete('workspace.mat')

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [satp] = satposin(t,eph)
%  SATPOSIN   Calculation of X,Y,Z coordinates in an INERTIAL
%	      reference frame at time t for given ephemeris eph

% Written by Kai Borre
% November 15, 1996

   GM = 3.986008e14;		% earth's universal gravitational parameter
				% m^3/s^2

   %  Units are either seconds, meters, or radians
   %  Assigning the local variables to eph
   svprn   =   eph(1);
   af2	   =   eph(2);
   M0	   =   eph(3);
   roota   =   eph(4);
   deltan  =   eph(5);
   ecc	   =   eph(6);
   omega   =   eph(7);
   cuc	   =   eph(8);
   cus	   =   eph(9);
   crc	   =  eph(10);
   crs	   =  eph(11);
   i0	   =  eph(12);
   idot    =  eph(13);
   cic	   =  eph(14);
   cis	   =  eph(15);
   Omega0  =  eph(16);
   Omegadot=  eph(17);
   toe	   =  eph(18);
   af0	   =  eph(19);
   af1	   =  eph(20);
   t0c	   =  eph(21);

   % Procedure for coordinate calculation
   A = roota*roota;
   t = check_t(t-t0c);
   n0 = sqrt(GM/A^3);
   n = n0+deltan;
   M = M0+n*t;
   M = rem(M,2*pi);
   E = M;
   E_old = E+.001;
   while abs(E-E_old) >= 1.e-12
      E_old = E;
      E = M+ecc*sin(E);
   end
   v = atan2(sqrt(1-ecc^2)*sin(E), cos(E)-ecc);
   phi = v+omega;
   phi = rem(phi,2*pi);
   u = phi		      + cuc*cos(2*phi)+cus*sin(2*phi);
   r = A*(1-ecc*cos(E))       + crc*cos(2*phi)+crs*sin(2*phi);
   i = i0+idot*t	      + cic*cos(2*phi)+cis*sin(2*phi);
   Omega = Omega0+Omegadot*t;
   x1 = cos(u)*r;
   y1 = sin(u)*r;
   satp(1,1) = x1*cos(Omega)-y1*cos(i)*sin(Omega);
   satp(2,1) = x1*sin(Omega)+y1*cos(i)*cos(Omega);
   satp(3,1) = y1*sin(i);

%%%%%%%%% end satposin.m %%%%%%%%%

author: yandld


text:
function F = flipud(F)
%FLIPUD   Flip the rows of a chebmatrix.

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

F.blocks = flipud(F.blocks);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [positionDecodingMaxCorr] = bz_positionDecodingMaxCorr(varargin)
% USAGE
%   [] = bz_positionDecodingMaxCorr()
%
% INPUTS
%
%   spikes
%
%   behavior
%
%   lfp
%
%   smoothingRange
%
%   plotting
%
%   saveMat
%
% OUTPUTS
%
%   positionDecodingMaxCorr     cellinfo struct with the following fields
%                      .UID
%                      .region
%                      .sessionName
%                      .results{1:nCells}
%                      .tau
%                      .condition
%
%
%
%
% SEE
%
% written by david tingley, 2017

%% parse inputs
p = inputParser();
addRequired(p,'spikes',@isstruct);
addRequired(p,'behavior',@isstruct);
addRequired(p,'lfp',@isstruct);
addRequired(p,'smoothingRange',@isvector);

addParameter(p,'plotting',@islogical);
addParameter(p,'saveMat',@islogical);
parse(p,varargin{:})

spikes = p.Results.spikes;
behavior = p.Results.behavior;
lfp = p.Results.lfp;
smoothingRange = p.Results.smoothingRange;
plotting = p.Results.plotting;
saveMat = p.Results.saveMat;

%% set up data format and popinfo struct that will be returned

positionDecodingMaxCorr.UID = spikes.UID;
positionDecodingMaxCorr.region = spikes.region;
positionDecodingMaxCorr.sessionName = spikes.sessionName; % session name

conditions = unique(behavior.events.trialConditions);
nCells = length(spikes.times);
positionSamplingRate = behavior.samplingRate;

% find a better way to get spike phase relationship...
[firingMaps] = bz_firingMap1D(spikes,behavior,5);
[phaseMaps] = bz_phaseMap1D(spikes,behavior,lfp,5);
[b a] = butter(3,[4/625 12/625],'bandpass');
theta_phases = angle(hilbert(FiltFiltM(b,a,double(lfp.data))));
% iterate through conditions and compile spike trains and spike-phase
% trains
for cond = conditions
    trials = find(behavior.events.trialConditions==cond);
    intervals = behavior.events.trialIntervals(trials,:);
    for t = 1:length(trials)
        trial = trials(t);
        spk_trains{cond}{t} = zeros(nCells,ceil((intervals(t,2)-intervals(t,1))*1000)); % assumes intervals are in seconds, rounds to nearest millisecond
        phase_trains{cond}{t} = zeros(nCells,ceil((intervals(t,2)-intervals(t,1))*1000));

        for cell = 1:nCells
            if ~isempty(phaseMaps.phaseMaps{cond}{cell})
                f = find(phaseMaps.phaseMaps{cond}{cell}(:,2)==t);
                if ~isempty(f)
                for s=1:length(f)
                    phase_trains{cond}{t}(cell,ceil(phaseMaps.phaseMaps{cond}{cell}(f(s),5)*1000)) = ...
                        phaseMaps.phaseMaps{cond}{cell}(f(s),end);
                end
                end
            end
            sp = find(InIntervals(spikes.times{cell},intervals(t,:)));
            spk_trains{cond}{t}(cell,ceil((spikes.times{cell}(sp)-intervals(t,1))*1000+.00001))=1;
        end

%         position{cond}{t} = interp1(1:length(behavior.events.trials{trial}.x)...
%             ,behavior.events.trials{trial}.mapping,1:positionSamplingRate/1000:length(...
%             behavior.events.trials{trial}.x));
        nBins = size(spk_trains{cond}{t},2);
        nPos = length(behavior.events.trials{trial}.x);
        position{cond}{t} = interp1(1:length(behavior.events.trials{trial}.x)...
            ,behavior.events.trials{trial}.mapping,1:(nPos-1)/nBins:length(...
            behavior.events.trials{trial}.x)); % the -1 gaurantees the length to be longer than the above spk/phase trains
        position{cond}{t} = position{cond}{t}(1:length(spk_trains{cond}{t}));
        take = FindInInterval(lfp.timestamps,intervals(t,:));
        theta{cond}{t} = makeLength(theta_phases(take(1):take(2)),length(position{cond}{t}));
    end
end

%% set up MaxCorr decoder
positionDecodingMaxCorr.results = table;

disp('running models...')
warning off % max cl model returns warning with NaNs and Inf values..
for cond = conditions
    r = randperm(length(phase_trains{cond}));
    for iter = 1:5%length(r)
        r = circshift(r,1);
    for wind = smoothingRange
        for cell = 1:nCells
            phase_trains_smooth=[];
            rates_trains_smooth = [];
            position_train = [];
            theta_train = [];
            for t = 1:length(phase_trains{cond})-1
                phase_trains_smooth=[phase_trains_smooth; circ_smoothTS(phase_trains{cond}{r(t)}(cell,:),wind,'method','mean','exclude',0)];
                rates_trains_smooth = [rates_trains_smooth; smooth(spk_trains{cond}{r(t)}(cell,:),wind)*wind];
                position_train = [position_train; position{cond}{r(t)}'];
                theta_train = [theta_train, theta{cond}{r(t)}];
            end
            phase_trains_smooth_train(cell,:) = phase_trains_smooth;
            rates_trains_smooth_train(cell,:) = rates_trains_smooth;
            phase_trains_smooth_test(cell,:)=[...
                circ_smoothTS(phase_trains{cond}{r(end)}(cell,:),wind,'method','mean','exclude',0)];
            rates_trains_smooth_test(cell,:) = [...
                smooth(spk_trains{cond}{r(end)}(cell,:),wind)*wind];

        end
        position_test = position{cond}{r(end)};
        theta_test = theta{cond}{r(end)};

        %% rate coding model
        cl = max_correlation_coefficient_CL;
        cl = train(cl,[rates_trains_smooth_train; theta_train],round(position_train));

        for ts = 1:length(position_test)
           yfit_rate(ts) = test(cl,[round(rates_trains_smooth_test(:,ts));theta_test(ts)]);
        end
        struct.mse_rate = mean((yfit_rate-position_test).^2);
        % chance rate
        rr = randperm(length(theta_train));
        rrr = randperm(length(theta_test));
        cl = max_correlation_coefficient_CL;
        cl = train(cl,[rates_trains_smooth_train(:,rr); theta_train],round(position_train));

        for ts = 1:length(position_test)
           yfit_chance_rate(ts) = test(cl,[round(rates_trains_smooth_test(:,rrr(ts)));theta_test(ts)]);
        end
        struct.mse_chance_rate = mean((yfit_chance_rate-position_test).^2);

        %% phase coding model

        % discretize phase_trains here...
        phase_trains_smooth_train_cos = cos(phase_trains_smooth_train);
        phase_trains_smooth_train_sin = sin(phase_trains_smooth_train);
        phase_trains_smooth_test_cos = cos(phase_trains_smooth_test);
        phase_trains_smooth_test_sin = sin(phase_trains_smooth_test);

        phase_trains_smooth_train(phase_trains_smooth_train==0)=nan;
        phase_trains_smooth_test(phase_trains_smooth_test==0)=nan;
        phase_trains_smooth_train = discretize(phase_trains_smooth_train,-pi:.1:pi);
        phase_trains_smooth_test = discretize(phase_trains_smooth_test,-pi:.1:pi);
        phase_trains_smooth_train(isnan(phase_trains_smooth_train))=0;
        phase_trains_smooth_test(isnan(phase_trains_smooth_test))=0;

        phase_trains_smooth_train_cos(phase_trains_smooth_train_cos==0)=nan;
        phase_trains_smooth_test_cos(phase_trains_smooth_test_cos==0)=nan;
        phase_trains_smooth_train_cos = discretize(phase_trains_smooth_train_cos,-1:.1:1);
        phase_trains_smooth_test_cos = discretize(phase_trains_smooth_test_cos,-1:.1:1);
        phase_trains_smooth_train_cos(isnan(phase_trains_smooth_train_cos))=0;
        phase_trains_smooth_test_cos(isnan(phase_trains_smooth_test_cos))=0;

        phase_trains_smooth_train_sin(phase_trains_smooth_train_sin==0)=nan;
        phase_trains_smooth_test_sin(phase_trains_smooth_test_sin==0)=nan;
        phase_trains_smooth_train_sin = discretize(phase_trains_smooth_train_sin,-1:.1:1);
        phase_trains_smooth_test_sin = discretize(phase_trains_smooth_test_sin,-1:.1:1);
        phase_trains_smooth_train_sin(isnan(phase_trains_smooth_train_sin))=0;
        phase_trains_smooth_test_sin(isnan(phase_trains_smooth_test_sin))=0;

        % non-transformed circular decoding
        cl = max_correlation_coefficient_CL;
        cl = train(cl,[phase_trains_smooth_train; theta_train],round(position_train));

        for ts = 1:length(position_test)
           yfit_circ(ts) = test(cl,[phase_trains_smooth_test(:,ts);theta_test(ts)]);
        end
        struct.mse_phase = mean((yfit_circ-position_test).^2);

        % cos and sin transformed circular decoding
        cl = max_correlation_coefficient_CL;
        cl = train(cl,[phase_trains_smooth_train_cos; theta_train],round(position_train));

        for ts = 1:length(position_test)
           yfit_cos(ts) = test(cl,[phase_trains_smooth_test_cos(:,ts);theta_test(ts)]);
        end
        struct.mse_phase_cos = mean((yfit_cos-position_test).^2);

        cl = max_correlation_coefficient_CL;
        cl = train(cl,[phase_trains_smooth_train_sin; theta_train],round(position_train));

        for ts = 1:length(position_test)
           yfit_sin(ts) = test(cl,[phase_trains_smooth_test_sin(:,ts);theta_test(ts)]);
        end
        struct.mse_phase_sin = mean((yfit_sin-position_test).^2);

        % all phase models in one...
        cl = max_correlation_coefficient_CL;
        all_phase_train = [phase_trains_smooth_train_cos;phase_trains_smooth_train_sin];
        all_phase_test =  [phase_trains_smooth_test_cos;phase_trains_smooth_test_sin];
        cl = train(cl,[all_phase_train; theta_train],round(position_train));

        for ts = 1:length(position_test)
           yfit_circ_all(ts) = test(cl,[all_phase_test(:,ts);theta_test(ts)]);
        end
        struct.mse_phase_all = mean((yfit_circ_all-position_test).^2);

        % chance phase
        cl = max_correlation_coefficient_CL;
        all_phase_train = [phase_trains_smooth_train_cos;phase_trains_smooth_train_sin];
        all_phase_test =  [phase_trains_smooth_test_cos;phase_trains_smooth_test_sin];
        cl = train(cl,[all_phase_train(:,rr); theta_train],round(position_train));

        for ts = 1:length(position_test)
           yfit_chance(ts) = test(cl,[all_phase_test(:,rrr(ts));theta_test(ts)]);
        end
        struct.mse_chance = mean((yfit_chance-position_test).^2);

        %% put data into struct/table
        struct.tau = wind;
        struct.condition = cond;
        struct.iter = iter;
        struct.trialOrder = r;
        positionDecodingMaxCorr.results = [positionDecodingMaxCorr.results;struct2table(struct)];
        if plotting
            clf
            subplot(2,2,1)
            t_rate = varfun(@mean,positionDecodingMaxCorr.results,'InputVariables','mse_rate',...
            'GroupingVariables',{'tau','condition'});
            t_phase = varfun(@mean,positionDecodingMaxCorr.results,'InputVariables','mse_phase',...
            'GroupingVariables',{'tau','condition'});
            t_phase_cos = varfun(@mean,positionDecodingMaxCorr.results,'InputVariables','mse_phase_cos',...
            'GroupingVariables',{'tau','condition'});
            t_phase_sin = varfun(@mean,positionDecodingMaxCorr.results,'InputVariables','mse_phase_sin',...
            'GroupingVariables',{'tau','condition'});
            t_chance = varfun(@mean,positionDecodingMaxCorr.results,'InputVariables','mse_chance',...
            'GroupingVariables',{'tau','condition'});
            t_chance_s = varfun(@std,positionDecodingMaxCorr.results,'InputVariables','mse_chance',...
            'GroupingVariables',{'tau','condition'});
            t_chance_rate = varfun(@mean,positionDecodingMaxCorr.results,'InputVariables','mse_chance_rate',...
            'GroupingVariables',{'tau','condition'});
            t_chance_s_rate = varfun(@std,positionDecodingMaxCorr.results,'InputVariables','mse_chance_rate',...
            'GroupingVariables',{'tau','condition'});
            t_rate_s = varfun(@std,positionDecodingMaxCorr.results,'InputVariables','mse_rate',...
            'GroupingVariables',{'tau','condition'});
            t_phase_s = varfun(@std,positionDecodingMaxCorr.results,'InputVariables','mse_phase',...
            'GroupingVariables',{'tau','condition'});
            t_phase_cos_s = varfun(@std,positionDecodingMaxCorr.results,'InputVariables','mse_phase_cos',...
            'GroupingVariables',{'tau','condition'});
            t_phase_sin_s = varfun(@std,positionDecodingMaxCorr.results,'InputVariables','mse_phase_sin',...
            'GroupingVariables',{'tau','condition'});
            tab = join(join(join(join(t_rate,t_phase),t_chance),t_phase_cos),t_phase_sin);
            tab_s = join(join(join(join(t_rate_s,t_phase_s),t_chance_s),t_phase_sin_s),t_phase_cos_s);
            rows = find(tab.condition==cond);

            title('MaxCorr decoding of pos, r-rate, g-phase')
%             iterRows = positionDecodingMaxCorr.iter == iter;
%             plot(positionDecodingMaxCorr.results.tau,...
%                 positionDecodingMaxCorr.results.mse_rate,'r')
%             hold on
%             plot(positionDecodingMaxCorr.results.tau,...
%                 positionDecodingMaxCorr.results.mse_phase_all,'g')
%             hold off
%             boundedline(tab.tau,t_chance_rate.mean_mse_chance_rate(rows),t_chance_s_rate.std_mse_chance_rate(rows),'k')
            boundedline(tab.tau,tab.mean_mse_chance(rows),tab_s.std_mse_chance(rows),'k')
            boundedline(tab.tau,tab.mean_mse_phase_cos(rows),tab_s.std_mse_phase_cos(rows),'g')
            boundedline(tab.tau,tab.mean_mse_rate(rows),tab_s.std_mse_rate(rows),'r')
%             set(gca,'xscale','log')
            subplot(2,2,2)
            boundedline(tab.tau,tab.mean_mse_phase(rows),tab_s.std_mse_phase(rows),'.g')
            boundedline(tab.tau,tab.mean_mse_phase_sin(rows),tab_s.std_mse_phase_sin(rows),'--g')
            boundedline(tab.tau,tab.mean_mse_phase_cos(rows),tab_s.std_mse_phase_cos(rows),'g')
%             subplot(2,2,3)
%             plot(positionDecodingMaxCorr.results.tau,...
%                 positionDecodingMaxCorr.results.mse_phase_sin,'g')
%             subplot(2,2,4)
%             plot(positionDecodingMaxCorr.results.tau,...
%                 positionDecodingMaxCorr.results.mse_phase,'g')
            pause(.001)
        end
        clear *train *test yfit*
        disp(['finished with window: ' num2str(wind) ' out of ' num2str(smoothingRange(end)) ' total'])
    end
    end
    disp(['finished with condition: ' num2str(cond) ' out of ' num2str(length(conditions)) ' total'])
    positionDecodingMaxCorr.dateRun = date;
    if saveMat
       save([spikes.sessionName '.positionDecodingMaxCorr.popinfo.mat'],'positionDecodingMaxCorr')
    end
end


text:
function exp_test ( )

%*****************************************************************************80
%
%% EXP_TEST tests R4_EXP and R8_EXP.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXP_TEST:\n' );
  fprintf ( 1, '  Test EXP_VALUES, R4_EXP, R8_EXP\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X          EXP(X)\n' );
  fprintf ( 1, '                     R4_EXP(X)        Diff\n' );
  fprintf ( 1, '                     R8_EXP(X)        Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = exp_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_exp ( single ( x ) );
    fx3 = r8_exp ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Cowell et al p72
G = zeros(10);
G(1,2)=1;
G(2,3)=1;
G(3,7)=1;
G(4,[5 8])=1;
G(5,6)=1;
G(6,7)=1;
G(7,[9 10])=1;
G(8,9)=1;

dsep(1, 4, [5 7], G)
dsep(1, 4, [7], G)
dsep(1, 4, [10 5], G)

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [check,box] = checkfp(ha,plr,box,ipx,ipy,fpx,fpy,dblbox,blank)
    if plr==1
        otherplr=2;
    elseif plr==2
        otherplr=1;
    end
    %at initial position there should be own piece
    %and at final there should be blank space
    if box(fpx,fpy)~=0
        check=0;
    %     back step is not allowed till the palyer has choose doubled piece
    elseif fpx<=ipx && plr==1 && dblbox(ipx,ipy)==0
        check=0;
    elseif fpx>=ipx && plr==2 && dblbox(ipx,ipy)==0
        check=0;
    %only diagonal movement is allowed
    elseif (fpx==ipx+1 && fpy==ipy+1) || (fpx==ipx-1 && fpy==ipy+1) || ...
            (fpx==ipx+1 && fpy==ipy-1) || (fpx==ipx-1 && fpy==ipy-1)
        check=1;
    %if approach to kill a single piece and opponent piece is
    %there....its a good move .. so let the player play it..killing will
    %include that specific box=0 and string=blank
    elseif (fpx==ipx+2 && fpy==ipy+2) && ( box(ipx+1,ipy+1)==otherplr )
            check=1;
            box(ipx+1,ipy+1)=0;
            xy=(ipx+1)*10+(ipy+1);
            set(ha(xy),'CData',blank);
    elseif (fpx==ipx-2 && fpy==ipy-2) && ( box(ipx-1,ipy-1)==otherplr )
            check=1;
            box(ipx-1,ipy-1)=0;
            xy=(ipx-1)*10+(ipy-1);
            set(ha(xy),'CData',blank);
    elseif (fpx==ipx+2 && fpy==ipy-2) && ( box(ipx+1,ipy-1)==otherplr )
            check=1;
            box(ipx+1,ipy-1)=0;
            xy=(ipx+1)*10+(ipy-1);
            set(ha(xy),'CData',blank);
    elseif (fpx==ipx-2 && fpy==ipy+2) && ( box(ipx-1,ipy+1)==otherplr )
            check=1;
            box(ipx-1,ipy+1)=0;
            xy=(ipx-1)*10+(ipy+1);
            set(ha(xy),'CData',blank);
    %if approach to kill two pieces together lets check
    elseif (fpx==ipx+4 && fpy==ipy+4) && ( box(ipx+1,ipy+1)==otherplr) && (box(ipx+2,ipy+2)==0) && (box(ipx+3,ipy+3)==otherplr)
        check=1;
        box(ipx+1,ipy+1)=0;
        box(ipx+3,ipy+3)=0;
        xy=(ipx+1)*10+(ipy+1);
        set(ha(xy),'CData',blank);
        xy=(ipx+3)*10+(ipy+3);
        set(ha(xy),'CData',blank);
    elseif (fpx==ipx+4 && fpy==ipy-4) && ( box(ipx+1,ipy-1)==otherplr) && (box(ipx+2,ipy-2)==0) && (box(ipx+3,ipy-3)==otherplr)
        check=1;
        box(ipx+1,ipy-1)=0;
        box(ipx+3,ipy-3)=0;
        xy=(ipx+1)*10+(ipy-1);
        set(ha(xy),'CData',blank);
        xy=(ipx+3)*10+(ipy-3);
        set(ha(xy),'CData',blank);
    elseif (fpx==ipx-4 && fpy==ipy+4) && ( box(ipx-1,ipy+1)==otherplr) && (box(ipx-2,ipy+2)==0) && (box(ipx-3,ipy+3)==otherplr)
        check=1;
        box(ipx-1,ipy+1)=0;
        box(ipx-3,ipy+3)=0;
        xy=(ipx-1)*10+(ipy+1);
        set(ha(xy),'CData',blank);
        xy=(ipx-3)*10+(ipy+3);
        set(ha(xy),'CData',blank);
    elseif (fpx==ipx-4 && fpy==ipy-4) && ( box(ipx-1,ipy-1)==otherplr) && (box(ipx-2,ipy-2)==0) && (box(ipx-3,ipy-3)==otherplr)
        check=1;
        box(ipx-1,ipy-1)=0;
        box(ipx-3,ipy-3)=0;
        xy=(ipx-1)*10+(ipy-1);
        set(ha(xy),'CData',blank);
        xy=(ipx-3)*10+(ipy-3);
        set(ha(xy),'CData',blank);
    elseif (fpx==ipx+4 && fpy==ipy) && ( box(ipx+1,ipy+1)==otherplr) && (box(ipx+2,ipy+2)==0) && (box(ipx+3,ipy+1)==otherplr)
        check=1;
        box(ipx+1,ipy+1)=0;
        box(ipx+3,ipy+1)=0;
        xy=(ipx+1)*10+(ipy+1);
        set(ha(xy),'CData',blank);
        xy=(ipx+3)*10+(ipy+1);
        set(ha(xy),'CData',blank);
    elseif (fpx==ipx+4 && fpy==ipy) && ( box(ipx+1,ipy-1)==otherplr) && (box(ipx+2,ipy-2)==0) && (box(ipx+3,ipy-1)==otherplr)
        check=1;
        box(ipx+1,ipy-1)=0;
        box(ipx+3,ipy-1)=0;
        xy=(ipx+1)*10+(ipy-1);
        set(ha(xy),'CData',blank);
        xy=(ipx+3)*10+(ipy-1);
        set(ha(xy),'CData',blank);
    elseif (fpx==ipx-4 && fpy==ipy) && ( box(ipx-1,ipy+1)==otherplr) && (box(ipx-2,ipy+2)==0) && (box(ipx-3,ipy+1)==otherplr)
        check=1;
        box(ipx-1,ipy+1)=0;
        box(ipx-3,ipy+1)=0;
        xy=(ipx-1)*10+(ipy+1);
        set(ha(xy),'CData',blank);
        xy=(ipx-3)*10+(ipy+1);
        set(ha(xy),'CData',blank);
    elseif (fpx==ipx-4 && fpy==ipy) && ( box(ipx-1,ipy-1)==otherplr) && (box(ipx-2,ipy-2)==0) && (box(ipx-3,ipy-1)==otherplr)
        check=1;
        box(ipx-1,ipy-1)=0;
        box(ipx-3,ipy-1)=0;
        xy=(ipx-1)*10+(ipy-1);
        set(ha(xy),'CData',blank);
        xy=(ipx-3)*10+(ipy-1);
        set(ha(xy),'CData',blank);
    elseif (fpx==ipx && fpy==ipy+4) && ( box(ipx+1,ipy+1)==otherplr) && (box(ipx+2,ipy+2)==0) && (box(ipx+1,ipy+3)==otherplr)
        check=1;
        box(ipx+1,ipy+1)=0;
        box(ipx+1,ipy+3)=0;
        xy=(ipx+1)*10+(ipy+1);
        set(ha(xy),'CData',blank);
        xy=(ipx+1)*10+(ipy+3);
        set(ha(xy),'CData',blank);
    elseif (fpx==ipx && fpy==ipy+4) && ( box(ipx-1,ipy+1)==otherplr) && (box(ipx-2,ipy+2)==0) && (box(ipx-1,ipy+3)==otherplr)
        check=1;
        box(ipx-1,ipy+1)=0;
        box(ipx-1,ipy+3)=0;
        xy=(ipx-1)*10+(ipy+1);
        set(ha(xy),'CData',blank);
        xy=(ipx-1)*10+(ipy+3);
        set(ha(xy),'CData',blank);
    elseif (fpx==ipx && fpy==ipy-4) && ( box(ipx+1,ipy-1)==otherplr) && (box(ipx+2,ipy-2)==0) && (box(ipx+1,ipy-3)==otherplr)
        check=1;
        box(ipx+1,ipy-1)=0;
        box(ipx+1,ipy-3)=0;
        xy=(ipx+1)*10+(ipy-1);
        set(ha(xy),'CData',blank);
        xy=(ipx+1)*10+(ipy-3);
        set(ha(xy),'CData',blank);
    elseif (fpx==ipx && fpy==ipy-4) && ( box(ipx-1,ipy-1)==otherplr) && (box(ipx-2,ipy-2)==0) && (box(ipx-1,ipy-3)==otherplr)
        check=1;
        box(ipx-1,ipy-1)=0;
        box(ipx-1,ipy-3)=0;
        xy=(ipx-1)*10+(ipy-1);
        set(ha(xy),'CData',blank);
        xy=(ipx-1)*10+(ipy-3);
        set(ha(xy),'CData',blank);
    else
        check=0;
    end
end


text:
% data1: sea-water Pc
% Written by Ali A. Eftekhari
% Last checked: June 2021
data{1}=[0.073 249
0.085 45
0.09 15
0.131 8
0.197 3.33
0.23 2.22
0.25 1.79
0.328 0.6
0.38 0
0.5 -0.1
0.6 -0.3
0.7 -0.7
0.8 -2.5
0.831 -5.8
0.84 -10
0.856 -20
0.863 -32.2
0.867 -43
0.869 -49.1
0.87 -54.3
0.882 -92.8
0.893 -126.9
0.895 -142.5
0.899 -167.7
0.901 -250];
% zero sulphate Pc, table 3
data{2}=[0.056 200
0.058 36.9
0.073 9.3
0.112 5.06
0.137 3.56
0.176 1.46
0.232 0.5
0.28 0
0.503 -0.8
0.606 -1.13
0.698 -2
0.768 -3
0.793 -4
0.81 -5
0.827 -13
0.84 -31.4
0.857 -79.3];
% JCR library for cretaceous chalk
data{3}=[0.079 0.1
0.1 0.02
0.3 0.015
0.54 0.0
0.68 -0.4
0.747 -0.7];
% fit model to data
pc1=data{3}(:,2)*1e5; % Pa
sw1=data{3}(:,1);
sw_plot=linspace(min(sw1),max(sw1),5000);
% pc_plot=interp1(sw, pc, sw_plot, 'linear', 'extrap');
% plot(sw, pc, 'o', sw_plot, pc_plot)
% plot(diff(pc_plot)./diff(sw_plot))
% log(pc)=log(pce)-(1/labda)*log((sw-sw0)/(1-sw0-sor))
% ind0=find(pc==0, 1);
% p1=polyfit(sw(1:ind0-1), log(pc(1:ind0-1)), 1)
% plot(sw(1:ind0-1), log(pc(1:ind0-1)), 'o')

pc=@(sw, cwi, coi, awi, aoi, swc, sor)(cwi./((sw-swc)/(1-swc)).^awi+coi./((1-sw-sor)/(1-sor)).^aoi);
swc0=0.0789;
sor0=0.2529;
labda=2.4;
f=@(x, sw)pc(sw, x(1), x(2), 1/labda, 1/labda, swc0, sor0);
fw=@(x, sw)pc(sw, x(1), 0.0 , x(3), x(4), swc0, sor0);
fo=@(x, sw)pc(sw, 0.0, x(2), x(3), x(4), swc0, sor0);
x=lsqcurvefit(f, [1e3, -1e3],sw1, pc1)
x=patternsearch(@(x)(sum(abs(f(x, sw1)-pc1))), x)
plot(sw1, pc1, 'o', sw_plot, f(x, sw_plot))%, ...
% sw_plot, fw(x, sw_plot), '--', ...
% sw_plot, fo(x, sw_plot), '-.');
grid;

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function AddShortcutCategories(categories)
%ADDSHORTCUTCATEGORIES  Adds a new category of shortcut.
%
% If CATEGORY is a character array, ADDSHORTCUTCATEGORIES(CATEGORY) adds
% a new category of shortcut named CATEGORY.
%
% If CATEGORY is a cell array of strings,
% ADDSHORTCUTCATEGORIES(CATEGORY) adds a new category of shortcut for
% each string in CATEGORY.
%
% Examples:
%
% Provide sample usage code here
%
% See also: List related files here

% $Author: rcotton $	$Date: 2010/10/01 13:58:56 $	$Revision: 1.1 $
% Copyright: Health and Safety Laboratory 2010

if ischar(categories)
   categories = cellstr(categories);
end

shortcutUtils = GetShortcutUtils();
cellfun(@(category) shortcutUtils.addNewCategory(category), categories);

end

author: Sable


text:
function M = spm_mesh_join(Ms)
% Join a list of surface meshes into a single one
% FORMAT M = spm_mesh_join(Ms)
% Ms       - a patch structure array
%
% M        - a patch structure
%
% See also spm_mesh_split
%__________________________________________________________________________
% Copyright (C) 2016 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_join.m 6912 2016-10-26 14:50:53Z guillaume $

fn = fieldnames(Ms(1));
M  = cell2struct(cell(numel(fn),1),fn);

for i=1:numel(Ms)
    if i==1
        M.faces   = Ms(i).faces;
    else
        M.faces   = [M.faces; Ms(i).faces+size(M.vertices,1)];
    end
    M.vertices    = [M.vertices; Ms(i).vertices];
    try, M.cdata  = [M.cdata; Ms(i).cdata]; end
    try
        if i==1, M.mat = Ms(i).mat; end
    end
    if isfield(M,'mat')
        if sum(sum(M.mat - Ms(i).mat)) > 10*eps
            error('Meshes have different orientation.');
        end
    end
end

author: spm


text:
% circlecoords - returns the pixel coordinates of a circle defined by the
%                radius and x, y coordinates of its centre.
%
% Usage:
% [x,y] = circlecoords(c, r, imgsize,nsides)
%
% Arguments:
%	c           - an array containing the centre coordinates of the circle
%	              [x,y]
%   r           - the radius of the circle
%   imgsize     - size of the image array to plot coordinates onto
%   nsides      - the circle is actually approximated by a polygon, this
%                 argument gives the number of sides used in this approximation. Default
%                 is 600.
%
% Output:
%	x		    - an array containing x coordinates of circle boundary
%	              points
%   y		    - an array containing y coordinates of circle boundary
%                 points
%
% Author:
% Libor Masek
% masekl01@csse.uwa.edu.au
% School of Computer Science & Software Engineering
% The University of Western Australia
% November 2003

function [x,y] = circlecoords(c, r, imgsize,nsides)

    if nargin == 3
	nsides = 600;
    end

    nsides = round(nsides);

    a = [0:pi/nsides:2*pi];
    xd = (double(r)*cos(a)+ double(c(1)) );
    yd = (double(r)*sin(a)+ double(c(2)) );

    xd = round(xd);
    yd = round(yd);

    %get rid of -ves
    %get rid of values larger than image
    xd2 = xd;
    coords = find(xd>imgsize(2));
    xd2(coords) = imgsize(2);
    coords = find(xd<=0);
    xd2(coords) = 1;

    yd2 = yd;
    coords = find(yd>imgsize(1));
    yd2(coords) = imgsize(1);
    coords = find(yd<=0);
    yd2(coords) = 1;

    x = int32(xd2);
    y = int32(yd2);

author: Qingbao

repo: iris

save_path: github-repos/MATLAB/Qingbao-iris