text:
function cvt_test11 ( )

%*****************************************************************************80
%
%% CVT_TEST11 tests CVT.
%
%  Discussion:
%
%    In this test, we initialize the generators to grid points; this is
%    an unstable CVT solution.  The data would "prefer" to be in a
%    different form.  However, even if we take 2000 steps of CVT iteration,
%    the data is still only slowly progressing towards that other
%    configuration.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 November 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST11\n' );
  fprintf ( 1, '  CVT computes a Centroidal Voronoi Tessellation.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this test, we initialize the generators to\n' );
  fprintf ( 1, '  grid points; this is an unstable CVT solution.\n' );

  dim_num = 2;
  n = 16;
  batch = 1000;
  init = 4;
  init_string = 'user initialization';
  it_max = 40;
  it_fixed = 1;
  sample = 0;
  sample_num = 1000;
  sample_string = 'uniform';
  seed = 123456789;

  seed_init = seed;
%
%  Initialize the tuple generator.
%
  rank = -1;
  ngrid = 4;
  tuple(1:dim_num) = tuple_next_fast ( ngrid, dim_num, rank );
%
%  Pick points on a grid.
%
  for rank = 0 : n-1
    tuple(1:dim_num) = tuple_next_fast ( ngrid, dim_num, rank );
    r(1:dim_num,rank+1) = ( 2 * tuple(1:dim_num)' - 1 ) / ( 2 * ngrid );
  end

  r8mat_transpose_print ( dim_num, n, r, '  Initial generators (rows):' );

  [ r, seed, it_num, it_diff, energy ] = cvt ( dim_num, n, batch, init, ...
    sample, sample_num, it_max, it_fixed, seed, r );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Dimension DIM_NUM =        %12d\n', dim_num );
  fprintf ( 1, '  Number of points N =       %12d\n', n );
  fprintf ( 1, '  Initial SEED =             %12d\n', seed_init );
  fprintf ( 1, '  Current SEED =             %12d\n', seed );
  fprintf ( 1, '  INIT =                    "%s".\n', init_string );
  fprintf ( 1, '  Max iterations IT_MAX =    %12d\n', it_max );
  fprintf ( 1, '  IT_FIXED (fixed samples) = %12d\n', it_fixed );
  fprintf ( 1, '  Iterations IT_NUM =        %12d\n', it_num );
  fprintf ( 1, '  Difference IT_DIFF =       %14f\n', it_diff );
  fprintf ( 1, '  CVT ENERGY =               %14f\n', energy );
  fprintf ( 1, '  SAMPLE =                  "%s".\n', sample_string );
  fprintf ( 1, '  Samples SAMPLE_NUM    =    %12d\n', sample_num );
  fprintf ( 1, '  Sampling BATCH size =      %12d\n', batch );
  fprintf ( 1, '  EPSILON (unit roundoff) =  %12e\n', eps );

  r8mat_transpose_print ( dim_num, n, r, '  Generators (rows):' );

  return
end


text:
function OUT = datatreat(DATA,vtreat,nlag)
% =======================================================================
% Treat a time series with the specified method. If changes are computed
% the function assumes 1 lag (unless otherwise specified)
% =======================================================================
% OUT = datatreat(DATA,vtreat,lag)
% -----------------------------------------------------------------------
% INPUT
%	- DATA: matrix DATA(T,N)
%   - vtreat: 0 No treatment
%             1	Log
%             2 Difference
%             3	Log Difference
%             4 Percent change
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%   - lag: number of lags for changes
% -----------------------------------------------------------------------
% OUTPUT
%	- OUT: matrix DATA(T,N) of treated data. The first "nlag"
%     observations are NaNs
% -----------------------------------------------------------------------
% EXAMPLE
%   x = [1 2; -3 4; 5 6; 7 8; 9 10];
%   OUT = datatreat(x,3)
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

% Check input
if ~exist('nlag','var')
    nlag = 1;
end

% Set matrices
[nobs, nvar] = size(DATA);
OUT = nan(nobs,nvar);

% No treatment
if vtreat==0
    OUT = DATA;

% Log
elseif vtreat==1
    if ~isempty(DATA<0), warning('Negative numbers set to NaN before taking logs'), end
    DATA(DATA<0) = NaN;
    OUT = log(DATA);

% Difference
elseif vtreat==2
    OUT(nlag+1:end,:) = DATA(1+nlag:end,:) - DATA(1:end-nlag,:);

% Log difference
elseif vtreat==3
    if ~isempty(DATA<0), warning('Negative numbers set to NaN before taking logs'), end
    DATA(DATA<0) = NaN;
    OUT(nlag+1:end,:) = log(DATA(1+nlag:end,:))-log(DATA(1:end-nlag,:));

% Percent change
elseif vtreat==4
    OUT(nlag+1:end,:) = DATA(1+nlag:end,:)./DATA(1:end-nlag,:)-1;
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
classdef TotalThrustConstraint < AbstractConstraint
    %TotalThrustConstraint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function obj = TotalThrustConstraint(event, lb, ub)
            obj.event = event;
            obj.lb = lb;
            obj.ub = ub;

             obj.id = rand();
        end

        function [lb, ub] = getBounds(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, celBodyData)
            type = obj.getConstraintType();

            switch obj.eventNode
                case ConstraintStateComparisonNodeEnum.FinalState
                    stateLogEntry = stateLog.getLastStateLogForEvent(obj.event);

                case ConstraintStateComparisonNodeEnum.InitialState
                    stateLogEntry = stateLog.getFirstStateLogForEvent(obj.event);

                otherwise
                    error('Unknown event node.');
            end

            totalThrust = TotalThrustConstraint.getTotalThrust(stateLogEntry);
            value = totalThrust;

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                switch obj.stateCompNode
                    case ConstraintStateComparisonNodeEnum.FinalState
                        stateLogEntryStateComp = stateLog.getLastStateLogForEvent(obj.stateCompEvent);

                    case ConstraintStateComparisonNodeEnum.InitialState
                        stateLogEntryStateComp = stateLog.getFirstStateLogForEvent(obj.stateCompEvent);

                    otherwise
                        error('Unknown event node.');
                end

                valueStateComp = TotalThrustConstraint.getTotalThrust(stateLogEntryStateComp); %it's probably not a bad idea to leave the state here in the original reference frame.
            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

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

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

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

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

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

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

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

        function tf = canUseSparseOutput(obj)
            tf = true;
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function type = getConstraintType(obj)
            type = 'Total Thrust';
        end

%         function name = getName(obj)
%             name = sprintf('%s - Event %i', obj.getConstraintType(), obj.event.getEventNum());
%         end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = 'kN';
            lbLim = 0;
            ubLim = Inf;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
%             addConstraintTf = lvd_EditGenericMAConstraintGUI(obj, lvdData);

            output = AppDesignerGUIOutput({false});
            lvd_EditGenericMAConstraintGUI_App(obj, lvdData, output);
            addConstraintTf = output.output{1};
        end
    end

    methods(Static)
        function totalThrust = getTotalThrust(stateLogEntry)
            ut = stateLogEntry.time;
            rVect = stateLogEntry.position;
            vVect = stateLogEntry.velocity;

            bodyInfo = stateLogEntry.centralBody;
            tankStates = stateLogEntry.getAllActiveTankStates();
            stageStates = stateLogEntry.stageStates;
            lvState = stateLogEntry.lvState;

            dryMass = stateLogEntry.getTotalVehicleDryMass();
            tankStatesMasses = [tankStates.tankMass]';

            throttleModel = stateLogEntry.throttleModel;
            steeringModel = stateLogEntry.steeringModel;

            attState = LaunchVehicleAttitudeState();
            attState.dcm = steeringModel.getBody2InertialDcmAtTime(ut, rVect, vVect, bodyInfo);

            altitude = norm(rVect) - bodyInfo.radius;
            pressure = getPressureAtAltitude(bodyInfo, altitude);

            powerStorageStates = stateLogEntry.getAllActivePwrStorageStates();
            storageSoCs = NaN(size(powerStorageStates));
            for(i=1:length(powerStorageStates)) %#ok<*NO4LP>
                storageSoCs(i) = powerStorageStates(i).getStateOfCharge();
            end

            throttle = throttleModel.getThrottleAtTime(ut, rVect, vVect, tankStatesMasses, dryMass, stageStates, lvState, tankStates, bodyInfo, storageSoCs, powerStorageStates);

            [~, totalThrust, ~] = LaunchVehicleStateLogEntry.getTankMassFlowRatesDueToEngines(tankStates, tankStatesMasses, stageStates, throttle, lvState, pressure, ut, rVect, vVect, bodyInfo, steeringModel, storageSoCs, powerStorageStates, attState);
        end
    end

    methods(Static)
        function constraint = getDefaultConstraint(~, ~)
            constraint = TotalThrustConstraint(LaunchVehicleEvent.empty(1,0),0,0);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%% SOBOL_TEST04 tests SOBOL.
%
global SOBOL_lastq;
global SOBOL_seed;

fprintf ( 1, '\n' );
fprintf ( 1, 'SOBOL_TEST04\n' );
fprintf ( 1, '  SOBOL returns the next element\n' );
fprintf ( 1, '  of a Sobol sequence.\n' );
fprintf ( 1, '\n' );
fprintf ( 1, '  In this test, we call Sobol repeatedly.\n' );

dim_max = 4;

for ( dim_num = 2 : dim_max )

  SOBOL_seed = 0;
  seed = 0;
  qs = prime_ge ( dim_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Using dimension DIM_NUM =   %d\n', dim_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Seed  Seed   Sobol\n' );
  fprintf ( 1, '  In    Out\n' );
  fprintf ( 1, '\n' );
  for ( i = 0 : 110 )
    [ r, seed_out ] = sobol ( dim_num, seed );
    if ( i <= 11 || 95 <= i )
      fprintf ( 1, '%6d %6d  ', seed, seed_out );
      for ( j = 1 : dim_num )
        fprintf ( 1, '%10f  ', r(j) );
      end
      fprintf ( 1, '\n' );
    elseif ( i == 12 )
      fprintf ( 1, '......................\n' );
    end
    seed = seed_out;
  end

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [CTTraceS, RTOGPBVectorsM, RTOGPBVectorsM_MC, PBMaskM, rowPBV, colPBV, xPBPosV, yPBPosV, beamlet_delta_x, beamlet_delta_y] = ...
   getPBRayData_cerr3(sourceS, numSamplePts, xySampleRate,planC, xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, gA);
% function [CTTraceS, RTOGPBVectorsM, RTOGPBVectorsM_MC, PBMaskM, rowPBV, colPBV, xPBPosV, yPBPosV, beamlet_delta_x, beamlet_delta_y] = ...
%    getPBRayData_cerr3(edgeS, sourceS, numSamplePts, xySampleRate,planC, xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, gA);
%JOD
% global planC
%
% indexS=planC{end};

water = 1000; %Assumes water equals 1000.

%-----------Get CT scan---------------------%

% [CTUniform3D, CTUniformInfoS] = getUniformizedCTScan(1,planC);
scanNum = 1;
[CTUniform3D, CTUniformInfoS] = getUniformizedCTScan(1,scanNum);

xOffset = CTUniformInfoS.xOffset;
yOffset = CTUniformInfoS.yOffset;

%-----------Fix source characteristics---------------------%

orgV   = [sourceS.isocenter.x, sourceS.isocenter.y, sourceS.isocenter.z];

%-----------Get ray parameters---------------------%

% [RTOGPBVectorsM, RTOGPBVectorsM_MC, PBMaskM, rowPBV, colPBV, xPBPosV, yPBPosV] = ...
%       getPBRays(edgeS, sourceS, xySampleRate);
% [RTOGPBVectorsM, RTOGPBVectorsM_MC, PBMaskM, rowPBV, colPBV, xPBPosV, yPBPosV, beamlet_delta_x, beamlet_delta_y] = ...
%       getPBRays(xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, gA);
[RTOGPBVectorsM, RTOGPBVectorsM_MC, PBMaskM, rowPBV, colPBV, xPBPosV, yPBPosV, beamlet_delta_x, beamlet_delta_y] = ...
      getPBRays_planChk(xPosV, yPosV, beamlet_delta_x, beamlet_delta_y, sourceS);

rayLength = 500;  %in cm.  This is the length of the rays which are passed
                  %through the CT matrix to determine cumulative CT densities.

numSlices = size(CTUniform3D,3);

zFirst = CTUniformInfoS.firstZValue;

sliceThickness = CTUniformInfoS.sliceThickness;

%minBox an maxBox containt coordinates of the corners of the CT box as used
%by the ray intersection routine.

maxBoxS.z = (numSlices - 1) * sliceThickness + zFirst + sliceThickness/2;

minBoxS.z = zFirst - sliceThickness/2;

delta_xy = CTUniformInfoS.grid1Units;

imageWidth = CTUniformInfoS.size(1);

minBoxS.x = - imageWidth/2 * delta_xy + xOffset;
maxBoxS.x =   imageWidth/2 * delta_xy + xOffset;

minBoxS.y = - imageWidth/2 * delta_xy + yOffset;
maxBoxS.y =   imageWidth/2 * delta_xy + yOffset;

imageWidth = CTUniformInfoS.sizeOfDimension1;

CTTraceS = struct('CTNumsRay',[],'CTCumNumsRay',[],'distSamplePts',[]);

for i = 1 : size(RTOGPBVectorsM,1)

  rayDeltaS.x = RTOGPBVectorsM(i,1) * rayLength;
  rayDeltaS.y = RTOGPBVectorsM(i,2) * rayLength;
  rayDeltaS.z = RTOGPBVectorsM(i,3) * rayLength;

  %are the components of the ray's direction and maximum length.

  deltaV = [rayDeltaS.x, rayDeltaS.y, rayDeltaS.z];

  t_entrance = rayBoxIntersection(sourceS,rayDeltaS,minBoxS,maxBoxS);

  %The entrance point is therefore
  entranceV = orgV + t_entrance * deltaV;

  %Reflect to find exit point (assume length of ray is long enough that ray does exit):

  if t_entrance ~= -1

    %find exit point
    %get end of ray
    rayOrgS2.xRel = sourceS.xRel + rayDeltaS.x;  %reflected source positions
    rayOrgS2.yRel = sourceS.yRel + rayDeltaS.y;
    rayOrgS2.zRel = sourceS.zRel + rayDeltaS.z;

    rayDeltaS2.x = - rayDeltaS.x;
    rayDeltaS2.y = - rayDeltaS.y;
    rayDeltaS2.z = - rayDeltaS.z;

    rayOrgS2.isocenter = sourceS.isocenter;

    t = rayBoxIntersection(rayOrgS2,rayDeltaS2,minBoxS,maxBoxS);

    t_exit = 1 - t;

    exitV = orgV + t_exit * deltaV;

else
    error('PB Ray does not intersect CT scan.');
end

  %Now produce a set of sampling points between the entrance and the exit

  nV = 1 : numSamplePts;

  delta_t  = (t_exit - t_entrance)/(numSamplePts - 1);

  tV = t_entrance + (nV - 1) * delta_t;

  CTTraceS(i).distSamplePts = tV * sum(deltaV.^2).^0.5;

  sampleV.x =  sourceS.xRel + tV * deltaV(1);

  sampleV.y =  sourceS.yRel + tV * deltaV(2);

  sampleV.z =  sourceS.zRel + tV * deltaV(3);

  sampleRTOGV.x = sampleV.x + sourceS.isocenter.x;
  sampleRTOGV.y = sampleV.y + sourceS.isocenter.y;
  sampleRTOGV.z = sampleV.z + sourceS.isocenter.z;

  %---------Sample CT densities----------%

  %To go from sample points in RTOG system to CT densities, we convert as follows:

  %What is the slice number?
  sliceV = 1 + (sampleRTOGV.z - zFirst)/sliceThickness;

  %Now do 3-D interpolation:
  zFieldV = [minBoxS.z + 0.5 * sliceThickness:sliceThickness:maxBoxS.z - 0.5 * sliceThickness];
  xFieldV = [minBoxS.x + 0.5 * delta_xy, delta_xy, maxBoxS.x - 0.5 * delta_xy];
  yFieldV = [minBoxS.y + 0.5 * delta_xy, delta_xy, maxBoxS.y - 0.5 * delta_xy];

  [CTNumsV] = finterp3(sampleRTOGV.x, sampleRTOGV.y, sampleRTOGV.z, CTUniform3D, xFieldV, yFieldV, zFieldV, 0);

%   [CTNumsV] = finterp3(sampleRTOGV.x, sampleRTOGV.y, sampleRTOGV.z, CTUniform3D, xFieldV, yFieldV, zFieldV, 0);

  CTTraceS(i).densityRay = delta_t * norm(deltaV) * CTNumsV/water;
  CTTraceS(i).cumDensityRay = cumsum(CTTraceS(i).densityRay);  %Account for sampling rate to convert to g/cm^2.

  if any(CTTraceS(i).cumDensityRay > 50)
      warning('Cumulative density ray appears to exceed maximum length.');
  end

end

%-----------fini---------------------%

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function cS = monazite

%
cs_Mnz = crystalSymmetry('2/m',[0.967 1 0.921], [90, 103.4, 90]*degree);
N = Miller({1,0,0},{0,1,0},{1,1,0},{0,1,1},{1,1,-1},{0,2,1},{0,1,2},{1,0,-1},{1,0,1},{3,-1,-1},{2,-1,0},{0,-1,-1},{1,-2,-1},cs_Mnz);
dist = [3.6, 4.95, 6.5, 10, 8.2, 10.65, 11.1, 6.75, 7.45, 14.65, 9.9, 6.15, 12.2];
cS = crystalShape(N./dist);

author: mtex-toolbox

repo: mtex

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


text:
clear variables;
EXPORT_DIR = '/home/mhughes/git/BPARHMM-NEW/figs/Mocap6/';

%jobID = 1121860; taskID = 'best';
%jobID = 2134873; taskID = 1;
jobID = 111; taskID = 1;

% ------------------------------------------------------- Load Ground truth
F{1} = zeros(6, 12 );
Q = loadSamplerInfo( jobID, taskID);
for ii = 1:Q.data.N
   activeFeatIDs =  unique( Q.data.zTrue(ii) );
   F{1}(ii, activeFeatIDs) = 1;
end

% ------------------------------------------------------- Load best BPHMM
%  Note that meanFfrac is aggregated across all post-burn-in samples
MIN_THR = 0.02;
% disregard feature if appeared active in less that 2%
%     of state sequence assignment zs when avg'd across all samples
X = loadSamplerOutput( jobID, taskID );
F{2} =  X.Summary.meanFfrac > 0.02;

if size( F{2}, 2 ) > 12
    if all( sum( F{2}(:,13:end), 1 ) == 0 )
        F{2} = F{2}(:,1:12);
    end
end

doTHR = 0;
MINITER = 10000;

keepiters = find( X.iters.Psi >= MINITER );

avgF = zeros( 6, 15 );
for id = keepiters
    K = size( X.A( id ).Ffrac, 2 );
    Fcur = zeros( size(avgF) );
    if doTHR
        Fcur(:,1:K) = X.A( id ).Ffrac > 0.02;
    else
        Fcur(:, 1:K) = X.A( id ).F;
    end
    avgF = avgF + Fcur;
end
avgF = avgF ./ length(keepiters );

kIDs = union( 1:12, find( sum( avgF, 1) > 0) );
F{2} = avgF(:, kIDs);

% ------------------------------------------------------- Load GMM/HMM
Q = load( '/data/liv/mhughes/data/MoCap6/Results/GMM_25trials_K2-20.mat' );
F{3} = Q.GMM(2,12).alignedF;
Q = load( '/data/liv/mhughes/data/MoCap6/Results/HMM_25trials_K2-20.mat' );
F{4} = Q.HMM(2,12).alignedF;

FigNames = {'GroundTruth','BPHMM','GMM','HMM'};
for ff = 1:length( F )
figure;
imagesc( F{ff} );
set( gcf, 'Name', ['Mocap6_FeatMatrix_' FigNames{ff} ]);
colormap( 'bone' );
set( gca, 'FontSize', 20 );

if ff==2
    rStr = num2str(doTHR);
else
    rStr = '';
end
set( gcf, 'InvertHardcopy','off','Color',[1 1 1]);

  export_fig( fullfile( EXPORT_DIR, ['Mocap6_FeatMatrix_' rStr FigNames{ff} ]), '-eps');
end

author: michaelchughes

repo: NPBayesHMM


text:
function [vw_nuisance, vw_nuisance_comps] = canlab_extract_ventricle_wm_timeseries(mask_image_dir, imgs, varargin)

% function [nuisance, nuisance_comps] = canlab_extract_ventricle_wm_timeseries(mask_image_dir)
%
% This function extracts data from ventricles and white matter for each
% time series. You can include vw_nuisance in your nuisance matrix in the
% first level model estimation.
%
% output: "vw_nuisance" returns the MEAN ventricle and white matter signal
%         1st column: ventricle signal, 2nd column: white matter signal,
%         rows: time series
%
%         "vw_nuisance_comps" returns first 5 components from PCA.
%
% input:
% mask_image_dir: a directory that has mask files for ventricle and white matter.
%        You can use "canlab_create_wm_ventricle_masks.m" to get these mask files.
%        The mask file names must be "ventricles.img' and 'white_matter.img'.
%        e.g) mask_image_dir = fullfile(subjdir, 'Structural/SPGR');
%
% imgs: image file names that you want to extract data from.
%        e.g) imgs = filenames(fullfile(subjdir, 'Functional/Preprocessed/run1/wra*nii'), 'absolute', 'char');
%
% option (varargin): 'noplot' - the default setting is plotting output
%                    variables. You can use 'noplot' option if you don't want the plot.
%
% 5/4/2012 by Tor Wager and Wani Woo

%% --------------------------------------------
% Extract data from ventricles
% ---------------------------------------------

doplot = true;

wh = strcmp(varargin, 'noplot');
if any(wh), doplot = false; varargin(wh) = []; end

vmask = filenames(fullfile(mask_image_dir, 'ventricles.img'), 'char', 'absolute');
wmask = filenames(fullfile(mask_image_dir, 'white_matter.img'), 'char', 'absolute');

vdat = fmri_data(imgs, vmask);

vw_nuisance(:, 1) = mean(vdat.dat, 1)'; % avg across voxels

% save the first 5 components

[~, nuisance_comps] = pca(vdat.dat', 'econ');

vw_nuisance_comps = nuisance_comps(:, 1:5);

% ---------------------------------------------
% Extract data from white matter
% ---------------------------------------------

wdat = fmri_data(imgs, wmask);

vw_nuisance(:, 2) = mean(wdat.dat, 1)'; % avg across voxels

% save the first 5 components

[~, n2] = pca(wdat.dat', 'econ');

vw_nuisance_comps = [vw_nuisance_comps n2(:, 1:5)];

%% plot

if doplot

    create_figure('vw_nuisance', 3, 1)

    plot(zscore(vw_nuisance));
    title('averages (z-scores)')

    subplot(3, 1, 2)

    plot(vw_nuisance_comps(:, 1:5))
    title('ventricle components')

    subplot(3, 1, 3)

    plot(vw_nuisance_comps(:, 6:10))
    title('white matter components')
end

return


text:
function features = LoadSpikeFeatures(filename,rate)

%LoadSpikeFeatures - Load spike features from file.
%
%  USAGE
%
%    features = LoadSpikeFeatures(filename,rate)
%
%    filename            spike file name (either .clu, .res or .fet)
%    rate                sampling rate
%
%  OUTPUT
%
%    The output is a list of (timestamp,group,cluster,features...) t-uples.
%
%  SEE
%
%    See also GetSpikeFeatures.

% Copyright (C) 2004-2011 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.

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

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

% Load .clu file
filename = [path '/' basename '.clu.' int2str(electrodeGroup)];
if ~exist(filename),
	error(['File ''' filename ''' not found.']);
end
clu = load(filename);
clu = clu(2:end);

% Load .fet file
filename = [path '/' basename '.fet.' int2str(electrodeGroup)];
if ~exist(filename),
	error(['File ''' filename ''' not found.']);
end
file = fopen(filename,'r');
if file == -1,
	error(['Cannot open file ''' filename '''.']);
end
nFeatures = fscanf(file,'%d',1);
fet = fscanf(file,'%f',[nFeatures,inf])';
fclose(file);

features = [fet(:,end)/rate electrodeGroup*ones(size(clu)) clu fet(:,1:end-1)];

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [varargout] = DES(input64,mode,key)
%DES: Data Encryption Standard
% Encrypt/Decrypt a 64-bit message using a 64-bit key using the Feistel Network
% -------------------------------------------------------------------------
% Inputs:
%        input64 = a 64-bit message
%           mode = either 'ENC' encryption or 'DEC' decryption (default 'ENC')
%            key = a 56/64-bit key (optional under 'ENC', but mandatory under 'DEC')
% Outputs:
%   varargout{1} = output64, a 64-bit message after encryption/decryption
%   varargout{2} = a 64-bit key, if a 64-bit key is not provided as an input
% -------------------------------------------------------------------------
% Demos:
%   plaintext = round(rand(1,64));
%   [ciphertext,key] = DES(plaintext);       % Encryption syntex 1
%   [ciphertext1,key] = DES(plaintext,'ENC'); % Encryption syntex 2
%   deciphertext1 = DES(ciphertext1,'DEC',key);% Decryption syntex
%
%   key56 = round(rand(1,56));
%   [ciphertext2,key64] = DES(plaintext,'ENC',key56);% Encryption syntex 3 (56-bit key)
%   deciphertext2 = DES(ciphertext2,'DEC',key64);     % Decryption syntex   (64-bit key)
%   ciphertext3 = DES(plaintext,'ENC',key64);       % Encryption syntex 3 (64-bit key)
%   deciphertext3 = DES(ciphertext3,'DEC',key56);     % Decryption syntex   (56-bit key)
%
%   % plot results
%   subplot(4,2,1),plot(plaintext),ylim([-.5,1.5]),xlim([1,64]),title('plaintext')
%   subplot(4,2,2),plot(ciphertext),ylim([-.5,1.5]),xlim([1,64]),title('ciphertext')
%   subplot(4,2,3),plot(deciphertext1),ylim([-.5,1.5]),xlim([1,64]),title('deciphertext1')
%   subplot(4,2,4),plot(ciphertext1),ylim([-.5,1.5]),xlim([1,64]),title('ciphertext1')
%   subplot(4,2,5),plot(deciphertext2),ylim([-.5,1.5]),xlim([1,64]),title('deciphertext2')
%   subplot(4,2,6),plot(ciphertext2),ylim([-.5,1.5]),xlim([1,64]),title('ciphertext2')
%   subplot(4,2,7),plot(deciphertext3),ylim([-.5,1.5]),xlim([1,64]),title('deciphertext3')
%   subplot(4,2,8),plot(ciphertext3),ylim([-.5,1.5]),xlim([1,64]),title('ciphertext3')
% -------------------------------------------------------------------------
% NOTE:
% 1. If a 64-bit key is provided, then its bit parities will be checked. If
%    a 56-bit key is provided, then it is automatically added 8 partity
%    checking bits. However, the 8 parity bits are never used in
%    DES encryption/decryption process. They are included just for the
%    completeness of a DES implementation.
% 2. Cipher modes are not provided in this simple script. If you are
%    interested or do not know what does cipher modes mean, please go to page
%    http://en.wikipedia.org/wiki/Block_cipher_modes_of_operation
%    for details. Please keep in mind that selecting an inappropriate working
%    mode may extremely weaken the security of your messages.
% 3. A general description of DES can be found at its wiki page:
%    http://en.wikipedia.org/wiki/Data_Encryption_Standard
%    The detailed cryptographical primitives can be found under the page:
%    http://en.wikipedia.org/wiki/DES_supplementary_material
%    If you want to speed-up the DES code here, you can simply store these
%    primitives in memory and call them when you need.
% -------------------------------------------------------------------------
% By Yue (Rex) Wu
% ECE Dept @ Tufts Univ.
% 08/18/2012
% If you find bugs, please email me via ywu03@ece.tufts.edu
% -------------------------------------------------------------------------

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                           0. Initialization                           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 0.1 check input
error(nargchk(1,3,nargin));
switch nargin
    case 1
        mode = 'ENC';
        K = round(rand(8,7));
        K(:,8) = mod(sum(K,2),2); % note these eight bits of key are never used in encryption
        K = reshape(K',1,64);
        varargout{2} = K;
    case 2
        switch mode
            case 'ENC'
                K = round(rand(8,7));
                K(:,8) = mod(sum(K,2),2); % note these eight bits of key are never used in encryption
                K = reshape(K',1,64);
                varargout{2} = K;
            case 'DEC'
                error('Key has to be provided in decryption mode (DEC)')
            otherwise
                error('WRONG working mode!!! Select either encrtyption mode: ENC or decryption mode: DEC !!!')
        end
    case 3
        if isempty(setdiff(unique(key),[0,1])) % check provided key type
            if numel(key) == 64  % check provided key parity
                keyParityCheck = @(k) (sum(mod(sum(reshape(k,8,8)),2))==0);
                if keyParityCheck(key) == 1
                    K = key(:)';
                else
                    error('Key parity check FAILED!!!')
                end
            elseif numel(key) == 56 % add parity bits
                K = reshape(key,7,8)';
                K(:,8) = mod(sum(K,2),2); % note these eight bits of key are never used in encryption
                K = reshape(K',1,64);
                varargout{2} = K;
                display('Key parity bits added')
            else
                error('Key has to be either 56 or 64-bit long!!!')
            end
        else
            error('Key has to be binary!!!')
        end
end

% 0.2 check message length and type
if numel(input64) == 64 && isempty(setdiff(unique(input64),[0,1]))
    P = input64;
else
    error('Message has to be a 64-bit message!!!')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                   1. Cryptographical primitives                       %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1.1 define splitting function
HALF_L = @(message) message(1:32);
HALF_R = @(message) message(33:64);
% 1.2 define expansion function
EF = @(halfMessage) [halfMessage([32,4:4:28])',(reshape(halfMessage,4,8))',halfMessage([5:4:29,1])'];
% 1.3 define key mixing (KM)
KM = @(expandedHalfMessage,rK) xor(expandedHalfMessage,reshape(rK,6,8)');
% 1.4 define eight substitution tables
% input: 0	1   2   3   4   5   6   7   8   9   10  11  12  13  14  15
st{1} = [14	4	13	1	2	15	11	8	3	10	6	12	5	9	0	7;...
         0  15	7	4	14	2	13	1	10	6	12	11	9	5	3	8;...
         4	1	14	8	13	6	2	11	15	12	9	7	3	10	5	0;...
         15	12	8	2	4	9	1	7	5	11	3	14	10	0	6	13];
st{2} = [15	1	8	14	6	11	3	4	9	7	2	13	12	0	5	10;...
    	3	13	4	7	15	2	8	14	12	0	1	10	6	9	11	5;...
		0	14	7	11	10	4	13	1	5	8	12	6	9	3	2	15;...
		13	8	10	1	3	15	4	2	11	6	7	12	0	5	14	9];
st{3} = [10	0	9	14	6	3	15	5	1	13	12	7	11	4	2	8;...
		13	7	0	9	3	4	6	10	2	8	5	14	12	11	15	1;...
		13	6	4	9	8	15	3	0	11	1	2	12	5	10	14	7;...
		1	10	13	0	6	9	8	7	4	15	14	3	11	5	2	12];
st{4} = [7	13	14	3	0	6	9	10	1	2	8	5	11	12	4	15;...
		13	8	11	5	6	15	0	3	4	7	2	12	1	10	14	9;...
		10	6	9	0	12	11	7	13	15	1	3	14	5	2	8	4;...
		3	15	0	6	10	1	13	8	9	4	5	11	12	7	2	14];
st{5} = [2	12	4	1	7	10	11	6	8	5	3	15	13	0	14	9;...
		14	11	2	12	4	7	13	1	5	0	15	10	3	9	8	6;...
		4	2	1	11	10	13	7	8	15	9	12	5	6	3	0	14;...
		11	8	12	7	1	14	2	13	6	15	0	9	10	4	5	3];
st{6} = [12	1	10	15	9	2	6	8	0	13	3	4	14	7	5	11;...
		10	15	4	2	7	12	9	5	6	1	13	14	0	11	3	8;...
		9	14	15	5	2	8	12	3	7	0	4	10	1	13	11	6;...
		4	3	2	12	9	5	15	10	11	14	1	7	6	0	8	13];
st{7} = [4	11	2	14	15	0	8	13	3	12	9	7	5	10	6	1;...
		13	0	11	7	4	9	1	10	14	3	5	12	2	15	8	6;...
		1	4	11	13	12	3	7	14	10	15	6	8	0	5	9	2;...
		6	11	13	8	1	4	10	7	9	5	0	15	14	2	3	12];
st{8} = [13	2	8	4	6	15	11	1	10	9	3	14	5	0	12	7;...
		1	15	13	8	10	3	7	4	12	5	6	11	0	14	9	2;...
		7	11	4	1	9	12	14	2	0	6	10	13	15	3	5	8;...
		2	1	14	7	4	10	8	13	15	12	9	0	3	5	6	11];
% the eight binary s-boxes
for i = 1:8
    ST{i} = mat2cell(blkproc(st{i},[1,1],@(x) de2bi(x,4,'left-msb')),ones(1,4),ones(1,16)*4);
end
% 1.5 define subsitution function (SBOX)
SUBS = @(expandedHalfMessage,blkNo) ST{blkNo}{bi2de(expandedHalfMessage(blkNo,[1,6]),'left-msb')+1,bi2de(expandedHalfMessage(blkNo,[2:5]),'left-msb')+1};
SBOX = @(expandedHalfMessage) [SUBS(expandedHalfMessage,1);SUBS(expandedHalfMessage,2);...
                               SUBS(expandedHalfMessage,3);SUBS(expandedHalfMessage,4);...
                               SUBS(expandedHalfMessage,5);SUBS(expandedHalfMessage,6);...
                               SUBS(expandedHalfMessage,7);SUBS(expandedHalfMessage,8)];
% 1.6 define permutation function (PBOX)
PBOX = @(halfMessage) halfMessage([16  7 20 21  29 12 28 17 ...
                                    1 15 23 26   5 18 31 10 ...
                                    2  8 24 14  32 27  3  9 ...
                                   19 13 30  6  22 11  4  25]);
% 1.7 define initial permutation (IP)
IP = @(message) message([58	50	42	34	26	18	10	2 ...
                        60	52	44	36	28	20	12	4 ...
                        62	54	46	38	30	22	14	6 ...
                        64	56	48	40	32	24	16	8 ...
                        57	49	41	33	25	17	9	1 ...
                        59	51	43	35	27	19	11	3 ...
                        61	53	45	37	29	21	13	5 ...
                        63	55	47	39	31	23	15	7]);
% 1.8 define final permutation (FP)
FP = @(message) message([40	8	48	16	56	24	64	32 ...
                        39	7	47	15	55	23	63	31 ...
                        38	6	46	14	54	22	62	30 ...
                        37	5	45	13	53	21	61	29 ...
                        36	4	44	12	52	20	60	28 ...
                        35	3	43	11	51	19	59	27 ...
                        34	2	42	10	50	18	58	26 ...
                        33	1	41	9	49	17	57	25]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                           2. key schedule                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2.1 define permuted choice 1 (PC1)
PC1L = @(key64) key64([57	49	41	33	25	17	9 ...
                    1	58	50	42	34	26	18 ...
                    10	2	59	51	43	35	27 ...
                    19	11	3	60	52	44	36]);
PC1R = @(key64) key64([63	55	47	39	31	23	15 ...
                    7	62	54	46	38	30	22 ...
                    14	6	61	53	45	37	29 ...
                    21	13	5	28	20	12	4]);
% 2.2 define permuted choice 2 (PC2)
PC2 = @(key56) key56([14 17	11	24	1	5	3	28 ...
                     15	6	21	10	23	19	12	4 ...
                     26	8	16	7	27	20	13	2 ...
                     41	52	31	37	47	55	30	40 ...
                     51	45	33	48	44	49	39	56 ...
                     34	53	46	42	50	36	29	32]);
% 2.3 define rotations in key-schedule (RK)
% round# 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6
   RK = [1 1 2 2 2 2 2 2 1 2 2 2 2 2 2 1];
% 2.4 define key shift function (KS)
KS = @(key28,s) [key28(s+1:end),key28(1:s)];
% 2.5 define sub-keys for each round
leftHKey = PC1L(K); % 28-bit half key
rightHKey = PC1R(K);% 28-bit half key
for i = 1:16
    leftHKey = KS(leftHKey,RK(i));
    rightHKey = KS(rightHKey,RK(i));
    key56 = [leftHKey ,rightHKey];
    subKeys(i,:) = PC2(key56(:));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                           3. DES main loop                            %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.1 initial permutation
C = IP(P);
switch mode
    case 'ENC' % if encryption, split 64 message to two halves
        L{1} = HALF_L(C); % left-half 32-bit
        R{1} = HALF_R(C); % right-half 32-bit
    case 'DEC' % if decryption, swapping two halves
        L{1} = HALF_R(C);
        R{1} = HALF_L(C);
end
% 3.2 cipher round 1 to 16
for i = 1:16
     L{i+1} = R{i}; % half key: 32-bit
     expended_R = EF(R{i}); % expended half key: 32-bit to 48-bit
     switch mode
        case 'ENC' % if encryption, apply sub-keys in the original order
            mixed_R = KM(expended_R,subKeys(i,:)); % mixed with sub-key: 48-bit
        case 'DEC' % if decryption, apply sub-keys in the reverse order
            mixed_R = KM(expended_R,subKeys(16-i+1,:)); % mixed with sub-key: 48-bit
     end
     substituted_R = SBOX(mixed_R); % substitution: 48-bit to 32-bit
     permuted_R = PBOX(reshape(substituted_R',1,32)); % permutation: 32-bit
     R{i+1} = xor(L{i},permuted_R); % Feistel function: 32-bit
end
% 3.3 final permutation
switch mode
    case 'ENC'
        C = [L{end},R{end}];
    case 'DEC'
        C = [R{end},L{end}];
end
output64 = FP(C);
varargout{1} = output64;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                   END                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: Sable

repo: mcbench-benchmarks


text:
%%********************************************************************
%% blkbarrier: calculate
%% [-v(p)*logdet(X{p}),   v(p)*logdet(Z{p}) + n*v(p)*(1-log(v(p)))]
%% [-v(p)*log(gam(X{p})), v(p)*log(gam(Z{p})) + v(p)]
%% [-v(p)*log(X{p}),      v(p)*log(Z{p}) + n*v(p)*(1-log(v(p)))]
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

  function objadd = blkbarrier(blk,X,Z,Xchol,Zchol,v);

   objadd = zeros(1,2); tmp = zeros(1,2);
   for p = 1:size(blk,1)
      pblk = blk(p,:);
      vp = v{p};
      idx = find(vp > 0);
      if ~isempty(idx)
         vpsub = vp(idx);
         if size(vpsub,1) < size(vpsub,2); vpsub = vpsub'; end
         if strcmp(pblk{1},'s')
            ss = [0, cumsum(pblk{2})];
            logdetX = 2*log(diag(Xchol{p}));
            logdetZ = 2*log(diag(Zchol{p}));
            logdetXsub = zeros(length(idx),1);
            logdetZsub = zeros(length(idx),1);
            for k = 1:length(idx)
               idxtmp = [ss(idx(k))+1:ss(idx(k)+1)];
               logdetXsub(k) = sum(logdetX(idxtmp));
               logdetZsub(k) = sum(logdetZ(idxtmp));
            end
            tmp(1) = -sum(vpsub.*logdetXsub);
            tmp(2) = sum(vpsub.*logdetZsub + (pblk{2}(idx)').*vpsub.*(1-log(vpsub)));
         elseif strcmp(pblk{1},'q')
            gamX = sqrt(qops(pblk,X{p},X{p},2));
            gamZ = sqrt(qops(pblk,Z{p},Z{p},2));
            tmp(1) = -sum(vpsub.*log(gamX(idx)));
            tmp(2) = sum(vpsub.*log(gamZ(idx)) + vpsub);
         elseif strcmp(pblk{1},'l')
            logX = log(X{p}); logZ = log(Z{p});
            tmp(1) = -sum(vpsub.*logX(idx));
            tmp(2) = sum(vpsub.*logZ(idx) + vpsub.*(1-log(vpsub)));
         end
         objadd = objadd + tmp;
      end
   end
%%********************************************************************

author: yu-jiang

repo: radpbook

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


text:
function [uwav]=u_wav(x,a_uwav,d_uwav,t_uwav,li)
l=li;
a=a_uwav
x=x-t_uwav;
b=(2*l)/d_uwav;
n=100;
u1=1/l
u2=0
for i = 1:n
    harm4=(((sin((pi/(2*b))*(b-(2*i))))/(b-(2*i))+(sin((pi/(2*b))*(b+(2*i))))/(b+(2*i)))*(2/pi))*cos((i*pi*x)/l);
    u2=u2+harm4;
end
uwav1=u1+u2;
uwav=a*uwav1;


text:
function [h,t,dat,d,m1,m2,sterr] = tor_fill_steplot(dat,color,varargin)
% :Usage:
% ::
%
%    [h,t] = tor_fill_steplot(dat,color,[robust flag],[p-thresh],[x vector],[covs no interest])
%
% Plots a mean vector (mean of each column of dat)
% surrounded by a fill with standard err bars
%
% If dat has 3 dimensions, then
% the diff between dat(:,:,1) and dat(:,:,2) is
% used as the difference for computing standard err
% (as in repeated measures)
%
% if behavior is entered as optional argument, removes it before plotting
% lines.  Also returns adjusted output in d, dat
%
% Optional: robust flag (1/0), robust IRLS
%
% :Examples:
% ::
%
%    tor_fig;
%    tor_fill_steplot(dat,{'b' 'r'},0,.05,secs);

dorobust = 0; pthresh = 0; x = 1:size(dat,2);
t = []; m1 = []; m2 = []; X = [];

if length(varargin) > 0, dorobust = varargin{1};,end
if length(varargin) > 1, pthresh = varargin{2};,end
if length(varargin) > 2, x = varargin{3};,end   % small x, time
if length(varargin) > 3, X = varargin{4};,end   % big X, model matrix of covs

if isempty(x), x = 1:size(dat,2);, end

if length(size(dat)) > 2

    % remove covariates of no interest, if any
    % center covs and fit without intercept to preserve mean in data
    if ~isempty(X)
        disp('Adjusting for covariates.');
        X = scale(X,1); % center
        W = X * pinv(X);
        y = squeeze(dat(:,:,1));
        dat(:,:,1) = y - W * y;

        y = squeeze(dat(:,:,2));
        dat(:,:,2) = y - W * y;
    end

    d = dat(:,:,1) - dat(:,:,2);

    if dorobust
        [m1] = robust_mean(dat(:,:,1));
        [m2] = robust_mean(dat(:,:,2));
        [dummy,t,p,sterr] = robust_mean(d);
    else

        m1 = nanmean(dat(:,:,1));
        m2 = nanmean(dat(:,:,2));
        md = nanmean(d);
        sterr = ste(d);
        t = md ./ sterr;
        [h,p,ci,stat] = ttest(d);
    end

    if ~iscell(color), error('For two groups, color should be cell, e.g.,  {''r'' ''b''})');,end

    hold on;
    h(1) = plot(x,m1,color{1},'LineWidth',2);
    h(2) = plot(x,m2,color{2},'LineWidth',2);

    drawnow

    fill_around_line(m1,sterr,color{1},x);
    fill_around_line(m2,sterr,color{2},x);

    drawnow

    if pthresh
        % significance markers at top of plot.
        yval = max([m1 m2]) + .04 * max([m1 m2]);
        k = 1;
        df = size(d,1) - k;
        %tthr = tinv(1 - pthresh,df);
        tsig = (p < pthresh) .* t;

        wh = yval * (tsig > 0);%double((t > tthr)); wh(wh==0) = NaN; wh(wh>0) = yval;
        %plot(x,wh,color{1},'LineWidth',3);
        for i = 1:length(wh)
            if wh(i)
                text(x(i),wh(i),'*','Color',color{1}(1),'FontSize',24);
            end
        end

        %wh = yval(find(tsig < 0)); % double((t < -tthr)); wh(wh==0) = NaN; wh(wh>0) = yval;
        %plot(x,wh,color{2},'LineWidth',3);

        wh = yval * (tsig < 0);%double((t > tthr)); wh(wh==0) = NaN; wh(wh>0) = yval;
        %plot(x,wh,color{1},'LineWidth',3);
        for i = 1:length(wh)
            if wh(i)
                text(x(i),wh(i),'*','Color',color{2}(1),'FontSize',24);
            end
        end

        text(x(5),yval+.04 * max([m1 m2]),'Significant','FontSize',16);
    end

    drawnow

else

    if dorobust
        [m1,t,p,sterr] = robust_mean(dat(:,:,1));
    else
        m1 = nanmean(dat(:,:,1), 1);
        sterr = ste(dat);
    end

    if ~iscell(color), tmp = color; color = []; color{1} = tmp; end

    hold on;
    if ischar(color{1})
        h = plot(x,m1,color{1},'LineWidth',2);
    else
        h = plot(x,m1,'o-', 'Color', color{1},'LineWidth',2);
    end

    fill_around_line(m1,sterr,color{1},x);

end

return


text:
function [f, tElapsed] = Complex_DIIVINE_feature(im)
% 82 features are extracted:
% Magnitude 1:16
% Relative Mag: 17:28
%Phase: 29:52
%CW-ssim: 53:82
% Magnitude: 1st scale -> 14 features (MGGD); 2nd scale -> 2 features (MGGD);
%        finest scale -> 12 features (GGD, relative magnitude, combined orientations)
% Fhase: 24 features (first 2 scales, relative phase, two orientations)
% Scale Correlation: 30 features (1&2, 2&3, 1&3, hpr & 1,2)

num_scale = 3;
num_ori = 6;
if(size(im,3)~=1)
    im = (double(rgb2gray(im)));
else
    im = double(im);
end

% complex steerable pyramid decomposition
[pyr, pind] = buildSCFpyr(im, num_scale, num_ori-1);
% f = [];
tStart_1 = tic;
[band_dnt ind] = complex_divisive_normalized(pyr,pind,num_scale,num_ori,1,1,3,3);

% Magnitude features
shifts = [0 1; 1 0; 1 1; -1 1];
para_diff = [];
para = [];
b_s = [];

for i = 1:num_ori
    b = abs(band_dnt{i});
    b_s = [b_s; b(:)];
    [alpha, beta] = MGGD_ParaEstimate(b(:));
    para = [para; log(alpha) beta];

    b_diff = b + circshift(b, shifts(3,:)) - circshift(b, shifts(1,:)) - circshift(b, shifts(2,:));
    [sigma, gam] = gaussian_para_esti(b_diff(:));
    para_diff = [para_diff; sigma gam];
end
[alpha_s, beta_s] = MGGD_ParaEstimate(b_s);
para = [para; log(alpha_s) beta_s];

for s_n = 2:num_scale-1
    b_s = [];
    for ori_n = 1:num_ori
        b_ori = abs(band_dnt{(s_n-1)*num_ori+ori_n});
        b_s = [b_s; b_ori(:)];
    end
    [alpha_s, beta_s] = MGGD_ParaEstimate(b_s);
    para = [para; log(alpha_s) beta_s];
end
tElapsed_1 = toc(tStart_1);

f = [para(:); para_diff(:)];

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% horizontal and vertical relative phase
tStart_2 = tic;
ph_rela = RelativePhase_HV(pyr, pind, num_ori);
cell_len = length(ph_rela);
wrap_cauchy = [];
for k = 1:cell_len
    re_phase = ph_rela{k};
    [mu, row, loop_num] = WrapCauchyEstimate(re_phase);
    wrap_cauchy = [wrap_cauchy; row];
end
tElapsed_2 = toc(tStart_2);
f = [f; wrap_cauchy];

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% CW-SSIM
tStart_3 = tic;
K = 0;
cw_ssim_pyr = Scale_CW_SSIM(pyr, pind, num_scale, num_ori, K);
cw_ssim = cw_ssim_pyr';
cw_ssim = cw_ssim(:);
tElapsed_3 = toc(tStart_3);

f = [f; cw_ssim];

tElapsed = [tElapsed_1 tElapsed_2 tElapsed_3];

% sub-function
function phase_relative = RelativePhase_HV(pyr, pind, NumOri)
m = size(pind,1);
shifts = [0 1; 1 0];
for bnum = 2:m-1*NumOri-1
    band = pyrBand(pyr,pind,bnum);
    phase = angle(band);
    for j = 1:size(shifts, 1)
        phase_shift_hori = shift(phase,shifts(j,:));
        phase_rela_hori = phase_shift_hori - phase;
        phase_rela_hori = phase_rela_hori(:,2:end);
        ind_1_hori = phase_rela_hori < -pi;
        ind_2_hori = phase_rela_hori > pi;
        ind_3_hori = (phase_rela_hori >= -pi)&(phase_rela_hori <= pi);
        ph_rela_hori = ind_1_hori .*(phase_rela_hori+2*pi) + ind_2_hori .*(phase_rela_hori-2*pi) + ind_3_hori .*phase_rela_hori;
        phase_relative(bnum - 1, j) = {[ph_rela_hori]};
    end
end
phase_relative = phase_relative(:);

function cw_ssim_pyr = Scale_CW_SSIM(pyr, pind, num_scale, num_ori, K)
m = size(pind,1);
for bnum = 1:m
    band = pyrBand(pyr,pind,bnum);
    suband(bnum) = {[band]};
end
for scale = 1:num_scale-1
    s1_ind = m - num_ori*scale;
    s2_ind = m - num_ori*(scale+1);

    for ori = 1:num_ori
        band1 = suband{1,(s1_ind + ori - 1)};
        band2 = suband{1,(s2_ind + ori - 1)};
        cw_ssim_pyr(scale,ori) = cwssim_index_new(imresize(band1,size(band2)), band2, K);
    end
end

scale = 1;
s1_ind = m - num_ori*scale;
s2_ind = m - num_ori*(scale+2);

for ori = 1:num_ori
    band1 = suband{1,(s1_ind + ori - 1)};
    band2 = suband{1,(s2_ind + ori - 1)};
    cw_ssim_ori(scale, ori) = cwssim_index_new(imresize(band1,size(band2)), band2, K);
end

st = 2; % only calculate HPR band and 1:num_scale+1-st scale subband
hp_band = suband{1,1};
for scale = st:num_scale
    s_ind = m - num_ori*scale;
    for ori = 1:num_ori
        bp_band =  suband{1,(s_ind + ori - 1)};
        cw_ssim_pyr(scale+num_scale-st,ori) = cwssim_index_new(imresize(bp_band,size(hp_band)), hp_band, K);
    end
end
cw_ssim_pyr = [cw_ssim_pyr; cw_ssim_ori];

function [mu, row, k] = WrapCauchyEstimate(AngleMatrix)
u1 = 0.3; u2 = 0.3;
e = .00001;
k = 0;
while (k < 1000)
    mu1 = u1;
    mu2 = u2;
    w = 1./(1 - mu1*cos(AngleMatrix) - mu2*sin(AngleMatrix));
    num1 = w.*cos(AngleMatrix);
    num2 = w.*sin(AngleMatrix);
    u1 = sum(num1(:))/sum(w(:));
    u2 = sum(num2(:))/sum(w(:));
    k = k + 1;
    if (abs(u1-mu1)< e)&&(abs(u2-mu2)< e)
        break;
    end
end
if k == 1000
    error('Cauchy data do not converge');
   %mu = 0; row = 0; k =0;
end
mu = atan(u2/u1);
row = (1 - sqrt(1 - u1^2 - u2^2))/sqrt(u1^2 + u2^2);

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL


text:
%
% In this file, we mex the C files used in this package.

clear, clc;
current_path=cd;

%% Output infor
%%
fprintf('\n ----------------------------------------------------------------------------');
fprintf('\n The program is mexing the C files. Please wait...');
fprintf('\n If you have problem with mex, you can refer to the help of Matlab.');
fprintf('\n If you cannot solve the problem, please contact with Jun Liu (j.liu@asu.edu)\n\n');

%% currently, this package uses the following C files
%%            (in the folder /SLEP/CFiles)

% files in the folder q1
cd([current_path '/SLEP/CFiles/q1']);
mex epp.c;
mex ep1R.c;
mex ep21d.c;
mex ep21R.c;
mex eplb.c;
mex eppMatrix.c;
mex eppVector.c;
mex eppVectorR.c;
mex epsgLasso.c;

% file in the folder flsa
cd([current_path '/SLEP/CFiles/flsa']);
mex flsa.c;

% file in the folder SpInvCoVa
cd([current_path '/SLEP/CFiles/SpInvCoVa']);
mex invCov.c;

% files in the folder tree
cd([current_path '/SLEP/CFiles/tree']);
mex altra.c;
mex altra_mt.c;
mex computeLambda2Max.c;
mex treeNorm.c;
mex findLambdaMax.c;
mex findLambdaMax_mt.c;

mex general_altra.c;
mex general_altra_mt.c;
mex general_treeNorm.c;
mex general_findLambdaMax.c;
mex general_findLambdaMax_mt.c;

% files in the folder order
cd([current_path '/SLEP/CFiles/order']);
mex orderTree.c;
mex orderTree_without_nonnegative.c;
mex orderTreeBinary.c;
mex orderTreeDepth1.c;
mex sequence_bottomup.c;
mex sequence_topdown.c;

%% Output infor
%%
fprintf('\n\n The C files in the folder CFiles have been successfully mexed.');
fprintf('\n\n You can now use the functions in the folder SLEP.');
fprintf('\n You are suggested to read the manual for better using the codes.');
fprintf('\n You are also suggested to run the examples in the folder Examples for these functions.');
fprintf('\n\n These codes are being developed by Jun Liu and Jieping Ye at Arizona State University.');
fprintf('\n If there is any problem, please contact with Jun Liu and Jieping Ye ({j.liu,jieping.ye}@asu.edu).');
fprintf('\n\n Thanks!');
fprintf('\n ----------------------------------------------------------------------------\n');

cd(current_path);

author: zzstefan

repo: BrainNetClass


text:
function [nlZ, dnlZ, HnlZ] = mcs_gp_optimizer(gpstruct, theta)
%MCS_GP_OPTIMIZER Wrapper function for GP optimization

fixed = gpstruct.fixed;
hyp = gpstruct.hyp;
inf = gpstruct.inf;
mean = gpstruct.mean;
cov = gpstruct.cov;
lik = gpstruct.lik;
x = gpstruct.x;
y = gpstruct.y;

newtheta = unwrap2vec(hyp);
newtheta(~fixed) = theta;
theta = rewrap(hyp, newtheta);

f = @() (feval(inf{:}, theta, mean, cov, lik, x, y));

if nargout <= 1
    [~, nlZ] = f();
    return;

elseif nargout == 2
    [~, nlZ, dnlZ] = f();

elseif nargout > 2
    [~, nlZ, dnlZ, ~, ~, HnlZ] = f();

    HnlZ = HnlZ.value;
    HnlZ(fixed,:) = [];
    HnlZ(:,fixed) = [];
end

dnlZ = unwrap2vec(dnlZ);
dnlZ(fixed) = [];

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function marginal = marginal_family(engine, i, t)
% MARGINAL_FAMILY Compute the marginal on the specified family (ff)
% marginal = marginal_family(engine, i, t)

if engine.filter
  error('can''t currently use marginal_family when filtering with ff');
end

if nargin < 3, t = 1; end

% The method is similar to the following HMM equation:
% xi(i,j,t) = normalise( alpha(i,t) * transmat(i,j) * obsmat(j,t+1) * beta(j,t+1) )
% where xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | y(1:T))

bnet = bnet_from_engine(engine);
ss = length(bnet.intra);

if myismember(i, engine.onodes)
  ps = parents(bnet.dag, i);
  p = ps(1);
  marginal = pot_to_marginal(engine.marginals{ps(1),t});
  fam = ([ps i]) + (t-1)*ss;
elseif t==1
  marginal = pot_to_marginal(engine.marginals{i,t});
  fam = i + (t-1)*ss;
else
  pot = engine.CPDpot{i,t};
  c = engine.obschild(i);
  if c>0
    pot = multiply_by_pot(pot, engine.CPDpot{c,t});
  end
  pot = multiply_by_pot(pot, engine.back{i,t});
  ps = parents(bnet.dag, i+ss);
  for p=ps(:)'
    pot = multiply_by_pot(pot, engine.fwd{p,t-1});
  end
  marginal = pot_to_marginal(normalize_pot(pot));
  fam = ([ps i+ss]) + (t-2)*ss;
end

% we convert the domain to the unrolled numbering system
% so that update_ess extracts the right evidence.
marginal.domain = fam;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [newSpikeTrain, newCalcium, newLL] = removeSpike(oldSpikeTrain,oldCalcium,oldLL,filters,tau,obsCalcium,timeToRemove,indx,Dt,A)

% Remove a given spike from the existing spike train.

% Inputs:
% oldSpikeTrain:        current spike train
% oldCalcium:           current noiseless calcium trace
% oldLL:                current value of the log-likelihood function
% filters:              exponential rise and decay kernels for calcium transient
% tau:                  continuous time rise and decay time constants
% obsCalcium:           observed fluorescence trace
% timetoRemove:         time of the spike to be removed
% indx:                 place where the spike to be removed is in the existing spike train vector
% Dt:                   time-bin width
% A:                    spike amplitude

% Outputs:
% newSpikeTrain:        new vector of spike times
% newCalcium:           new noiseless calcium trace
% newLL:                new value of the log-likelihood function

% Author: Eftychios A. Pnevmatikakis and Josh Merel

    tau_h = tau(1);
    tau_d = tau(2);

    ef_h = filters{1,1};
    ef_d = filters{1,2};
    ef_nh = filters{2,1};
    ef_nd = filters{2,2};

    newSpikeTrain = oldSpikeTrain;
    newSpikeTrain(indx) = [];

    %use infinite precision to scale the precomputed FIR approximation to the calcium transient
    wk_h = A*exp((timeToRemove - Dt*ceil(timeToRemove/Dt))/tau_h);
    wk_d = A*exp((timeToRemove - Dt*ceil(timeToRemove/Dt))/tau_d);

    %%%%%%%%%%%%%%%%%
    %handle ef_h first
    newCalcium = oldCalcium;
    tmp = 1+ (floor(timeToRemove):min((length(ef_h)+floor(timeToRemove)-1),length(newCalcium)-1));
    wef_h = wk_h*ef_h(1:length(tmp));
    newCalcium(tmp) = newCalcium(tmp) - wef_h;

    relevantResidual = obsCalcium(tmp)-oldCalcium(tmp);
    relevantResidual(isnan(relevantResidual)) = 0;
    newLL = oldLL - ( wk_h^2*ef_nh(length(tmp)) + 2*relevantResidual*wef_h(:));
    oldCalcium = newCalcium;
    oldLL = newLL;
    %%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%
    %handle ef_d next
    newCalcium = oldCalcium;
    tmp = 1+ (floor(timeToRemove):min((length(ef_d)+floor(timeToRemove)-1),length(newCalcium)-1));
    wef_d = wk_d*ef_d(1:length(tmp));
    newCalcium(tmp) = newCalcium(tmp) - wef_d;

    relevantResidual = obsCalcium(tmp)-oldCalcium(tmp);
    relevantResidual(isnan(relevantResidual)) = 0;
    newLL = oldLL - ( wk_d^2*ef_nd(length(tmp)) + 2*relevantResidual*wef_d(:));
    %%%%%%%%%%%%%%%%

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [Xhat] = aprxMAPGMM(Y,patchSize,noiseSD,imsize,GS,excludeList,SigmaNoise)
% approximate GMM MAP estimation - a single iteration of the "hard version"
% EM MAP procedure (see paper for a reference)
%
% Inputs:
%   Y - the noisy patches (in columns)
%   noiseSD - noise standard deviation
%   imsize - size of the original image (not used in this case, but may be
%   used for non local priors)
%   GS - the gaussian mixture model structure
%   excludeList - used only for inpainting, misleading name - it's a list
%   of patch indices to use for estimation, the rest are just ignored
%   SigmaNoise - if the noise is non-white, this is the noise covariance
%   matrix
%
% Outputs:
%   Xhat - the restore patches

% handle exclusion list - used for inpainting
if ~exist('excludeList','var')
    excludeList = [];
end

% Supports general noise covariance matrices
if (~exist('SigmaNoise','var'))
    SigmaNoise = noiseSD^2*eye(patchSize^2);
end

if ~isempty(excludeList)
    T = Y;
    Y = Y(:,excludeList);
end

% remove DC component
meanY = mean(Y);
Y = bsxfun(@minus,Y,meanY);

% calculate assignment probabilities for each mixture component for all
% patches
GS2 = GS;
PYZ = zeros(GS.nmodels,size(Y,2));
for i=1:GS.nmodels
    GS2.covs(:,:,i) = GS.covs(:,:,i) + SigmaNoise;
    PYZ(i,:) = log(GS.mixweights(i)) + loggausspdf2(Y,GS2.covs(:,:,i));
end

% find the most likely component for each patch
[~,ks] = max(PYZ);

% and now perform weiner filtering
Xhat = zeros(size(Y));
for i=1:GS.nmodels
    inds = find(ks==i);
    Xhat(:,inds) = ((GS.covs(:,:,i)+SigmaNoise)\(GS.covs(:,:,i)*Y(:,inds) + SigmaNoise*repmat(GS.means(:,i),1,length(inds))));
end

% handle exclusion list stuff (inpainting only)
if ~isempty(excludeList)
    tt = T;
    tt(:,excludeList) = bsxfun(@plus,Xhat,meanY);
    Xhat = tt;
else
    Xhat = bsxfun(@plus,Xhat,meanY);
end

author: lbasek

repo: image-denoising-benchmark


text:
classdef SimplAllTestExplicitVsImplicit < handle

     properties (Access = protected)
        rho
        matInterpExplicit
        matInterpImplicit
        materialPropertiesSettings
        nElem
        nGaus
        dim
     end

     methods (Access = public)

        function obj = SimplAllTestExplicitVsImplicit(dim)
            obj.dim = dim;
            obj.init();
            obj.createProperties();
            obj.computeExplicitMaterialInterpolation();
            obj.computeImplicitMaterialInterpolation();
        end

        function error = computeError(obj)
            mE = obj.matInterpExplicit;
            mI = obj.matInterpImplicit;
            err(1) = norm(mE.mu(:) - mI.mu(:));
            err(2) = norm(mE.kappa(:) - mI.kappa(:));
            err(3) = norm(mE.dmu(:) - mI.dmu(:));
            err(4) = norm(mE.dkappa(:) - mI.dkappa(:));
            error = max(err);
        end
     end

     methods (Access = protected)

        function init(obj)
            obj.nElem = 800;
            obj.nGaus = 4;
            obj.rho = rand(obj.nElem,obj.nGaus);
        end

        function createProperties(obj)
            sC = SettingsConstitutiveProperties();
            sC.E_plus   = rand(1);
            sC.E_minus  = rand(1);
            sC.nu_minus = rand(1);
            sC.nu_plus  = rand(1);
            obj.materialPropertiesSettings = sC;
        end

        function computeExplicitMaterialInterpolation(obj)
            type = 'EXPLICIT';
            mI = obj.computeMaterialInterpolation(type);
            obj.matInterpExplicit = mI;
        end

        function computeImplicitMaterialInterpolation(obj)
            type = 'IMPLICIT';
            mI = obj.computeMaterialInterpolation(type);
            obj.matInterpImplicit = mI;
        end

        function p = computeMaterialInterpolation(obj,type)
            sC = obj.materialPropertiesSettings;
            s = SettingsInterpolation();
            s.constitutiveProperties = sC;
            s.nElem = obj.nElem;
            s.simpAllType = type;
            s.dim = obj.dim;
            mI = MaterialInterpolation.create(s);
            p = mI.computeMatProp(obj.rho);
            obj.matInterpExplicit = p;
        end

     end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%  internal function

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)

function save_nii_hdr(hdr, fid)

   if ~isequal(hdr.hk.sizeof_hdr,348),
      error('hdr.hk.sizeof_hdr must be 348.');
   end

   write_header(hdr, fid);

   return;					% save_nii_hdr

%---------------------------------------------------------------------
function write_header(hdr, fid)

        %  Original header structures
	%  struct dsr				/* dsr = hdr */
	%       {
	%       struct header_key hk;            /*   0 +  40       */
	%       struct image_dimension dime;     /*  40 + 108       */
	%       struct data_history hist;        /* 148 + 200       */
	%       };                               /* total= 348 bytes*/

   header_key(fid, hdr.hk);
   image_dimension(fid, hdr.dime);
   data_history(fid, hdr.hist);

   %  check the file size is 348 bytes
   %
   fbytes = ftell(fid);

   if ~isequal(fbytes,348),
      msg = sprintf('Header size is not 348 bytes.');
      warning(msg);
   end

   return;					% write_header

%---------------------------------------------------------------------
function header_key(fid, hk)

   fseek(fid,0,'bof');

	%  Original header structures
	%  struct header_key                      /* header key      */
	%       {                                /* off + size      */
	%       int sizeof_hdr                   /*  0 +  4         */
	%       char data_type[10];              /*  4 + 10         */
	%       char db_name[18];                /* 14 + 18         */
	%       int extents;                     /* 32 +  4         */
	%       short int session_error;         /* 36 +  2         */
	%       char regular;                    /* 38 +  1         */
	%       char hkey_un0;                   /* 39 +  1 */
	%       };                               /* total=40 bytes  */

   fwrite(fid, hk.sizeof_hdr(1),    'int32');	% must be 348.

   % data_type = sprintf('%-10s',hk.data_type);	% ensure it is 10 chars from left
   % fwrite(fid, data_type(1:10), 'uchar');
   pad = zeros(1, 10-length(hk.data_type));
   hk.data_type = [hk.data_type  char(pad)];
   fwrite(fid, hk.data_type(1:10), 'uchar');

   % db_name   = sprintf('%-18s', hk.db_name);	% ensure it is 18 chars from left
   % fwrite(fid, db_name(1:18), 'uchar');
   pad = zeros(1, 18-length(hk.db_name));
   hk.db_name = [hk.db_name  char(pad)];
   fwrite(fid, hk.db_name(1:18), 'uchar');

   fwrite(fid, hk.extents(1),       'int32');
   fwrite(fid, hk.session_error(1), 'int16');
   fwrite(fid, hk.regular(1),       'uchar');

   fwrite(fid, hk.hkey_un0(1),    'uchar');

   return;					% header_key

%---------------------------------------------------------------------
function image_dimension(fid, dime)

	%struct image_dimension
	%       {                                /* off + size      */
	%       short int dim[8];                /* 0 + 16          */
	%       char vox_units[4];               /* 16 + 4          */
	%       char cal_units[8];               /* 20 + 8          */
	%       short int unused1;               /* 28 + 2          */
	%       short int datatype;              /* 30 + 2          */
	%       short int bitpix;                /* 32 + 2          */
	%       short int dim_un0;               /* 34 + 2          */
	%       float pixdim[8];                 /* 36 + 32         */
	%			/*
	%				pixdim[] specifies the voxel dimensions:
	%				pixdim[1] - voxel width
	%				pixdim[2] - voxel height
	%				pixdim[3] - interslice distance
	%					..etc
	%			*/
	%       float vox_offset;                /* 68 + 4          */
	%       float roi_scale;                 /* 72 + 4          */
	%       float funused1;                  /* 76 + 4          */
	%       float funused2;                  /* 80 + 4          */
	%       float cal_max;                   /* 84 + 4          */
	%       float cal_min;                   /* 88 + 4          */
	%       int compressed;                  /* 92 + 4          */
	%       int verified;                    /* 96 + 4          */
	%       int glmax;                       /* 100 + 4         */
	%       int glmin;                       /* 104 + 4         */
	%       };                               /* total=108 bytes */

   fwrite(fid, dime.dim(1:8),      'int16');

   pad = zeros(1, 4-length(dime.vox_units));
   dime.vox_units = [dime.vox_units  char(pad)];
   fwrite(fid, dime.vox_units(1:4),  'uchar');

   pad = zeros(1, 8-length(dime.cal_units));
   dime.cal_units = [dime.cal_units  char(pad)];
   fwrite(fid, dime.cal_units(1:8),  'uchar');

   fwrite(fid, dime.unused1(1),    'int16');
   fwrite(fid, dime.datatype(1),   'int16');
   fwrite(fid, dime.bitpix(1),     'int16');
   fwrite(fid, dime.dim_un0(1),    'int16');
   fwrite(fid, dime.pixdim(1:8),   'float32');
   fwrite(fid, dime.vox_offset(1), 'float32');
   fwrite(fid, dime.roi_scale(1),  'float32');
   fwrite(fid, dime.funused1(1),   'float32');
   fwrite(fid, dime.funused2(1),   'float32');
   fwrite(fid, dime.cal_max(1),    'float32');
   fwrite(fid, dime.cal_min(1),    'float32');
   fwrite(fid, dime.compressed(1), 'int32');
   fwrite(fid, dime.verified(1),   'int32');
   fwrite(fid, dime.glmax(1),      'int32');
   fwrite(fid, dime.glmin(1),      'int32');

   return;					% image_dimension

%---------------------------------------------------------------------
function data_history(fid, hist)

	% Original header structures - ANALYZE 7.5
	%struct data_history
	%       {                                /* off + size      */
	%       char descrip[80];                /* 0 + 80          */
	%       char aux_file[24];               /* 80 + 24         */
	%       char orient;                     /* 104 + 1         */
	%       char originator[10];             /* 105 + 10        */
	%       char generated[10];              /* 115 + 10        */
	%       char scannum[10];                /* 125 + 10        */
	%       char patient_id[10];             /* 135 + 10        */
	%       char exp_date[10];               /* 145 + 10        */
	%       char exp_time[10];               /* 155 + 10        */
	%       char hist_un0[3];                /* 165 + 3         */
	%       int views                        /* 168 + 4         */
	%       int vols_added;                  /* 172 + 4         */
	%       int start_field;                 /* 176 + 4         */
	%       int field_skip;                  /* 180 + 4         */
	%       int omax;                        /* 184 + 4         */
	%       int omin;                        /* 188 + 4         */
	%       int smax;                        /* 192 + 4         */
	%       int smin;                        /* 196 + 4         */
	%       };                               /* total=200 bytes */

   % descrip     = sprintf('%-80s', hist.descrip);     % 80 chars from left
   % fwrite(fid, descrip(1:80),    'uchar');
   pad = zeros(1, 80-length(hist.descrip));
   hist.descrip = [hist.descrip  char(pad)];
   fwrite(fid, hist.descrip(1:80), 'uchar');

   % aux_file    = sprintf('%-24s', hist.aux_file);    % 24 chars from left
   % fwrite(fid, aux_file(1:24),   'uchar');
   pad = zeros(1, 24-length(hist.aux_file));
   hist.aux_file = [hist.aux_file  char(pad)];
   fwrite(fid, hist.aux_file(1:24), 'uchar');

   fwrite(fid, hist.orient(1),      'uchar');
   fwrite(fid, hist.originator(1:5), 'int16');

   pad = zeros(1, 10-length(hist.generated));
   hist.generated = [hist.generated  char(pad)];
   fwrite(fid, hist.generated(1:10),  'uchar');

   pad = zeros(1, 10-length(hist.scannum));
   hist.scannum = [hist.scannum  char(pad)];
   fwrite(fid, hist.scannum(1:10),  'uchar');

   pad = zeros(1, 10-length(hist.patient_id));
   hist.patient_id = [hist.patient_id  char(pad)];
   fwrite(fid, hist.patient_id(1:10),  'uchar');

   pad = zeros(1, 10-length(hist.exp_date));
   hist.exp_date = [hist.exp_date  char(pad)];
   fwrite(fid, hist.exp_date(1:10),  'uchar');

   pad = zeros(1, 10-length(hist.exp_time));
   hist.exp_time = [hist.exp_time  char(pad)];
   fwrite(fid, hist.exp_time(1:10),  'uchar');

   pad = zeros(1, 3-length(hist.hist_un0));
   hist.hist_un0 = [hist.hist_un0  char(pad)];
   fwrite(fid, hist.hist_un0(1:3),  'uchar');

   fwrite(fid, hist.views(1),      'int32');
   fwrite(fid, hist.vols_added(1), 'int32');
   fwrite(fid, hist.start_field(1),'int32');
   fwrite(fid, hist.field_skip(1), 'int32');
   fwrite(fid, hist.omax(1),       'int32');
   fwrite(fid, hist.omin(1),       'int32');
   fwrite(fid, hist.smax(1),       'int32');
   fwrite(fid, hist.smin(1),       'int32');

   return;					% data_history

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function angles3d(varargin)
%ANGLES3D Conventions for manipulating angles in 3D.
%
%   The library uses both radians and degrees angles;
%   Results of angle computation between shapes usually returns angles in
%   radians.
%   Representation of 3D shapes use angles in degrees (easier to manipulate
%   and to save).
%
%   Contrary to the plane, there are no oriented angles in 3D. Angles
%   between lines or between planes are comprised between 0 and PI.
%
%   Spherical angles
%   Spherical angles are defined by 2 angles:
%   * THETA, the colatitude, representing angle with Oz axis (between 0 and
%       PI)
%   * PHI, the azimut, representing angle with Ox axis of horizontal
%       projection of the direction (between 0 and 2*PI)
%
%   Spherical coordinates can be represented by THETA, PHI, and the
%   distance RHO to the origin.
%
%   Euler angles
%   Some functions for creating rotations use Euler angles. They follow the
%   ZYX convention in the global reference system, that is eqivalent to the
%   XYZ convention ine a local reference system.
%   Euler angles are given by a triplet of angles [PHI THETA PSI] that
%   represents the succession of 3 rotations:
%   * rotation around X by angle PSI    ("roll")
%   * rotation around Y by angle THETA  ("pitch")
%   * rotation around Z by angle PHI    ("yaw")
%
%   In this library, euler angles are given in degrees. The functions that
%   use euler angles use the keyword 'Euler' in their name.
%
%
%   See also
%   cart2sph2, sph2cart2, cart2sph2d, sph2cart2d
%   anglePoints3d, angleSort3d, sphericalAngle, randomAngle3d
%   dihedralAngle, polygon3dNormalAngle, eulerAnglesToRotation3d
%   rotation3dAxisAndAngle, rotation3dToEulerAngles
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2008-10-13,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.


text:
function tc = tc_addMenus(tc,hfig);
% tc = tc_addMenus(tc,[hfig]): add menus to time course UI figure.
%
% tc: the analysis struct used by the UI.
%
% hfig: handle to the figure (default is gcf).
%
% 07/04 ras.
% 09/05 ras: made each top-level menu a sub-function, to
% make 'em easier to get at.
if ~exist('hfig')   hfig = gcf;     end

tc.ui.fig = hfig;

% turn off standard menus for now:
set(hfig,'Menubar','none');

% add each top-level menu individually
tc.ui = tc_plotMenu(tc, tc.ui);
tc.ui = tc_viewMenu(tc, tc.ui);
tc.ui = tc_analMenu(tc, tc.ui);
tc.ui = tc_settingsMenu(tc, tc.ui);
tc.ui = tc_condMenu(tc, tc.ui);
tc.ui = tc_exportMenu(tc, tc.ui);
tc.ui = tc_helpMenu(tc, tc.ui);

% add standard MATLAB menus
tc.ui.spacerMenu = uimenu('Label',' ');
% eval(get(tc.ui.figMenuToggle,'Callback'));

return
% /---------------------------------------------------------------/ %

% /---------------------------------------------------------------/ %
function h = tc_plotMenu(tc,h);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add menus (1): Plotting options for time course
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h.plotMenu = uimenu('ForegroundColor',[0 0 0],'Label','Plot','Separator','on');

% callback for all plotting options:
%   TMP = get(gcf, 'UserData');
%   TMP_H = get(gcbo, 'UserData');
%   set(TMP_H, 'Checked', 'off');
%   set(gcbo, 'Checked', 'on');
%   TMP.plotType = find(TMP_H==gcbo);
%   set(gcf, 'UserData', TMP);
%   timeCourseUI;
%   clear TMP_H TMP;
cb=['TMP_H = get(gcbo, ''UserData''); ' ...
    'TMP = get(gcf, ''UserData''); ' ...
    'set(TMP_H, ''Checked'', ''off''); ' ...
    'set(gcbo, ''Checked'', ''on''); ' ...
    'TMP.plotType = find(TMP_H==gcbo); ' ...
    'set(gcf, ''UserData'', TMP); ' ...
    'timeCourseUI; clear TMP_H TMP; ' ];

h.plotWholeTc = uimenu(h.plotMenu,'Label','Plot Whole Time Course','Separator','off',...
    'Accelerator','K','Callback',cb);

h.sparklineWholeTc = uimenu(h.plotMenu, 'Label', 'Plot Whole Time Course (sparkline)',...
                        'Callback',cb);

h.plotMeanTrials = uimenu(h.plotMenu,'Label','Plot Mean Time Course','Separator','off',...
    'Checked','off','Accelerator','M','Callback',cb);

h.meanSubplots = uimenu(h.plotMenu,'Label','Mean Time Course Subplots','Separator','off',...
    'Checked','off','Accelerator','B','Callback',cb);

h.plotAllTrials = uimenu(h.plotMenu,'Label','Plot All Trials','Separator','off',...
    'Accelerator','I','Callback',cb);

h.plotRelAmps = uimenu(h.plotMenu,'Label','Plot Relative Amplitudes','Separator','on',...
    'Accelerator','J','Callback',cb);

h.plotMeanAmps = uimenu(h.plotMenu,'Label','Plot Mean Amplitudes','Separator','off',...
    'Accelerator','U','Callback',cb);

h.plotSummary = uimenu(h.plotMenu,'Label','Mean Time Course Plus Amplitudes','Separator','on',...
    'Checked','on','Accelerator','D','Callback',cb);

h.plotFFT = uimenu(h.plotMenu,'Label','FFT of mean tSeries','Separator','on',...
    'Checked','off','Accelerator','F','Callback',cb);

h.plotGlm = uimenu(h.plotMenu,'Label','Visualize GLM Results','Separator','on',...
    'Checked','off','Accelerator','G','Callback',cb);

h.plotCorAnal = uimenu(h.plotMenu,'Label','Visualize Cor Anal','Separator','off',...
    'Checked','off','Accelerator','T','Callback',cb);

% this code allows the UI to exclusively check the selected
% plot option, and deselect the others:
handlesArray = [h.plotWholeTc h.plotAllTrials h.plotMeanTrials ...
                h.plotRelAmps h.plotMeanAmps h.meanSubplots ...
                h.plotSummary h.plotFFT h.plotGlm h.plotCorAnal ...
                h.sparklineWholeTc];
set(h.plotWholeTc,'UserData',handlesArray);
set(h.plotAllTrials,'UserData',handlesArray);
set(h.plotMeanTrials,'UserData',handlesArray);
set(h.plotRelAmps,'UserData',handlesArray);
set(h.plotMeanAmps,'UserData',handlesArray);
set(h.meanSubplots,'UserData',handlesArray);
set(h.plotSummary,'UserData',handlesArray);
set(h.plotFFT,'UserData',handlesArray);
set(h.plotGlm,'UserData',handlesArray);
set(h.plotCorAnal,'UserData',handlesArray);
set(h.sparklineWholeTc,'UserData',handlesArray);

return
% /---------------------------------------------------------------/ %

% /---------------------------------------------------------------/ %
function h = tc_viewMenu(tc,h);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add menus (2): View options for time course
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h.viewMenu = uimenu('ForegroundColor',[0 0 0],'Label','View','Separator','on');

% spawn new UI callback:
% tc_openFig(get(gcf,'UserData'));
h.newTCUI = uimenu(h.viewMenu,'Label','New UI Window','Separator','off',...
    'Callback',' tc_openFig(get(gcf,''UserData'')); timeCourseUI;');

% view/edit parfiles callback:
% tc = get(gcf, 'UserData');
% for p = tc.trials.parfiles, edit(p{1}); end
cb = 'tc = get(gcf, ''UserData''); ';
cb = [cb 'for p = tc.trials.parfiles, edit(p{1}); end '];
h.editParfiles = uimenu(h.viewMenu, 'Label', 'View/Edit Parfiles', ...
    'Separator', 'on', 'Callback', cb);

% mark individual trials in All Trials plot:
%   umtoggle(gcbo);
%   TMP = get(gcf,'UserData');
%   TMP.params.markEachTrial = isequal(get(gcbo,'Checked'),'on');
%   set(gcf,'UserData',TMP);
%   timeCourseUI;
%   clear TMP
cb = ['umtoggle(gcbo); '...
       'TMP = get(gcf,''UserData''); '...
       'TMP.params.markEachTrial = isequal(get(gcbo,''Checked''),''on''); '...
       'set(gcf, ''UserData'', TMP); '...
       'timeCourseUI; clear TMP;'];
h.markEachTrial = uimenu(h.viewMenu,'Label','Color Single Trials (All Trials Plot)', ....
                'Separator','on', 'Checked', 'on', 'Callback', cb);

% show grid:
%   umtoggle(gcbo);
%   TMP = get(gcf,'UserData');
%   TMP.params.grid = isequal(get(gcbo,'Checked'),'on');
%   set(gcf,'UserData',TMP);
%   timeCourseUI;
%   clear TMP
cb = ['umtoggle(gcbo); '...
       'TMP = get(gcf,''UserData''); '...
       'TMP.params.grid = isequal(get(gcbo,''Checked''),''on''); '...
       'set(gcf, ''UserData'', TMP); '...
       'timeCourseUI; clear TMP;'];
h.grid = uimenu(h.viewMenu, 'Label', 'Grid',  'Separator', 'off', ...
                'Checked', 'off', 'Callback', cb);

% show peak and baseline periods callback:
%   umtoggle(gcbo);
%   TMP = get(gcf,'UserData');
%   TMP.params.showPkBsl = isequal(get(gcbo,'Checked'),'on');
%   set(gcf,'UserData',TMP);
%   timeCourseUI;
%   clear TMP
cb = ['umtoggle(gcbo); '...
       'TMP = get(gcf,''UserData''); '...
       'TMP.params.showPkBsl = isequal(get(gcbo,''Checked''),''on''); '...
       'set(gcf, ''UserData'', TMP); '...
       'timeCourseUI; clear TMP;'];
h.showPkBsl = uimenu(h.viewMenu,'Label','Show Peak/Baseline Periods', ....
                'Separator','off', 'Checked', 'on', 'Callback', cb);

% legend callback:
% TMP = get(gcf, 'UserData');
% TMP.params.legend = mrvPanelToggle(TMP.ui.legend, gcbo);
% set(gcf, 'UserData', TMP);
% clear TMP
cb = ['TMP = get(gcf, ''UserData''); ' ...
      'TMP.params.legend = mrvPanelToggle(TMP.ui.legend, gcbo); ' ...
      'set(gcf, ''UserData'', TMP); ' ...
      'clear TMP; '];
h.showLegend = uimenu(h.viewMenu, 'Label', 'Show Legend', ...
    'Separator', 'off', 'Checked', 'on', 'Accelerator', 'L', ...
    'Callback', cb);

h.figMenuToggle = addFigMenuToggle(h.viewMenu);

return
% /---------------------------------------------------------------/ %

% /---------------------------------------------------------------/ %
function h = tc_analMenu(tc,h);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add menus (3): Analysis options for time course
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h.analMenu = uimenu('ForegroundColor',[0 0 0],'Label','Analysis','Separator','on');

% calc SNR (peak v bsl)
cb = 'tc = get(gcf,''UserData''); ';
cb = [cb 'title(sprintf(''SNR Peak Vs Baseline: %2.2f (%2.2f dB)'',tc.SNR,tc.SNRdb),''FontSize'',14);'];
h.calcSNR = uimenu(h.analMenu,'Label','Calculate SNR (peak vs bsl)','Separator','off','Callback',cb);

% blur time course
cb = 'tc_blurTimeCourse;';
h.blurTc = uimenu(h.analMenu,'Label','Blur Time Course','Separator','off','Callback',cb);

% detrend time course
cb = 'tc_detrendTimeCourse;';
h.blurTc = uimenu(h.analMenu,'Label','Detrend Time Course','Separator','off','Callback',cb);

% remove outlier data points
cb = 'tc_removeOutliers;';
h.blurTc = uimenu(h.analMenu,'Label','Remove Outlier Time Points',...
    'Separator','on','Callback',cb);

% apply corAnal
cb = 'tc_applyCorAnal; tc_visualizeCorAnal; ';
h.applyCorAnal = uimenu(h.analMenu, 'Label', 'Apply Cor Anal',...
    'Separator', 'on', 'Callback', cb);

% apply GLM
cb = 'tc_applyGlm; tc_visualizeGlm; ';
h.applyGLM = uimenu(h.analMenu, 'Label', 'Apply GLM',...
    'Separator', 'off', 'Callback', cb);

% Apply statistical contrast callback:
% tc_contrast;
h.contrast = uimenu(h.analMenu,'Label', 'Apply Statistical Contrast',...
    'Separator', 'on', 'Callback', ' tc_contrast;');

% % for a certain design type, with alternating null periods b/w active
% % conditions, allow the user to automatically merge each null condition
% % with the previous trial, (maybe) saving a parfile along the way
%  if all(tc.trials.cond(1:2:end)==0) | all(tc.trials.cond(2:2:end)==0)
%     uimenu(h.analMenu,'Label','Remove Null Conditions...','Separator','on',...
%                      'Callback','tc_mergeNullTrials;');
%  end

% save HRF callback:
% tc_saveHrf;
h.saveHrf = uimenu(h.analMenu,'Label','Save HRF','Separator','on',...
    'Callback',' tc_saveHrf; timeCourseUI;');

% dump data to workspace callback:
% tc_dumpDataToWorkspace;
h.dumpToWs = uimenu(h.analMenu,'Label','Dump Data to Workspace','Separator','on',...
    'Callback','tc_dumpDataToWorkspace;');

return
% /---------------------------------------------------------------/ %

% /---------------------------------------------------------------/ %
function h = tc_settingsMenu(tc,h);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add menus (3): Settings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h.settingsMenu = uimenu('ForegroundColor',[0 0 0],'Label','Settings','Separator','on');

% recompute tc callback:
% tc = get(gcf,'UserData');
% tc_recomputeTc(tc);
cb = sprintf('tc = get(gcf,''UserData'');\n');
cb = [cb sprintf('tc = tc_recomputeTc(tc);')];
h.recompute = uimenu(h.settingsMenu,'Label','Change chopTSeries params...','Separator','off','Callback',cb);

% set zero %signal submenu
h.zeroSubMenu = uimenu(h.settingsMenu,'Label','Set Zero % Value...',...
                        'Separator','on');

% set mean of data to 0% callback:
% tc_setZeroPt(1);
h.setMean2Null = uimenu(h.zeroSubMenu,'Label','Set Mean of All Data as Zero',...
                    'Separator','off','Callback','tc_setZeroPt([],-1);');

% set null to 0% callback:
% tc_setZeroPt(1);
h.setZero2Null = uimenu(h.zeroSubMenu,'Label','Set Null Condition as Zero',...
                    'Separator','off','Callback','tc_setZeroPt([],-2);');

% set other condition as 0% callback:
% val = input('What condition should be set as the 0% signal condition? ');
% tc_setZeroPt(val);
h.setOther2Null = uimenu(h.zeroSubMenu,'Label','Set Other Condition as Zero',...
                    'Separator','off','Callback','tc_setZeroPt([],-3);');

% set null to 0% callback:
% tc_setZeroPt(1);
h.setGLMNull = uimenu(h.zeroSubMenu, 'Label', 'Estimate Null from GLM',...
                    'Separator', 'off', 'Callback', 'tc_setZeroPt([],-4);');

% set axis bounds submenu
h.axisSubmenu = uimenu(h.settingsMenu,'Label','Axis Bounds...','Separator','on');

% manually reset
uimenu(h.axisSubmenu,'Label','Set Manually','Separator','off',...
                     'Callback','tc_setAxisBounds;');
% reset to nice value
cb = 'tc_setAxisBounds([],''auto'');';
uimenu(h.axisSubmenu,'Label','Reset to nice values','Separator','off',...
                     'Callback',cb);

% add also a way to change the assigned colors for each condition
uimenu(h.settingsMenu,'Label','Assign Condition Colors...',...
                     'Separator','off','Callback','tc_assignColors;');

% add also a way to change the assigned names for each condition
uimenu(h.settingsMenu,'Label','Assign Condition Names...','Separator','off',...
                     'Callback','tc_assignNames;');

% allow reordering of condition #s
uimenu(h.settingsMenu,'Label','Reorder Condition #s...','Separator','on',...
                     'Callback','tc_reorderCondNums;');

% group several conditions together
uimenu(h.settingsMenu,'Label','Group Conditions...','Separator','off',...
                     'Callback','tc_groupConditions([], ''dialog'');');

% % allow saving of the parfile info (will soon add color info, not yet)
% uimenu(h.settingsMenu,'Label','Save parfiles...','Separator','on',...
%                      'Callback','tc = get(gcf,''UserData''); tc_saveParfiles(tc);');
return
% /---------------------------------------------------------------/ %

% /---------------------------------------------------------------/ %
function h = tc_condMenu(tc,h);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add menus (4): condition menu w/ conditions toggle buttons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Callback for all menu items:
% umtoggle(gcbo);
% tc_legend(get(gcf, 'UserData'));
% timeCourseUI;
h.condMenu = uimenu('ForegroundColor',[0 0 0],'Label','Conditions','Separator','on');
accelChars = '0123456789-=|';
cb = ['umtoggle(gcbo); tc_legend(get(gcf,''UserData'')); timeCourseUI; '];
for i = 1:length(tc.trials.condNames)
    if i < length(accelChars)
        accel = accelChars(i);
    else
        accel = '';
    end

    if isempty(tc.trials.condNames{i})
        tc.trials.condNames{i} = num2str(i);
    end

    hc(i) = uimenu(h.condMenu, 'Label', tc.trials.condNames{i}, ...
                 'Separator', 'off', 'Checked', 'on', ...
                 'Accelerator', accel, ...
                 'Tag', tc.trials.condNames{i}, ...
                 'UserData', tc.trials.condNums(i), ...
                 'Callback', cb);
end

% unselect the null condition if there is one
if any(tc.trials.condNums==0)
    null = find(tc.trials.condNums==0);
    set(hc(null), 'Checked', 'off');
end
h.condMenuHandles = hc;

return
% /---------------------------------------------------------------/ %

% /---------------------------------------------------------------/ %
function h = tc_exportMenu(tc,h);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add menus (6): export options menu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h.exportMenu = uimenu('Label', 'Export', 'Separator', 'on');

h.exportOptions = uimenu(h.exportMenu, 'Label', 'Export Options', ...
							'Callback', 'exportsetupdlg;');

% dump data to workspace:
% (replicated from analysis menu, will disable it once people get used to
% this new menu)
h.dumpToWs = uimenu(h.exportMenu, 'Label', 'Dump Data to Workspace', ...
					'Callback', 'tc_dumpDataToWorkspace;');

% export parfile info (moved over from analysis menu)
uimenu(h.exportMenu, 'Label', 'Save parfiles as...', ...
		 'Callback', 'tc_saveParfiles(get(gcf,''UserData''));');

return
% /---------------------------------------------------------------/ %

% /---------------------------------------------------------------/ %
function h = tc_helpMenu(tc,h);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add menus (6): help menu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h.helpMenu = uimenu('ForegroundColor',[.6 .2 0],'Label','Help','Separator','on');

% web page callback:
% web http://white.stanford.edu/newlm/index.php/Time_Course_UI
cb = 'web http://white.stanford.edu/newlm/index.php/Time_Course_UI';
h.webHelp = uimenu(h.helpMenu,'Label','Time Course UI page','Separator','off',...
    'Callback',cb);

% web page (external browser) callback:
% web http://white.stanford.edu/newlm/index.php/Time_Course_UI -browser
cb = 'web http://white.stanford.edu/newlm/index.php/Time_Course_UI -browser';
cb = 'web http://white/newlm/index.php/Main_Page -browser';
h.webHelp2 = uimenu(h.helpMenu,'Label','Time Course UI page (external browser)','Separator','off',...
    'Callback',cb);

% mrVista web page callback:
% web web http://white.stanford.edu/newlm/index.php/Main_Page
cb = 'web http://white.stanford.edu/newlm/index.php/Main_Page';
h.webWiki = uimenu(h.helpMenu, 'Label', 'mrVista wiki', ...
    'Separator', 'off', 'Callback', cb);

% mrVista web page (external browser) callback:
% web web http://white.stanford.edu/newlm/index.php/Main_Page
% -browser
h.webWiki2 = uimenu(h.helpMenu, 'Label', 'mrVista wiki (external browser)', ...
    'Separator', 'off', 'Callback', cb);

% identify callback for menu item
cb = 'helpFindCallback; ';
uimenu(h.helpMenu, 'Label', 'Identify Callback for a menu item', ...
        'Separator', 'on', 'Callback', cb);

% add handles to tc struct
fields = fieldnames(h);
for i = 1:length(fields)
    tc.ui.(fields{i}) = h.(fields{i});
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft