text:
function sp = mdp_span(W)

% mdp_span  Returns the span of W
%           sp(W) = max W(s) - min W(s)

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% Redistribution and use in source and binary forms, with or without modification,
% are permitted provided that the following conditions are met:
%    * Redistributions of source code must retain the above copyright notice,
%      this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright notice,
%      this list of conditions and the following disclaimer in the documentation
%      and/or other materials provided with the distribution.
%    * Neither the name of the <ORGANIZATION> nor the names of its contributors
%      may be used to endorse or promote products derived from this software
%      without specific prior written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
% OF THE POSSIBILITY OF SUCH DAMAGE.

sp = max(W) - min(W);


text:
function a = conex3_inverse ( n )

%*****************************************************************************80
%
%% CONEX3_INVERSE returns the inverse of the CONEX3 matrix.
%
%  Example:
%
%    N = 5
%
%     1  0  0  0  0
%     1  1  0  0  0
%     2  1  1  0  0
%     4  2  1  1  0
%    -8 -4 -2 -1 -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alan Cline, RK Rew,
%    A set of counterexamples to three condition number estimators,
%    SIAM Journal on Scientific and Statistical Computing,
%    Volume 4, 1983, pages 602-611.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      if ( i < n )

        if ( j < i )
          a(i,j) = 2.0^(i-j-1);
        elseif ( i == j )
          a(i,j) = 1.0;
        else
          a(i,j) = 0.0;
        end

      elseif ( i == n )

        if ( j < i )
          a(i,j) = - 2.0^(i-j-1);
        else
          a(i,j) = -1.0;
        end

      end

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear, clc;

% This is an example for running the function nnLogisticC
%
%  Problem:
%
%  min  f(x,c) = - weight_i * log (p_i) + 1/2 * rsL2 * ||x||_2^2
%  s.t. \|x\|_1 <=z, x>=0
%
%  a_i denotes a training sample,
%      and a_i' corresponds to the i-th row of the data matrix A
%
%  y_i (either 1 or -1) is the response
%
%  p_i= 1/ (1+ exp(-y_i (x' * a_i + c) ) ) denotes the probability
%
%  weight_i denotes the weight for the i-th sample
%
% For detailed description of the function, please refer to the Manual.
%
%% Related papers
%
% [1]  Jun Liu and Jieping Ye, Efficient Euclidean Projections
%      in Linear Time, ICML 2009.
%
% [2]  Jun Liu and Jieping Ye, Sparse Learning with Efficient Euclidean
%      Projections onto the L1 Ball, Technical Report ASU, 2008.
%
% [3]  Jun Liu, Jianhui Chen, and Jieping Ye,
%      Large-Scale Sparse Logistic Regression, KDD, 2009.
%
%% ------------   History --------------------
%
% First version on August 10, 2009.
%
% September 5, 2009: adaptive line search is added
%
% For any problem, please contact Jun Liu (j.liu@asu.edu)

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/L1;

m=1000;  n=1000;    % The data matrix is of size m x n

randNum=1;          % a random number

% ---------------------- generate random data ----------------------
randn('state',(randNum-1)*3+1);
A=randn(m,n);         % the data matrix

y=[ones(n/2,1);...
    -ones(n/2, 1)];  % the response

z=40;                % the radius of the L1 ball

%----------------------- Set optional items -----------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% Termination
opts.tFlag=5;       % run .maxIter iterations
opts.maxIter=40;    % maximum number of iterations

% Normalization
opts.nFlag=0;       % without normalization

% Regularization
opts.rsL2=0;        % the squared two norm term

% Group Property
opts.sWeight=[1,1]; % set the weight for positive and negative samples

%----------------------- Run the code LogisticC -----------------------
[x, c, funVal]=nnLogisticC(A, y, z, opts);

figure;
plot(funVal);
xlabel('Iteration (i)');
ylabel('The objective function value');

% --------------------- compute the pathwise solutions ----------------
opts.fName='nnLogisticC';      % set the function name to 'LogisticC'
Z=[10, 20, 30, 40];          % set the parameters

% run the function pathSolutionLogistic
fprintf('\n Compute the pathwise solutions, please wait...');
[X,C]=pathSolutionLogistic(A, y, Z, opts);

author: zzstefan

repo: BrainNetClass


text:
% predicting numerical values using Linear Regression
clear all;
format long
disp('===== Linear Regression ====');
disp('Reading featur vector');

Indices = crossvalind('Kfold', 26548, 10);
figure;

subplot(3,2,1);

for feat = 1:3
      MSEarray =[];
    elapsedarray =[];
    for crossvalidateIter = 1:10
        (fprintf('%d',crossvalidateIter));
    featurs = csvread('data\forWeka_featuresonly.csv');
    num_data = size(featurs,1); %5000;
    %disp(sprintf('Number of datapoints %d',num_data))

    possiblefeaturizations =  {'bernouli', 'tfidf','multinomial'};
    %featurization = 'bernouli'%'tfidf'%'tfidf'%'multinomial'%'tfidf' %'multinomial'; % 'bernouli', 'tfidf'
    featurization  = possiblefeaturizations{feat};

    featurs = featurs(:,2:size(featurs,2));
    if strcmp(featurization,'multinomial')
        %just pass
    elseif strcmp(featurization,'bernouli')
        featurs = bernoulli(featurs);
    elseif strcmp(featurization,'tfidf')
        featurs = tfidf(featurs);
    end

     size_training = floor(.9*num_data);

        trainingset = featurs(Indices~=crossvalidateIter,:);
        testset = featurs(Indices==crossvalidateIter,:);

    %disp('Splitting up data into training/test sets');
    [num,txt,raw] = xlsread('data\final104.xls');

    % reading the description of each shoe
    descriptions = raw(2:size(raw,1),2);
    style_ratings = num(1:size(num,1),1);
    comfort_ratings = num(1:size(num,1),4);
    overal_ratings = num(1:size(num,1),5);

    % only take m data points
    m=num_data;
    descriptions = descriptions(1:m);
    style_ratings = style_ratings(1:m);
    comfort_ratings = comfort_ratings(1:m);
    overal_ratings = overal_ratings(1:m);

    responsevals = [style_ratings, comfort_ratings, overal_ratings];

        responsevals_training = responsevals(Indices~=crossvalidateIter,:);
        responsevals_test = responsevals(Indices==crossvalidateIter,:);

    %disp('Linear Regression');
    %

    tic;

    predictions = [];
    actual = [];
    for i =1:3
        a = responsevals_training(:,i);
        b = responsevals_test(:,i)';
        regresscoeff = regress(a, trainingset);
        C2 = (regresscoeff'*(testset'));
        predictions = [predictions, C2'];
        actual = [actual, b'];
    end

    MSE = mean(sum(((predictions-actual).^2)'))
     elapsed = toc;
        MSEarray = [MSEarray MSE];
        elapsedarray = [elapsedarray elapsed];

    end
    featurization  = possiblefeaturizations{feat}
    fprintf('Avergae MSE for Linear Regression is %0.10f\n', mean(MSEarray));
    subplot(3,2,feat*2-1)
    plot(MSEarray,'-o');
    title(strcat('MSE for Linear Regression (',featurization,') ', sprintf(' avergae MSE = %0.10f\n', mean(MSEarray))));
    xlabel('10 fold cross-validation (iteration no)')
    ylabel('MSE')
    subplot(3,2,feat*2)
    plot(elapsedarray,'-o');
    title(strcat('Elapsed time for Linear Regression (',featurization, ') ' , sprintf(' avergae elapsed time = %0.10f\n', mean(elapsedarray))));
    xlabel('10 fold cross-validation (iteration no)')
    ylabel('Elapsed Time')
    drawnow;
end

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function [alpha]=learningBasedMatting(imdata,mask)
% [alpha]=learningBasedMatting(im,mask) computes matte given image and
% scribble mask
%
% Input arguments:
% imdata:   MxNxd image data matrix. Image size is MxN, and the number of features is d. Value
%           range is within [0 255]
% mask:     MxN matrix specifying scribbles, with 1 foreground, -1 background
%           and 0 otherwise
%
% Output arguments:
% alpha:    (MxN) matrix of alpha solution
%
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% @InProceedings{ZhengICCV09,
%   author = {Yuanjie Zheng and Chandra Kambhamettu},
%   title = {Learning Based Digital Matting},
%   booktitle = {The 20th IEEE International Conference on Computer Vision},
%   year = {2009},
%   month = {September--October}
% }
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% Copyright Yuanjie Zheng @ PICSL @ UPenn on May 9th, 2011
% zheng.vision@gmail.com
% http://sites.google.com/site/zhengvision/
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

% get laplacian matrix
winsz=3;
c=800;
lambda=0.0000001;

L=getLap(imdata,winsz,mask,lambda);

% get regularization matrix
C=getC(mask,c);

% get already known values
alpha_star=getAlpha_star(mask);

% solve alpha
alpha=solveQurdOpt(L,C,alpha_star);

author: Sable

repo: mcbench-benchmarks


text:
function [rec,prec,ap] = VOCevalcls(VOCopts,id,cls,draw)

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

% load results
[ids,confidence]=textread(sprintf(VOCopts.clsrespath,id,cls),'%s %f');

% map results to ground truth images
out=ones(size(gt))*-inf;
tic;
for i=1:length(ids)
    % display progress
    if toc>1
        fprintf('%s: pr: %d/%d\n',cls,i,length(ids));
        drawnow;
        tic;
    end

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

% compute precision/recall

[so,si]=sort(-out);
tp=gt(si)>0;
fp=gt(si)<0;

fp=cumsum(fp);
tp=cumsum(tp);
rec=tp/sum(gt>0);
prec=tp./(fp+tp);

% compute average precision

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

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

author: CrazyStoneonRoad

repo: TGRS-HRRSD-Dataset


text:
%% housekeeping
gentle_clear()
close all
clc

%% load and transform the data
db=xlsread('usmacro2.xlsx');

db=ts('1955Q1',db(:,4:end),{'Y','P','R','C','N','I','E'});

db=pages2struct(db);

%% Rise the madel

m=rise('model10');

%% priors
priors=struct();
priors.beta={0.98,0.1,0.999};
priors.kappa={0.1,0.001,5};
priors.psi={0.5,0.001,1};
priors.rhou={0.5,-0.999,0.999};
priors.rhog={0.5,-0.999,0.999};
priors.rhoe={0.5,-0.999,0.999};
priors.sigu={0.1,0.001,15};
priors.sigg={0.1,0.001,15};
priors.sige={0.1,0.001,15};
%% estimate model

mest=estimate(m,'data_demean',true,'data',db,'priors',priors,...
    'estim_start_date','1973Q2','estim_end_date','2015Q4');

author: jmaih


text:
% this program is designed to correct the monthly variability for weather
% generator generated daily precip according to the monthly precip generated
% by FFT using a simple linear relationship

function [monthly_corrected_precip]=monthly_precip_correction(filenameout)
% load WG generated data
load(filenameout)
n=size(gP,1);

% load FFT generated monthly precip
monthly_FFT=zeros(n,12);
load('Pnew1')
monthly_FFT(:,1)=Pnew';
load('Pnew2')
monthly_FFT(:,2)=Pnew';
load('Pnew3')
monthly_FFT(:,3)=Pnew';
load('Pnew4')
monthly_FFT(:,4)=Pnew';
load('Pnew5')
monthly_FFT(:,5)=Pnew';
load('Pnew6')
monthly_FFT(:,6)=Pnew';
load('Pnew7')
monthly_FFT(:,7)=Pnew';
load('Pnew8')
monthly_FFT(:,8)=Pnew';
load('Pnew9')
monthly_FFT(:,9)=Pnew';
load('Pnew10')
monthly_FFT(:,10)=Pnew';
load('Pnew11')
monthly_FFT(:,11)=Pnew';
load('Pnew12')
monthly_FFT(:,12)=Pnew';

% calculate the monthly precip for the generated data
monthly_generated=zeros(n,12);
for i=1:n
    monthly_generated(i,1)=sum(gP(i,1:31));
    monthly_generated(i,2)=sum(gP(i,32:59));
    monthly_generated(i,3)=sum(gP(i,60:90));
    monthly_generated(i,4)=sum(gP(i,91:120));
    monthly_generated(i,5)=sum(gP(i,121:151));
    monthly_generated(i,6)=sum(gP(i,152:181));
    monthly_generated(i,7)=sum(gP(i,182:212));
    monthly_generated(i,8)=sum(gP(i,213:243));
    monthly_generated(i,9)=sum(gP(i,244:273));
    monthly_generated(i,10)=sum(gP(i,274:304));
    monthly_generated(i,11)=sum(gP(i,305:334));
    monthly_generated(i,12)=sum(gP(i,335:365));
end

% calculate the ratio of FFT generated data to WG generated data (monthly ratio)
monthly_ratio=zeros(n,12);
for i=1:n
    for j=1:12
        if monthly_generated(i,j)==0
            monthly_generated(i,j)=monthly_FFT(i,j);
        else monthly_ratio(i,j)=monthly_FFT(i,j)/monthly_generated(i,j);
        end
    end
end
monthly_ratio=monthly_ratio';
monthly_ratio=reshape(monthly_ratio,[],1);

% extend the above monthly ratio to daily scale,the data in each month are the same
monthly_ratio2=zeros(365*n,1);
j=0;
for i=0:365:365*n-1
    monthly_ratio2(i+1:i+31,1)=monthly_ratio(j+1,1);
    monthly_ratio2(i+32:i+59,1)=monthly_ratio(j+2,1);
    monthly_ratio2(i+60:i+90,1)=monthly_ratio(j+3,1);
    monthly_ratio2(i+91:i+120,1)=monthly_ratio(j+4,1);
    monthly_ratio2(i+121:i+151,1)=monthly_ratio(j+5,1);
    monthly_ratio2(i+152:i+181,1)=monthly_ratio(j+6,1);
    monthly_ratio2(i+182:i+212,1)=monthly_ratio(j+7,1);
    monthly_ratio2(i+213:i+243,1)=monthly_ratio(j+8,1);
    monthly_ratio2(i+244:i+273,1)=monthly_ratio(j+9,1);
    monthly_ratio2(i+274:i+304,1)=monthly_ratio(j+10,1);
    monthly_ratio2(i+305:i+334,1)=monthly_ratio(j+11,1);
    monthly_ratio2(i+335:i+365,1)=monthly_ratio(j+12,1);
    j=j+12;
end

% generate the years, months and days
daily_generated=zeros(365*n+1,6);
for i=1:366
    daily_generated(i,:)=datevec(i);
end
% delete the 60th row because it is the Feb 29th
daily_generated(60,:)=[];
% delete the fifthly and sixthly columns
daily_generated(:,6)=[];
daily_generated(:,5)=[];
% extend the date of the first to all years
j=1;
for i=366:365:365*n
    daily_generated(i:i+365-1,:)=daily_generated(1:365,:);
    daily_generated(i:i+365-1,1)=j;
    j=j+1;
end

% correct the WG generatd data
[n,m]=size(gP);
gP=gP';
j=1;
for i=1:365:m*n
    daily_generated(i:i+364,4)=gP(:,j);
    j=j+1;
end

monthly_adjust=zeros(size(daily_generated));
monthly_adjust(:,1:3)=daily_generated(:,1:3);
for i=1:365*n
    monthly_adjust(i,4)=daily_generated(i,4)*monthly_ratio2(i,1);
end
monthly_corrected_precip=monthly_adjust;
save('monthly_corrected_precip','monthly_corrected_precip')


text:
function [out] = soilmoisture_1(S1,S1max,S2,S2max)
%soilmoisture_1

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Water rebalance to equal relative storage (2 stores)
% Constraints:  -
% @(Inputs):    S1    - current storage in S1 [mm]
%               S1max - maximum storage in S1 [mm]
%               S2    - current storage in S2 [mm]
%               S2max - maximum storage in S2 [mm]

out = ((S2.*S1max-S1.*S2max)/(S1max+S2max)).*smoothThreshold_storage_logistic(S1./S1max,S2./S2max);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
%compute velmag_centralhead
function [data,units]=compute_velmagcentralhead(trx,n)
larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
velmagcentralhead=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    velmagcentralhead{1,i}=bsxfun(@hypot,trx(larva).dxcentralhead_mm,trx(larva).dycentralhead_mm);
end

units=parseunits('rad');
data=velmagcentralhead;

author: kristinbranson

repo: JAABA


text:
function b = r8sd_mxv ( n, ndiag, offset, a, x )

%*****************************************************************************80
%
%% R8SD_MXV multiplies a R8SD matrix times a vector.
%
%  Discussion:
%
%    The R8SD storage format is for symmetric matrices whose only nonzero entries
%    occur along a few diagonals, but for which these diagonals are not all
%    close enough to the main diagonal for band storage to be efficient.
%
%    In that case, we assign the main diagonal the offset value 0, and
%    each successive superdiagonal gets an offset value 1 higher, until
%    the highest superdiagonal (the A(1,N) entry) is assigned the offset N-1.
%
%    Assuming there are NDIAG nonzero diagonals (ignoring subdiagonals!),
%    we then create an array B that has N rows and NDIAG columns, and simply
%    "collapse" the matrix A to the left:
%
%  Example:
%
%    The "offset" value is printed above each column.
%
%    Original matrix               New Matrix
%
%       0   1   2   3   4   5       0   1   3   5
%
%      11  12   0  14   0  16      11  12  14  16
%      21  22  23   0  25   0      22  23  25  --
%       0  32  33  34   0  36      33  34  36  --
%      41   0  43  44  45   0      44  45  --  --
%       0  52   0  54  55  56      55  56  --  --
%      61   0  63   0  65  66      66  --  --  --
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer NDIAG, the number of diagonals that are stored.
%    NDIAG must be at least 1 and no more than N.
%
%    Input, integer OFFSET(NDIAG), the offsets for the diagonal storage.
%
%    Input, real A(N,NDIAG), the R8SD matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(N), the product A * x.
%
  b(1:n) = 0.0E+00;

  for i = 1 : n
    for jdiag = 1 : ndiag
      j = i + offset(jdiag);
      if ( 1 <= j & j <= n )
        b(i) = b(i) + a(i,jdiag) * x(j);
        if ( offset(jdiag) ~= 0 )
          b(j) = b(j) + a(i,jdiag) * x(i);
        end
      end
    end
  end

  return
end


text:
function results = test_fbcca(eeg, list_freqs, fs, num_harms, num_fbs)
% Steady-state visual evoked potentials (SSVEPs) detection using the filter
% bank canonical correlation analysis (FBCCA)-based method [1].
%
% function results = test_fbcca(eeg, list_freqs, fs, num_harms, num_fbs)
%
% Input:
%   eeg             : Input eeg data
%                     (# of targets, # of channels, Data length [sample])
%   list_freqs      : List for stimulus frequencies
%   fs              : Sampling frequency
%   num_harms       : # of harmonics
%   num_fbs         : # of filters in filterbank analysis
%
% Output:
%   results         : The target estimated by this method
%
% Reference:
%   [1] X. Chen, Y. Wang, S. Gao, T. -P. Jung and X. Gao,
%       "Filter bank canonical correlation analysis for implementing a
%        high-speed SSVEP-based brain-computer interface",
%       J. Neural Eng., vol.12, 046008, 2015.
%
% Masaki Nakanishi, 22-Dec-2017
% Swartz Center for Computational Neuroscience, Institute for Neural
% Computation, University of California San Diego
% E-mail: masaki@sccn.ucsd.edu

if nargin < 3
    error('stats:test_fbcca:LackOfInput', 'Not enough input arguments.');
end

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

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

fb_coefs = [1:num_fbs].^(-1.25)+0.25;

[num_targs, ~, num_smpls] = size(eeg);
y_ref = cca_reference(list_freqs, fs, num_smpls, num_harms);
for targ_i = 1:1:num_targs
     test_tmp = squeeze(eeg(targ_i, :, :));
     for fb_i = 1:1:num_fbs
         testdata = filterbank(test_tmp, fs, fb_i);
         for class_i = 1:1:num_targs
             refdata = squeeze(y_ref(class_i, :, :));
             [~,~,r_tmp] = canoncorr(testdata', refdata');
             r(fb_i, class_i) = r_tmp(1,1);
         end % class_i
     end % fb_i
     rho = fb_coefs*r;
    [~, tau] = max(rho);
    results(targ_i) = tau;
end % targ_i

function [ y_ref ] = cca_reference(list_freqs, fs, num_smpls, num_harms)
% Generate reference signals for the canonical correlation analysis (CCA)
% -based steady-state visual evoked potentials (SSVEPs) detection [1, 2].
%
% function [ y_ref ] = cca_reference(listFreq, fs,  nSmpls, nHarms)
%
% Input:
%   listFreq        : List for stimulus frequencies
%   fs              : Sampling frequency
%   nSmpls          : # of samples in an epoch
%   nHarms          : # of harmonics
%
% Output:
%   y_ref           : Generated reference signals
%                    (# of targets, 2*# of channels, Data length [sample])
%
% Reference:
%   [1] Z. Lin, C. Zhang, W. Wu, and X. Gao,
%       "Frequency Recognition Based on Canonical Correlation Analysis for
%        SSVEP-Based BCI",
%       IEEE Trans. Biomed. Eng., 54(6), 1172-1176, 2007.
%   [2] G. Bin, X. Gao, Z. Yan, B. Hong, and S. Gao,
%       "An online multi-channel SSVEP-based brain-computer interface using
%        a canonical correlation analysis method",
%       J. Neural Eng., 6 (2009) 046002 (6pp).
%
% Masaki Nakanishi, 28-Jul-2016
% Swartz Center for Computational Neuroscience, Institute for Neural
% Computation, University of California San Diego
% E-mail: masaki@sccn.ucsd.edu

if nargin < 3
    error('stats:cca_reference:LackOfInput',...
        'Not enough input arguments.');
end

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

num_freqs = length(list_freqs);
tidx = (1:num_smpls)/fs;
for freq_i = 1:1:num_freqs
    tmp = [];
    for harm_i = 1:1:num_harms
        stim_freq = list_freqs(freq_i);
        tmp = [tmp;...
            sin(2*pi*tidx*harm_i*stim_freq);...
            cos(2*pi*tidx*harm_i*stim_freq)];
    end % harm_i
    y_ref(freq_i, 1:2*num_harms, 1:num_smpls) = tmp;
end % freq_i

author: mnakanishi

repo: TRCA-SSVEP


text:
function [OutputError,ima,frame_nr_display] = PIVlab_capture_OPTOcam_calibration_image(img_amount,exposure_time,ROI_OPTOcam)
OutputError=0;
hgui=getappdata(0,'hgui');
%% Prepare camera
try
    delete(imaqfind); %clears all previous videoinputs
    hwinf = imaqhwinfo;
    %imaqreset
catch
    errordlg('Error: Image Acquisition Toolbox not available!','Error!','modal')
    disp('Error: Image Acquisition Toolbox not available!')
end

info = imaqhwinfo(hwinf.InstalledAdaptors{1});

if strcmp(info.AdaptorName,'gentl')
    disp('gentl adaptor found.')
else
    disp('ERROR: gentl adaptor not found. Please got to Matlab file exchange and search for "GenICam Interface " to install it.')
	errordlg('ERROR: gentl adaptor not found. Please got to Matlab file exchange and search for "GenICam Interface " to install it.','Error!','modal')
end

try
	OPTOcam_name = info.DeviceInfo.DeviceName;
catch
	errordlg('Error: Camera not found! Is it connected?','Error!','modal')
end
disp(['Found camera: ' OPTOcam_name])
OPTOcam_supported_formats = info.DeviceInfo.SupportedFormats;
OPTOcam_vid = videoinput(info.AdaptorName,info.DeviceInfo.DeviceID,'Mono12'); %calibration image in 12 bit always.

OPTOcam_settings = get(OPTOcam_vid);
OPTOcam_settings.Source.DeviceLinkThroughputLimitMode = 'off';
OPTOcam_settings.PreviewFullBitDepth='On';
OPTOcam_vid.PreviewFullBitDepth='On';

triggerconfig(OPTOcam_vid, 'manual');
OPTOcam_settings.TriggerMode ='manual';
OPTOcam_settings.Source.TriggerMode ='Off';
OPTOcam_settings.Source.ExposureMode ='Timed';
OPTOcam_settings.Source.ExposureTime =exposure_time;

ROI_OPTOcam=[ROI_OPTOcam(1)-1,ROI_OPTOcam(2)-1,ROI_OPTOcam(3),ROI_OPTOcam(4)]; %unfortunaletly different definitions of ROI in pco and OPTOcam.
OPTOcam_vid.ROIPosition=ROI_OPTOcam;

OPTOcam_settings.Source.ReverseX = 'True';
OPTOcam_settings.Source.ReverseY = 'True';
OPTOcam_gain = getappdata(hgui,'OPTOcam_gain');
if isempty (OPTOcam_gain)
    OPTOcam_gain=0;
end
OPTOcam_settings.Source.Gain = OPTOcam_gain;

%% prapare axis

crosshair_enabled = getappdata(hgui,'crosshair_enabled');
sharpness_enabled = getappdata(hgui,'sharpness_enabled');
PIVlab_axis = findobj(hgui,'Type','Axes');

%image_handle_OPTOcam=imagesc(zeros(OPTOcam_settings.VideoResolution(2),OPTOcam_settings.VideoResolution(1)),'Parent',PIVlab_axis,[0 2^8]);

image_handle_OPTOcam=imagesc(zeros(ROI_OPTOcam(4),ROI_OPTOcam(3)),'Parent',PIVlab_axis,[0 2^12]);

setappdata(hgui,'image_handle_OPTOcam',image_handle_OPTOcam);

frame_nr_display=text(100,100,'Initializing...','Color',[1 1 0]);
colormap default %reset colormap steps
new_map=colormap('gray');
new_map(1:3,:)=[0 0.2 0;0 0.2 0;0 0.2 0];
new_map(end-2:end,:)=[1 0.7 0.7;1 0.7 0.7;1 0.7 0.7];
colormap(new_map);axis image;
set(gca,'ytick',[])
set(gca,'xtick',[])
colorbar

%% get images
OPTOcam_vid.FramesPerTrigger = 1;
set(frame_nr_display,'String','');
preview(OPTOcam_vid,image_handle_OPTOcam)
clim([0 2^12]); %seems to be a workaround to force preview to show full data range...
displayed_img_amount=0;
while getappdata(hgui,'cancel_capture') ~=1 && displayed_img_amount < img_amount
    ima = image_handle_OPTOcam.CData;%*16; %stretch 12 bit to 16 bit
    %% sharpness indicator
    sharpness_enabled = getappdata(hgui,'sharpness_enabled');
    if sharpness_enabled == 1 % sharpness indicator
        textx=1240;
        texty=950;
        [~,~] = PIVlab_capture_sharpness_indicator (ima,textx,texty);
    else
        delete(findobj('tag','sharpness_display_text'));
    end
    crosshair_enabled = getappdata(hgui,'crosshair_enabled');
    if crosshair_enabled == 1 %cross-hair
        %% cross-hair
        locations=[0.15 0.5 0.85];
        half_thickness=1;
        brightness_incr=101;
        ima_ed=ima;
        old_max=max(ima(:));
        for loca=locations
            %vertical
            ima_ed(:,round(size(ima,2)*loca)-half_thickness:round(size(ima,2)*loca)+half_thickness)=ima_ed(:,round(size(ima,2)*loca)-half_thickness:round(size(ima,2)*loca)+half_thickness)+brightness_incr;
            %horizontal
            ima_ed(round(size(ima,1)*loca)-half_thickness:round(size(ima,1)*loca)+half_thickness,:)=ima_ed(round(size(ima,1)*loca)-half_thickness:round(size(ima,1)*loca)+half_thickness,:)+brightness_incr;
        end
        ima_ed(ima_ed>old_max)=old_max;
        set(image_handle_OPTOcam,'CData',ima_ed);
    end
    %% HISTOGRAM
    if getappdata(hgui,'hist_enabled')==1
        if isvalid(image_handle_OPTOcam)
            hist_fig=findobj('tag','hist_fig');
            if isempty(hist_fig)
                hist_fig=figure('numbertitle','off','MenuBar','none','DockControls','off','Name','Live histogram','Toolbar','none','tag','hist_fig','CloseRequestFcn', @HistWindow_CloseRequestFcn);
            end
            if ~exist ('old_hist_y_limits','var')
                old_hist_y_limits =[0 35000];
            else
                if isvalid(hist_obj)
                    old_hist_y_limits=get(hist_obj.Parent,'YLim');
                end
            end
            hist_obj=histogram(ima(1:2:end,1:2:end),'Parent',hist_fig,'binlimits',[0 2^12]);
        end
        %lowpass hist y limits for better visibility
        if ~exist ('new_hist_y_limits','var')
            new_hist_y_limits =[0 35000];
        end
        new_hist_y_limits=get(hist_obj.Parent,'YLim');
        set(hist_obj.Parent,'YLim',(new_hist_y_limits*0.5 + old_hist_y_limits*0.5))
    else
        hist_fig=findobj('tag','hist_fig');
        if ~isempty(hist_fig)
            close(hist_fig)
        end
    end
    drawnow limitrate;
    %% Autofocus
    %% Lens control
    %Sowieso machen: Nicht lineare schritte für die anzufahrenden fokuspositionen. Diese Liste vorher ausrechnen und dann nur index anspringen

    autofocus_enabled = getappdata(hgui,'autofocus_enabled');

    if autofocus_enabled == 1
        delaycounter=delaycounter+1;
    else
        delaycounter=0;
        delaycounter2=0;
        delay_time_1=tic;
    end
    %immer mehrere Bilder abfragen nachdem fokus verstellt wurde.... nicht nur eins, sondern z.B. drei Davon nur das letzte per sharpness beurteilen

    delay_time= 0.5; %1 seconds delay between measurements %350000 / exposure_time;
    if autofocus_enabled == 1
        if delaycounter>10 %wait 10 images before starting autofocus. Needed so that servo can reach target position
            focus_start = getappdata(hgui,'focus_servo_lower_limit');
            focus_end = getappdata(hgui,'focus_servo_upper_limit');
            amount_of_raw_steps=20;
            fine_step_resolution_increase = 8;
            focus_step_raw=round(abs(focus_end - focus_start)/amount_of_raw_steps);% in microseconds)
            focus_step_fine=round(1/fine_step_resolution_increase*(abs(focus_end - focus_start)/amount_of_raw_steps));% in microseconds)
            if ~exist('sharpness_focus_table','var') || isempty(sharpness_focus_table) || isempty(sharp_loop_cnt)
                sharpness_focus_table=zeros(1,2);
                sharp_loop_cnt=0;
                focus=focus_start;
                raw_finished=0;
                aperture=getappdata(hgui,'aperture');
                lighting=getappdata(hgui,'lighting');
                PIVlab_capture_lensctrl(focus,aperture,lighting)
            end
            if raw_finished==0
                if focus < focus_end % maxialer focus = endanschlag. Bis zu dem wert wird von null gefahren
                    if toc(delay_time_1)>=delay_time %only every second image is taken for analysis. This gives more time to the servo to reach position
                        delay_time_1=tic;
                        sharp_loop_cnt=sharp_loop_cnt+1;
                        [sharpness,~] = PIVlab_capture_sharpness_indicator (ima,[],[]);
                        sharpness_focus_table(sharp_loop_cnt,1)=focus;
                        sharpness_focus_table(sharp_loop_cnt,2)=sharpness;
                        focus=focus+focus_step_raw;
                        PIVlab_capture_lensctrl(focus,aperture,lighting)		%kann steuern und aktuelle position ausgeben
						autofocus_notification(1)
                    else
                        %do nothing
                    end
                else
                    %assignin('base','sharpness_focus_table',sharpness_focus_table)
                    %find best focus
                    [r,~]=find(sharpness_focus_table == max(sharpness_focus_table(:,2)));
                    focus_peak=sharpness_focus_table(r(1),1);
                    disp(['Best raw focus: ' num2str(focus_peak)])
                    raw_finished=1;
                    %focus vs. distance is not linear!
                    focus_start_fine=focus_peak-6*focus_step_raw; %start of finer focussearch
                    focus_end_fine=focus_peak+3*focus_step_raw;
                    if focus_start_fine < focus_start
                        focus_start_fine = focus_start;
                    end
                    if focus_end_fine > focus_end
                        focus_end_fine = focus_end;
                    end
                    %original focus=focus_end_fine;
                    focus=focus_start_fine;
                    PIVlab_capture_lensctrl(focus,aperture,lighting)
                    sharp_loop_cnt=0;
                    raw_data=[sharpness_focus_table(:,1),normalize(sharpness_focus_table(:,2),'range')];
                    sharpness_focus_table=zeros(1,2);
                end
            end

            if raw_finished == 1
                delaycounter2=delaycounter2+1;
            else
                delaycounter2=0;
            end

            if raw_finished == 1
                delay_time= 0.35;
                if delaycounter2>10
                    %repeat with finer steps
                    %original if focus > focus_start_fine % maxialer focus = endanschlag. Bis zu dem wert wird von null gefahren
                    if focus < focus_end_fine % maxialer focus = endanschlag. Bis zu dem wert wird von null gefahren
                        if toc(delay_time_1)>=delay_time %only every second image is taken for analysis. This gives more time to the servo to reach position
                            delay_time_1=tic;
                            sharp_loop_cnt=sharp_loop_cnt+1;
                            [sharpness,~] = PIVlab_capture_sharpness_indicator (ima,[],[]);
                            sharpness_focus_table(sharp_loop_cnt,1)=focus;
                            sharpness_focus_table(sharp_loop_cnt,2)=sharpness;
                            %original focus=focus-focus_step_fine;
                            focus=focus+focus_step_fine;
                            PIVlab_capture_lensctrl(focus,aperture,lighting)		%kann steuern und aktuelle position ausgeben
							autofocus_notification(1)
                        else
                            %do nothing
                        end
                    else %fine focus search finished
                        %assignin('base','sharpness_focus_table',sharpness_focus_table)
                        %find best focus
                        [r,~]=find(sharpness_focus_table == max(sharpness_focus_table(:,2)));
                        focus_peak=sharpness_focus_table(r(1),1);
                        disp(['Best fine focus: ' num2str(focus_peak)])
                        PIVlab_capture_lensctrl(focus_end_fine,aperture,lighting)%backlash compensation
                        pause(0.5)
                        PIVlab_capture_lensctrl(focus_start_fine,aperture,lighting) %backlash compensation
                        pause(0.5)
                        PIVlab_capture_lensctrl(focus_peak,aperture,lighting) %set to best focus

                        setappdata(hgui,'autofocus_enabled',0); %autofocus am ende ausschalten

                        lens_control_window = getappdata(0,'hlens');
                        focus_edit_field=getappdata(lens_control_window,'handle_to_focus_edit_field');
                        set(focus_edit_field,'String',num2str(focus_peak)); %update
                        %setappdata(hgui,'cancel_capture',1); %stop recording....?
                        figure;plot(raw_data(:,1),raw_data(:,2))
                        hold on;plot(sharpness_focus_table(:,1),normalize(sharpness_focus_table(:,2),'range'));hold off
                        title('Focus search')
                        xlabel('Pulsewidth us')
                        ylabel('Sharpness')
                        legend('Coarse search','Fine search')
                        grid on

                    end
                end
            end
        end
	else
		autofocus_notification(0)
        sharpness_focus_table=[];
        sharp_loop_cnt=[];
    end

	if img_amount == 1
		if sum(ima(1:10,1,1)) ~=10 %check if the display was updated, if there is real camera data. I didnt find a more elegant way...
			displayed_img_amount=displayed_img_amount+1;
		end
	end

end
stoppreview(OPTOcam_vid)

function autofocus_notification(running)
auto_focus_active_hint=findobj('tag', 'auto_focus_active');
if running == 1

	hgui=getappdata(0,'hgui');
	PIVlab_axis = findobj(hgui,'Type','Axes');
	%image_handle_OPTOcam=getappdata(hgui,'image_handle_OPTOcam');
	postix=get(PIVlab_axis,'XLim');
	postiy=get(PIVlab_axis,'YLim');
	bg_col=get(auto_focus_active_hint,'BackgroundColor'); % Toggle background color while autofocus is active

	if ~isempty(bg_col)
		if  sum(bg_col)==0.75 %hint is currently displayed
			bg_col = [0.05 0.05 0.05];
		else
			bg_col = [0.25 0.25 0.25];
		end
		set(auto_focus_active_hint,'BackgroundColor',bg_col);
	else
		bg_col= [0.25 0.25 0.25];
		axes(PIVlab_axis);
		text(postix(2)/2,postiy(2)/2,'Autofocus running, please wait...','HorizontalAlignment','center','VerticalAlignment','middle','color','y','fontsize',24, 'BackgroundColor', bg_col,'tag','auto_focus_active','margin',10,'Clipping','on');

	end
else
	delete(auto_focus_active_hint);
end

function HistWindow_CloseRequestFcn(hObject,~)
hgui=getappdata(0,'hgui');
setappdata(hgui,'hist_enabled',0);
try
    delete(hObject);
catch
    delete(gcf);
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function [X, imageDim] = slaveLoadTinyBatch(tinyFileName, batchStart, batchSize)

  idx=batchStart:batchStart+batchSize-1;
  X = loadTinyImages(idx, tinyFileName);
  X = double(reshape(X, 32*32*3, length(idx))');

  imageDim = [32 32 3];
  %X = slaveRef(X);

author: huashiyiqike

repo: LSTM-MATLAB


text:
function hole_num = p14_hole_num ( )

%*****************************************************************************80
%
%% P14_HOLE_NUM counts the holes in problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer HOLE_NUM, the number of holes.
%
  hole_num = 1;

  return
end


text:
function [ bdry ]  = seg2bdry(seg, fmt)
if nargin<2, fmt = 'doubleSize'; end;

if ~strcmp(fmt,'imageSize') && ~strcmp(fmt,'doubleSize'),
    error('possible values for fmt are: imageSize and doubleSize');
end

[tx, ty, nch] = size(seg);

if nch ~=1,
    error('seg must be a scalar image');
end

bdry = zeros(2*tx+1, 2*ty+1);

edgels_v = ( seg(1:end-1, :) ~= seg(2:end, :) );
edgels_v(end+1, :) = 0;
edgels_h = ( seg(:, 1:end-1) ~= seg(:, 2:end) );
edgels_h(:, end+1) = 0;

bdry(3:2:end, 2:2:end) = edgels_v;
bdry(2:2:end, 3:2:end) = edgels_h;
bdry(3:2:end-1, 3:2:end-1)= max ( max(edgels_h(1:end-1, 1:end-1), edgels_h(2:end, 1:end-1)), max(edgels_v(1:end-1,1:end-1), edgels_v(1:end-1,2:end)) );

bdry(1, :) = bdry(2, :);
bdry(:, 1) = bdry(:, 2);
bdry(end, :) = bdry(end-1, :);
bdry(:, end) = bdry(:, end-1);

if strcmp(fmt,'imageSize'),
    bdry = bdry(3:2:end, 3:2:end);
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
% this is an experimentaal spiral sequence

fov=256e-3; Nx=96; Ny=Nx;  % Define FOV and resolution
sliceThickness=3e-3;             % slice thinckness
Nslices=1;
Oversampling=2; % by looking at the periphery of the spiral I would say it needs to be at least 2
phi=pi/2; % orientation of the readout e.g. for interleaving

% Set system limits
sys = mr.opts('MaxGrad',30,'GradUnit','mT/m',...
    'MaxSlew',120,'SlewUnit','T/m/s',...
    'rfRingdownTime', 30e-6, 'rfDeadtime', 100e-6, 'adcDeadTime', 10e-6);
seq=mr.Sequence(sys);          % Create a new sequence object
warning('OFF', 'mr:restoreShape'); % restore shape is not compatible with spirals and will throw a warning from each plot() or calcKspace() call

% Create fat-sat pulse
% (in Siemens interpreter from January 2019 duration is limited to 8.192 ms, and although product EPI uses 10.24 ms, 8 ms seems to be sufficient)
B0=2.89; % 1.5 2.89 3.0
sat_ppm=-3.45;
sat_freq=sat_ppm*1e-6*B0*sys.gamma;
rf_fs = mr.makeGaussPulse(110*pi/180,'system',sys,'Duration',8e-3,...
    'bandwidth',abs(sat_freq),'freqOffset',sat_freq);
gz_fs = mr.makeTrapezoid('z',sys,'delay',mr.calcDuration(rf_fs),'Area',1/1e-4); % spoil up to 0.1mm

% Create 90 degree slice selection pulse and gradient
[rf, gz] = mr.makeSincPulse(pi/2,'system',sys,'Duration',3e-3,...
    'SliceThickness',sliceThickness,'apodization',0.5,'timeBwProduct',4);

% define k-space parameters
deltak=1/fov;
kRadius = round(Nx/2);
kSamples=round(2*pi*kRadius)*Oversampling;
readoutTime = 4.2e-4;

% calculate a raw Archimedian spiral trajectory
clear ka;
ka(kRadius*kSamples+1)=1i; % init as complex
for c=0:kRadius*kSamples
    r=deltak*c/kSamples;
    a=mod(c,kSamples)*2*pi/kSamples;
    ka(c+1)=r*exp(1i*a);
end
ka=[real(ka); imag(ka)];
% calculate gradients and slew rates
[ga, sa]=mr.traj2grad(ka);

% limit analysis
safety_magrin=0.94; % we need that  otherwise we just about violate the slew rate due to the rounding errors
dt_gcomp=abs(ga)/(sys.maxGrad*safety_magrin)*sys.gradRasterTime;
dt_gabs=abs(ga(1,:)+1i*ga(2,:))/(sys.maxGrad*safety_magrin)*sys.gradRasterTime;
dt_scomp=sqrt(abs(sa)/(sys.maxSlew*safety_magrin))*sys.gradRasterTime;
dt_sabs=sqrt(abs(sa(1,:)+1i*sa(2,:))/(sys.maxSlew*safety_magrin))*sys.gradRasterTime;

figure;plot([dt_gabs; max(dt_gcomp); dt_sabs; max(dt_scomp)]');title('time stepping defined by gradient and slew-rate');

dt_smooth=max([dt_gabs;dt_sabs]);
dt_rough=max([dt_gcomp;dt_scomp]);

% apply the lower limit not to lose the trajectory detail
dt_min=4*sys.gradRasterTime/kSamples; % we want at least 4 points per revolution
dt_smooth0=dt_smooth;
dt_rough0=dt_rough;
dt_smooth(dt_smooth<dt_min)=dt_min;
dt_rough(dt_rough<dt_min)=dt_min;

figure;plot([dt_smooth0; dt_smooth; dt_rough0; dt_rough]');title('combined time stepping');

t_smooth=[0 cumsum(dt_smooth,2)];
t_rough=[0 cumsum(dt_rough,2)];

kopt_smooth=interp1(t_smooth, ka', (0:floor(t_smooth(end)/sys.gradRasterTime))*sys.gradRasterTime)';
kopt_rough=interp1(t_rough, ka', (0:floor(t_rough(end)/sys.gradRasterTime))*sys.gradRasterTime)';

% analyze what we've got
fprintf('duration orig %d us\n', round(1e6*sys.gradRasterTime*length(ka)));
fprintf('duration smooth %d us\n', round(1e6*sys.gradRasterTime*length(kopt_smooth)));
fprintf('duration rough %d us\n', round(1e6*sys.gradRasterTime*length(kopt_rough)));

[gos, sos]=mr.traj2grad(kopt_smooth);
[gor, sor]=mr.traj2grad(kopt_rough);

figure;plot([gos;abs(gos(1,:)+1i*gos(2,:))]');title('gradient with smooth (abs) constraint')
figure;plot([gor;abs(gor(1,:)+1i*gor(2,:))]');title('gradient with rough (component) constraint')

figure;plot([sos;abs(sos(1,:)+1i*sos(2,:))]');title('slew rate with smooth (abs) constraint')
figure;plot([sor;abs(sor(1,:)+1i*sor(2,:))]');title('slew rate with rough (component) constraint')

% Define gradients and ADC events
spiral_grad_shape=gos;
% Create 90 degree slice selection pulse and gradient
[rf, gz] = mr.makeSincPulse(pi/2,'system',sys,'Duration',3e-3,...
    'SliceThickness',sliceThickness,'apodization',0.5,'timeBwProduct',4);
gzReph = mr.makeTrapezoid('z',sys,'Area',-gz.area/2);

% calculate ADC
% round-down dwell time to 10 ns
adcTime = sys.gradRasterTime*size(spiral_grad_shape,2);
% actually it is trickier than that: the (Siemens) interpreter sequence
% per default will try to split the trajectory into segments <=1000 samples
% and every of these segments will have to have duration aligned to the
% gradient raster time
adcSamplesPerSegment=1000; % you may need to play with this number to fill the entire trajectory
adcSamplesDesired=kRadius*kSamples;
adcSegments=round(adcSamplesDesired/adcSamplesPerSegment);
adcSamples=adcSegments*adcSamplesPerSegment;
adcDwell=round(adcTime/adcSamples/100e-9)*100e-9; % on Siemens adcDwell needs to be aligned to 100ns (if my memory serves me right)
adcSegmentDuration=adcSamplesPerSegment*adcDwell; % with the 100 samples above and the 100ns alignment we automatically fullfill the segment alignment requirement
if mod(adcSegmentDuration, sys.gradRasterTime)>eps
    error('ADC segmentation model results in incorrect segment duration');
end
% update segment count
adcSegments=floor(adcTime/adcSegmentDuration);
adcSamples=adcSegments*adcSamplesPerSegment;
adc = mr.makeAdc(adcSamples,'Dwell',adcDwell,'Delay',mr.calcDuration(gzReph));%lims.adcDeadTime);

% extend spiral_grad_shape by repeating the last sample
% this is needed to accomodate for the ADC tuning delay
spiral_grad_shape = [spiral_grad_shape spiral_grad_shape(:,end)];

% readout grad
gx = mr.makeArbitraryGrad('x',spiral_grad_shape(1,:),'Delay',mr.calcDuration(gzReph));
gy = mr.makeArbitraryGrad('y',spiral_grad_shape(2,:),'Delay',mr.calcDuration(gzReph));

% spoilers
gz_spoil=mr.makeTrapezoid('z',sys,'Area',deltak*Nx*4);
gx_spoil=mr.makeExtendedTrapezoid('x','times',[0 mr.calcDuration(gz_spoil)],'amplitudes',[spiral_grad_shape(1,end),0]); %todo: make a really good spoiler
gy_spoil=mr.makeExtendedTrapezoid('y','times',[0 mr.calcDuration(gz_spoil)],'amplitudes',[spiral_grad_shape(2,end),0]); %todo: make a really good spoiler

% because of the ADC alignment requirements the sampling window possibly
% extends past the end of the trajectory (these points will have to be
% discarded in the reconstruction, which is no problem). However, the
% ramp-down parts and the Z-spoiler now have to be added to the readout
% block otherwise there will be a gap inbetween
% gz_spoil.delay=mr.calcDuration(gx);
% gx_spoil.delay=gz_spoil.delay;
% gy_spoil.delay=gz_spoil.delay;
% gx_combined=mr.addGradients([gx,gx_spoil], lims);
% gy_combined=mr.addGradients([gy,gy_spoil], lims);
% gz_combined=mr.addGradients([gzReph,gz_spoil], lims);

% Define sequence blocks
for s=1:Nslices
    seq.addBlock(rf_fs,gz_fs); % fat-sat
    rf.freqOffset=gz.amplitude*sliceThickness*(s-1-(Nslices-1)/2);
    seq.addBlock(rf,gz);
    seq.addBlock(mr.rotate('z',phi,gzReph,gx,gy,adc));
    seq.addBlock(mr.rotate('z',phi,gx_spoil,gy_spoil,gz_spoil));
    %seq.addBlock(gx_combined,gy_combined,gz_combined,adc);
end

% check whether the timing of the sequence is correct
[ok, error_report]=seq.checkTiming;

if (ok)
    fprintf('Timing check passed successfully\n');
else
    fprintf('Timing check failed! Error listing follows:\n');
    fprintf([error_report{:}]);
    fprintf('\n');
end

%
seq.setDefinition('FOV', [fov fov sliceThickness]);
seq.setDefinition('Name', 'spiral');
seq.setDefinition('MaxAdcSegmentLength', adcSamplesPerSegment);

seq.write('spiral.seq');   % Output sequence for scanner

% the sequence is ready, so let's see what we got
seq.plot();             % Plot sequence waveforms

%% k-space trajectory calculation
[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing] = seq.calculateKspacePP();

% plot k-spaces
figure; plot(t_ktraj, ktraj'); title('k-space components as functions of time'); % plot the entire k-space trajectory
figure; plot(ktraj(1,:),ktraj(2,:),'b'); % a 2D plot
hold;plot(ktraj_adc(1,:),ktraj_adc(2,:),'r.'); title('2D k-space');

% seq.install('siemens');

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
% DEMROBOTWIRELESS2 Wireless Robot data from University of Washington, without dynamics and without back constraints.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'robotWireless';
experimentNo = 2;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('ftc');
options.back = 'mlp';
options.backOptions = mlpOptions;

latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
capName = dataSetName;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model');

% Load the results and display dynamically.
fgplvmResultsDynamic(dataSetName, experimentNo, 'vector')

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef NodesCalculatorFactory < handle
    methods (Access = public, Static)
        function obj = create(cParams)
            switch cParams.nvert
                case 4
                    obj = QuadrilateralNodesCalculator(cParams);
                case 6
                    obj = HexagonalNodesCalculator(cParams);
            end
        end
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function F  = ...
    check_frequency_response(x,frequency_20,frequency_100, method)
% Checks the result of optimization by plotting phase characteristics at
% optimal parameters and computing the frequencies at -pi/2 phase shift.
% The frequency responses of a nonlinear system are obtained by processing
% the pulse transient characteristics with the FFT algorithm.
% Copyright 2010 MathWorks, Inc.

% frequency_20   - frequency (Hz) at phase shift in -pi/2 at 20% input signal
% frequency_100  - frequency (Hz) at phase shift in -pi/2 at 100% input signal
% method - method to check the results, can be set to 'FFT' or 'direct'.

if strcmpi(method, 'FFT')
    model = 'actuator_freq_testrig_pulse_FFT_method';
elseif strcmpi(method, 'direct')
    model = 'actuator_freq_testrig_direct_method';
else
    error('Unspecified method');
end

load_system(model);

assignin('base','gain', x(1));
assignin('base','time_const', x(2));
assignin('base','saturation', x(3));

sim(model);

if strcmpi(method, 'FFT')

    y_20 = yout(:,2);               % Pulse transient characteristic at 20% input
    y_100 = yout(:,1);              % Pulse transient characteristic at 100% input
    fs = 1000;                      % Sampling frequency
    n = length(y_20);               % Window length = Transform length
    y_20_fft = fft(y_20,n);         % Discrete Fourier Transform
    y_100_fft = fft(y_100,n);       % Discrete Fourier Transform
    f0 = (0:n/2-1)*(fs/n);          % Shifted frequency range, positive region
    y_20_0 = fftshift(y_20_fft);    % Shifted DFT at 20% input
    y_100_0 = fftshift(y_100_fft);  % Shifted DFT at 100% input
    % Phase characteristic at 20% input for positive frequencies after unwrap
    phase_20 = unwrap(angle(y_20_0(257:end)));
    % Phase characteristic at 100% input for positive frequencies after unwrap
    phase_100 = unwrap(angle(y_100_0(257:end)));

    % Computing frequency at 90 deg phase shift by interpolation of phase
    % characteristics
    frq_20 = interp1(phase_20,f0,-pi/2);
    frq_100 = interp1(phase_100,f0,-pi/2);

    % Computing frequency at phase shift in -pi/2 by interpolation of phase
    % characteristics
    frq_20 = interp1(phase_20,f0,-pi/2);
    frq_100 = interp1(phase_100,f0,-pi/2);
    % Computing errors
    error_20 = (frequency_20 - frq_20)/frequency_20 * 100;
    error_100 = (frequency_100 - frq_100)/frequency_100 * 100;

    disp(['******************************************************************']);
    disp(['     System Characteristics with FFT at Obtained Parameters']);
    disp(['     Frequencies at -90 deg phase shift and 20% input:']);
    disp(['Target:  ', num2str(frequency_20),' Hz', ...
        '  Found:  ',num2str(frq_20),' Hz', ...
        '   Error :   ',num2str(error_20),'%']);
    disp(['     Frequencies at -90 deg phase shift and 100% input']);
    disp(['Target:  ', num2str(frequency_100),' Hz', ...
        ' Found:     ',num2str(frq_100),' Hz', ...
        '   Error:   ',num2str(error_100),'%']);

    % Plotting the phase frequency characterictics
    semilogx(f0,phase_20*180/pi,'LineWidth',3);
    hold on
    semilogx(f0,phase_100*180/pi,'r--','LineWidth',3), grid on;
    xlabel('Frequency (Hz)','FontSize',14)
    ylabel('Phase (degrees)','FontSize',14)
    title('Frequency Response (Phase)','FontWeight','Bold','FontSize',16)
    legend({'20% Input' '100% Input'},'FontSize',12);
    axis([2 40 -100 0]);
    plot(frq_100,-90,'ro','MarkerSize',10,'MarkerFaceColor','r')
    plot(frq_20,-90,'bo','MarkerSize',10,'MarkerFaceColor','b')
    hold off

else    % Direct measurement. yout contains input and output values of
        % both signals

    % Phase angle at specified frequency and input signal in 100%
    phase_100 = phase_computation(yout(:,[1,2]), tout, frequency_100,4);
    % Phase angle at specified frequency and input signal in 20%
    phase_20 = phase_computation(yout(:,[3,4]), tout, frequency_20,6);

    error_20 = (phase_20 + 90)/90 * 100;
    error_100 = (phase_100 +90)/90 * 100;

    disp(['******************************************************************']);
    disp(['     System Characteristics after Direct Measurement at Obtained Parameters']);
    disp(['     Phase shift at ',num2str(frequency_20),'Hz',' and 20% input:']);
    disp(['Target:  ',' -90 deg', ...
        '  Found:  ',num2str(phase_20),' deg', ...
        '   Error :   ',num2str(error_20),'%']);
    disp(['     Phase shift at ',num2str(frequency_100),'Hz',' and 100% input']);
    disp(['Target:  ',' -90 deg', ...
        '  Found:   ',num2str(phase_100),' deg', ...
        '   Error:   ',num2str(error_100),'%']);
    disp(['******************************************************************']);

end
end

function [phase] = phase_computation(yout, tout, freq_at_90, N)
% Function computes phase angle by processing sinusoids stored in yout

% yout(:,1)- input; yout(:,2) - output
% tout - simulation step time
% freq_at_90 - frequency at which phase angle equals -90 deg
% N - number of periods to settle the transients

% Measurement time. Time in N periods is assumed to be enough to settle
% down the transients. The input signal at this time is equal zero

t_s = 1/freq_at_90 * N;             % Start time
t_end = t_s + 1/freq_at_90;         % End time

% yout index at measurement time
for j = 1:length(tout)- 1
    if t_s >= tout(j) && t_s < tout(j+1)
        k_start = j;                % First array index after start time
        t_corr_1 = tout(j) - t_s;   % Correction due to discretization
    end
    if t_end >= tout(j) && t_end < tout(j+1)
        k_end = j-1;                  % Last array index within period
    end
end
% Computing time the output signal crosses zero
for j = k_start : k_end
    if yout(j,2) <= 0 && yout(j+1,2) > 0
        t_tab_cross = tout(j);   % First value in the table after crossing
        % Approximating real crossing time
        t_corr_2 = -yout(j,2) / (yout(j+1,2) - yout(j,2)) * ...
            (tout(j+1) - tout(j));
    end
end
% Crossing time for the output signal
t_cross = t_tab_cross + t_corr_1 + t_corr_2;
% Phase angle
phase = -(t_cross - t_s) / (1/freq_at_90) * 360;

end

% EOF


text:
function pass = test_end(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

%% test end in infinite dimension:
% Scalar:
f = chebfun(@(x) sin(pi*x));
out = f(end);
pass(1) = isnumeric(out) && length(out) == 1 && abs(out) < eps;

% Array:
f = chebfun(@(x) [sin(pi*x), cos(pi*x), x]);
out = f(end);
pass(2) = isnumeric(out) && all(size(out) == [1, 3]) && ...
    norm(out - [0 -1 1], inf) < eps;

out = f(end,2);
pass(3) = isnumeric(out) && numel(out) == 1 && abs(out + 1) < eps;

% transpose
g = f.';
out = g(end);
pass(4) = isnumeric(out) && all(size(out) == [3, 1]) && ...
    norm(out - [0 -1 1].', inf) < eps;

%% test end in scalar dimension:
f = chebfun(@(x) [sin(pi*x), cos(pi*x), x]);
out = f(0,end);
pass(5) = isnumeric(out) && length(out) == 1&& abs(out) < eps;

g = f.';
out = f(0,end);
pass(6) = isnumeric(out) && length(out) == 1&& abs(out) < eps;

%% test use of ':'
f = chebfun(@(x) [sin(pi*x), cos(pi*x), x]);
out = f(end, :);
pass(7) = isnumeric(out) && all(size(out) == [1, 3]) && ...
    norm(out - [0 -1 1], inf) < eps;

out = f(:, end);
x = chebfun(@(x) x);
err = normest(out - x);
pass(8) = err < eps;

% Transpose:
g = f.';
out = g(end, :);
err = normest(out - x.');
pass(9) = err < eps;

out = g(:, end);
pass(10) = isnumeric(out) && all(size(out) == [3, 1]) && ...
    norm(out - [0 -1 1].', inf) < eps;

%% Test on singular function:

% Set a domain:
dom = [-2 7];

pow = -0.5;
op = @(x) (x - dom(2)).^pow.*sin(100*x);
f = chebfun(op, dom, 'exps', [0 pow], 'splitting', 'on');
out = f(end);
pass(11) = ( isnumeric(out) ) && all( size(out) == ones(1,2) );

%% Test for function defined on unbounded domain:

% Set a domain:
dom = [0 Inf];

op = @(x) 0.75+sin(10*x)./exp(x);
f = chebfun(op, dom, 'splitting', 'on');
out = f(end);
pass(12) = ( isnumeric(out) ) && all(size(out) == ones(1, 2) ) && ...
    ( abs(out - 0.75) < 1e2*eps );

end

function out = normest(f, dom)

% Generate a few random points to use as test values.
seedRNG(6178);
if ( nargin == 1 )
    x = 2 * rand(100, 1) - 1;
else
    x = sum(dom) * rand(10, 1) - dom(1);
end

out = norm(feval(f, x), inf);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks, deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=getINITModel(refModel, tissue, celltype, hpaData, arrayData, metabolomicsData, taskFile, useScoresForTasks, printReport, taskStructure, params, paramsFT)
% getINITModel_legacy
%   Generates a model using the INIT algorithm, based on proteomics and/or
%   transcriptomics and/or metabolomics and/or metabolic tasks. This is the original
%   implementation of tINIT, which is replaced by ftINIT.
%
%   Input:
%   refModel            a model structure. The model should be in the
%                       closed form (no exchange reactions open). Import
%                       using import(filename,false). If the model is not
%                       loaded using importModel, it might be that there
%                       is no "unconstrained" field. In that case,
%                       manually add the field like:
%                       model.unconstrained=false(numel(model.mets),1);
%   tissue              tissue to score for. Should exist in either
%                       hpaData.tissues or arrayData.tissues
%   celltype            cell type to score for. Should exist in either
%                       hpaData.celltypes or arrayData.celltypes for this
%                       tissue (opt, default is to use the best values
%                       among all the cell types for the tissue. Use [] if
%                       you want to supply more arguments)
%   hpaData             HPA data structure from parseHPA (opt if arrayData is
%                       supplied, default [])
%   arrayData           gene expression data structure (opt if hpaData is
%                       supplied, default [])
%       genes           cell array with the unique gene names
%       tissues         cell array with the tissue names. The list may not be
%                       unique, as there can be multiple cell types per tissue
%       celltypes       cell array with the cell type names for each tissue
%       levels          GENESxTISSUES array with the expression level for
%                       each gene in each tissue/celltype. NaN should be
%                       used when no measurement was performed
%       threshold       a single value or a vector of gene expression
%                       thresholds, above which genes are considered to be
%                       "expressed". (opt, by default, the mean expression
%                       levels of each gene across all tissues in arrayData
%                       will be used as the threshold values)
%       singleCells     binary value selecting whether to use the
%                       single-cell algorithm to identify expressed genes.
%                       If used, specify cell subpopulations in CELLTYPES
%                       (opt, default [])
%       plotResults     true if single cell probability distributions
%                       should be plotted (opt, default = False)
%   metabolomicsData    cell array with metabolite names that the model
%                       should produce (opt, default [])
%   taskFile            a task list in Excel format. See parseTaskList for
%                       details (opt, default [])
%   useScoresForTasks   true if the calculated reaction scored should be used as
%                       weights in the fitting to tasks (opt, default true)
%   printReport         true if a report should be printed to the screen
%                       (opt, default true)
%   taskStructure       task structure as from parseTaskList. Can be used
%                       as an alternative way to define tasks when Excel
%                       sheets are not suitable. Overrides taskFile (opt,
%                       default [])
%   params              parameter structure as used by getMILPParams. This is
%                       for the INIT algorithm. For the the MILP problems
%                       solved to fit tasks, see paramsFT (opt, default [])
%   paramsFT            parameter structure as used by getMILPParams. This is
%                       for the fitTasks step. For the INIT algorithm, see
%                       params (opt, default [])
%
%
%   Output:
%   model                   the resulting model structure
%   metProduction           array that indicates which of the
%                           metabolites in metabolomicsData that could be
%                           produced. Note that this is before the
%                           gap-filling process to enable defined tasks. To
%                           see which metabolites that can be produced in
%                           the final model, use canProduce.
%                           -2: metabolite name not found in model
%                           -1: metabolite found, but it could not be produced
%                           1: metabolite could be produced
%   essentialRxnsForTasks   cell array of the reactions which were
%                           essential to perform the tasks
%   addedRxnsForTasks       cell array of the reactions which were added in
%                           order to perform the tasks
%   deletedDeadEndRxns      cell array of reactions deleted because they
%                           could not carry flux (INIT requires a
%                           functional input model)
%   deletedRxnsInINIT       cell array of the reactions which were deleted by
%                           the INIT algorithm
%   taskReport              structure with the results for each task
%   	id                  cell array with the id of the task
%       description         cell array with the description of the task
%       ok                  boolean array with true if the task was successful
%       essential           cell array with cell arrays of essential
%                           reactions for the task
%       gapfill             cell array of cell arrays of reactions included
%                           in the gap-filling for the task
%
%   This is the main function for automatic reconstruction of models based
%   on the (t)INIT algorithm (PLoS Comput Biol. 2012;8(5):e1002518,
%   Mol Syst Biol. 2014;10:721). Not all settings are possible using this
%   function, and you may want to call the functions scoreModel, runINIT
%   and fitTasks individually instead.
%
%   NOTE: Exchange metabolites should normally not be removed from the model
%   when using this approach, since checkTasks/fitTasks rely on putting specific
%   constraints for each task. The INIT algorithm will remove exchange metabolites
%   if any are present. Use importModel(file,false) to import a model with
%   exchange metabolites remaining.
%
%   Usage: [model, metProduction, essentialRxnsForTasks, addedRxnsForTasks,...
%               deletedDeadEndRxns, deletedRxnsInINIT, taskReport]=...
%               getINITModel(refModel, tissue, celltype, hpaData, arrayData,...
%               metabolomicsData, taskFile, useScoresForTasks, printReport,...
%               taskStructure, params, paramsFT)

if nargin<3
    celltype=[];
else
    celltype=char(celltype);
end
if nargin<4
    hpaData=[];
end
if nargin<5
    arrayData=[];
end
if nargin<6
    metabolomicsData=[];
end
if nargin<7
    taskFile=[];
else
    taskFile=char(taskFile);
end
if nargin<8 || isempty(useScoresForTasks)
    useScoresForTasks=true;
end
if nargin<9 || isempty(printReport)
    printReport=true;
end
if nargin<10
    taskStructure=[];
end
if nargin<11
    params=[];
end
if nargin<12
    paramsFT=[];
end

%Check that the model is in the closed form
if ~isfield(refModel,'unconstrained')
    EM='Exchange metabolites should normally not be removed from the model when using getINITModel. Use importModel(file,false) to import a model with exchange metabolites remaining (see the documentation for details)';
    dispEM(EM);
end

%Create the task structure if not supplied
if any(taskFile) && isempty(taskStructure)
    taskStructure=parseTaskList(taskFile);
end

% sc-tINIT to identify confidence levels of gene expression
if ~isempty(arrayData) && isfield(arrayData,'singleCells')
    if arrayData.singleCells == 1
        % Check to ensure cell type is defined
        if ~isfield(arrayData,'celltypes')
            dispEM('arrayData must contain cell type information if sc-tINIT is to be used','false');
        end
        if ~ismember(upper(celltype),upper(arrayData.celltypes))
            dispEM('The cell type name does not match');
        end

        % Analyze only cell type of interest
        J= strcmpi(arrayData.celltypes,celltype);

        % Analyze only genes included in the reference model
        I=ismember(arrayData.genes,refModel.genes);

        % Convert expression data to population fractions
        binary_levels = arrayData.levels(I,J)~=0; % Classify each gene as detected (1) or not (0) in each cell
        cell_count_levels = sum(binary_levels,2); % Number of cells expressing each transcript
        cell_frac_levels = cell_count_levels/size(binary_levels,2); % Number of cells expressing each transcript

        % Bin cell_frac_counts manually
        x = 0:.01:1;
        for(i = 1:length(x))
            cell_frac_count(i) = sum(round(cell_frac_levels,2)==x(i));
        end

        % Fit four beta distributions
        cell_frac_count(cell_frac_count==0) = NaN; % Remove zeros from optimization
        cell_frac_count(1) = NaN; % Remove non-expressed genes from optimization
        x_lim = 1; % Somewhat arbitrary parameter to fit left tail of distr.
        myfun = @(par) nansum((cell_frac_count(1:find(x>=x_lim,1)) - ...
            abs(par(1))*betapdf(x(1:find(x>=x_lim,1)),abs(par(2)),abs(par(3))) - ...
            abs(par(4))*betapdf(x(1:find(x>=x_lim,1)),abs(par(5)),abs(par(6))) - ...
            abs(par(7))*betapdf(x(1:find(x>=x_lim,1)),abs(par(8)),abs(par(9))) - ...
            abs(par(10))*betapdf(x(1:find(x>=x_lim,1)),abs(par(11)),abs(par(12)))).^2);

        par0 = [4,2,100,7,2,30,7,5,20,5,15,20];
        opts = optimset('Display','off');
        [par,f_val] = fminsearch(myfun,par0,opts);
        par = abs(par);

        % Plot results
        if (isfield(arrayData,'plotResults'))
            if arrayData.plotResults == true
                figure(); hold on; plot(x,cell_frac_count,'ko','MarkerSize',5);
                plot(x,abs(par(1))*betapdf(x,abs(par(2)),abs(par(3))),'b-','LineWidth',1)
                plot(x,abs(par(4))*betapdf(x,abs(par(5)),abs(par(6))),'b-','LineWidth',1)
                plot(x,abs(par(7))*betapdf(x,abs(par(8)),abs(par(9))),'b-','LineWidth',1)
                plot(x,abs(par(10))*betapdf(x,abs(par(11)),abs(par(12))),'b-','LineWidth',1)
                plot(x,abs(par(1))*betapdf(x,abs(par(2)),abs(par(3))) + ...
                    abs(par(4))*betapdf(x,abs(par(5)),abs(par(6))) + ...
                    abs(par(7))*betapdf(x,abs(par(8)),abs(par(9))) + ...
                    abs(par(10))*betapdf(x,abs(par(11)),abs(par(12))),'-','Color',[.5 .5 .5],'LineWidth',2)
                xlabel('Expression Probability');ylabel('# of genes');set(gca,'FontSize',14,'LineWidth',1.25);
                title('Expression prediction','FontSize',18,'FontWeight','bold')
            end
        end

        % Score genes based on population expression (p = .05)
        exprs_cutoff_1 = find(cdf('beta',x,par(2),par(3)) >.95,1)-1; % Find index of no confidence genes
        exprs_cutoff_2 = find(cdf('beta',x,par(5),par(6)) >.95,1)-1; % Find index of low confidence genes
        exprs_cutoff_3 = find(cdf('beta',x,par(8),par(9)) >.95,1)-1; % Find index of low confidence genes
        exprs_cutoffs = sort([exprs_cutoff_1,exprs_cutoff_2,exprs_cutoff_3]);
        gene_scores = cell_frac_levels*0;
        gene_scores(cell_frac_levels <= x(exprs_cutoffs(1))) = 4; % Not detected
        gene_scores(logical((cell_frac_levels >= x(exprs_cutoffs(1))).*(cell_frac_levels < x(exprs_cutoffs(2))))) = 3; % Low detection
        gene_scores(logical((cell_frac_levels >= x(exprs_cutoffs(2))).*(cell_frac_levels < x(exprs_cutoffs(3))))) = 2; % Medium detection
        gene_scores(cell_frac_levels > x(exprs_cutoffs(3))) = 1; % High detection

        % Replace hpaData with singleCellData
        if printReport==true
            dispEM('Single cell data is not currently compatible with HPA data. \n         Replacing hpaData with single cell-based scoring.',false);
        end
        hpaData.genes = arrayData.genes;
        hpaData.tissues = arrayData.tissues;
        hpaData.celltypes = arrayData.celltypes;
        hpaData.levels = [{'High'},{'Medium'},{'Low'},{'None'}];
        hpaData.gene2Level = zeros(length(arrayData.genes),length(arrayData.celltypes));
        for i = 1:length(find(J))
            find_var = find(J,i);
            hpaData.gene2Level(I,find_var(end)) = gene_scores;
        end

        % Remove arrayData from the analysis (Might be a bad idea)
        clear arrayData
        arrayData=[];
    end
end

if printReport==true
    if any(celltype)
        fprintf(['***Generating model for: ' tissue ' - ' celltype '\n']);
    else
        fprintf(['***Generating model for: ' tissue '\n']);
    end
    if ~isempty(hpaData)
        fprintf('-Using HPA data\n');
    end
    if ~isempty(arrayData)
        fprintf('-Using array data\n');
    end
    if ~isempty(metabolomicsData)
        fprintf('-Using metabolomics data\n');
    end
    if ~isempty(taskFile) || ~isempty(taskStructure)
        fprintf('-Using metabolic tasks\n');
    end
    fprintf('\n');

    printScores(refModel,'Reference model statistics',hpaData,arrayData,tissue,celltype);
end

%Remove dead-end reactions to speed up the optimization and to
%differentiate between reactions removed by INIT and those that are
%dead-end
[~, deletedDeadEndRxns]=simplifyModel(refModel,true,false,true,true,true);
cModel=removeReactions(refModel,deletedDeadEndRxns,false,true);

%Store the connected model like this to keep track of stuff
if printReport==true
    printScores(cModel,'Pruned model statistics',hpaData,arrayData,tissue,celltype);
end

%If tasks have been defined, then go through them and get essential
%reactions
if ~isempty(taskStructure)
    [taskReport, essentialRxnMat]=checkTasks(cModel,[],printReport,true,true,taskStructure);

    essentialRxnsForTasks=cModel.rxns(any(essentialRxnMat,2));

    %Remove tasks that cannot be performed
    taskStructure(taskReport.ok==false)=[];
    if printReport==true
        printScores(removeReactions(cModel,setdiff(cModel.rxns,essentialRxnsForTasks),true,true),'Reactions essential for tasks',hpaData,arrayData,tissue,celltype);
    end
else
    essentialRxnsForTasks={};
end

%Score the connected model
[rxnScores, geneScores]=scoreModel(cModel,hpaData,arrayData,tissue,celltype);

%Run the INIT algorithm. The exchange reactions that are used in the final
%reactions will be open, which doesn't fit with the last step. Therefore
%delete reactions from the original model instead of taking the output. The
%default implementation does not constrain reversible reactions to only
%carry flux in one direction. Runs without the constraints on reversibility
%and with all output allowed. This is to reduce the complexity of the
%problem.
[~, deletedRxnsInINIT, metProduction]=runINIT(simplifyModel(cModel),rxnScores,metabolomicsData,essentialRxnsForTasks,0,true,false,params);
initModel=removeReactions(cModel,deletedRxnsInINIT,true,true);
if printReport==true
    printScores(initModel,'INIT model statistics',hpaData,arrayData,tissue,celltype);
    printScores(removeReactions(cModel,setdiff(cModel.rxns,deletedRxnsInINIT),true,true),'Reactions deleted by INIT',hpaData,arrayData,tissue,celltype);
end

%The full model has exchange reactions in it. fitTasks calls on fillGaps,
%which automatically removes exchange metabolites (because it assumes that
%the reactions are constrained when appropriate). In this case the
%uptakes/outputs are retrieved from the task sheet instead. To prevent
%exchange reactions being used to fill gaps, they are delete from the
%reference model here.
initModel.id='INITModel';

%If gaps in the model should be filled using a task list
if ~isempty(taskStructure)
    %Remove exchange reactions and reactions already included in the INIT
    %model
    refModelNoExc=removeReactions(refModel,union(initModel.rxns,getExchangeRxns(refModel)),true,true);

    %At this stage the model is fully connected and most of the genes with
    %good scores should have been included. The final gap-filling should
    %take the scores of the genes into account, so that "rather bad"
    %reactions are preferred to "very bad" reactions. However, reactions
    %with positive scores will be included even if they are not connected
    %in the current formulation. Therefore, such reactions will have to be
    %assigned a small negative score instead.
    if useScoresForTasks==true
        refRxnScores=scoreModel(refModelNoExc,hpaData,arrayData,tissue,celltype);
        [outModel, addedRxnMat]=fitTasks(initModel,refModelNoExc,[],true,min(refRxnScores,-0.1),taskStructure,paramsFT);
    else
        [outModel, addedRxnMat]=fitTasks(initModel,refModelNoExc,[],true,[],taskStructure,paramsFT);
    end
    if printReport==true
        printScores(outModel,'Functional model statistics',hpaData,arrayData,tissue,celltype);
        printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,arrayData,tissue,celltype);
    end

    addedRxnsForTasks=refModelNoExc.rxns(any(addedRxnMat,2));
else
    outModel=initModel;
    addedRxnMat=[];
    addedRxnsForTasks={};
end

%The model can now perform all the tasks defined in the task list. The
%algorithm cannot deal with gene-complexes at the moment. It is therefore
%ok to remove bad genes from a reaction (as long as at least one gene is
%kept)
model=outModel;

[~, I]=ismember(model.genes,cModel.genes); %All should be found
%This is a little weird way to make sure that only one bad gene is included
%if there are no good ones (since all -Inf==max(-Inf))
geneScores(isinf(geneScores))=-1000+rand(sum(isinf(geneScores)),1);

model.grRules(:)={''};
for i=1:numel(model.rxns)
    ids=find(model.rxnGeneMat(i,:));
    if numel(ids)>1
        scores=geneScores(I(ids));
        %Only keep the positive ones if possible
        model.rxnGeneMat(i,ids(~(scores>0 | scores==max(scores))))=0;
    end
    %Rewrite the grRules to be only OR
    if isfield(model,'grRules')
        J=find(model.rxnGeneMat(i,:));
        for j=1:numel(J)
            model.grRules{i}=[model.grRules{i} '(' model.genes{J(j)} ')'];
            if j<numel(J)
                model.grRules{i}=[model.grRules{i} ' or '];
            end
        end
    end
end

%Find all genes that are not used and delete them
I=sum(model.rxnGeneMat)==0;
model.genes(I)=[];
model.rxnGeneMat(:,I)=[];
if isfield(model,'geneShortNames')
    model.geneShortNames(I)=[];
end
if isfield(model,'geneMiriams')
    model.geneMiriams(I)=[];
end
if isfield(model,'geneFrom')
    model.geneFrom(I)=[];
end
if isfield(model,'geneComps')
    model.geneComps(I)=[];
end

%At this stage the model will contain some exchange reactions but probably
%not all (and maybe zero). This can be inconvenient, so all exchange
%reactions from the reference model are added, except for those which
%involve metabolites that are not in the model.

%First delete and included exchange reactions in order to prevent the order
%from changing
model=removeReactions(model,getExchangeRxns(model));

%Create a model with only the exchange reactions in refModel
excModel=removeReactions(refModel,setdiff(refModel.rxns,getExchangeRxns(refModel)),true,true);

%Find the metabolites there which are not exchange metabolites and which do
%not exist in the output model
I=~ismember(excModel.mets,model.mets) & excModel.unconstrained==0;

%Then find those reactions and delete them
[~, J]=find(excModel.S(I,:));
excModel=removeReactions(excModel,J,true,true);

%Merge with the output model
model=mergeModels({model;excModel},'metNames');
model.id='INITModel';
model.name=['Automatically generated model for ' tissue];
if any(celltype)
    model.name=[model.name ' - ' celltype];
end

if printReport==true
    printScores(model,'Final model statistics',hpaData,arrayData,tissue,celltype);
end

%Add information about essential reactions and reactions included for
%gap-filling and return a taskReport
if ~isempty(taskStructure)
    I=find(taskReport.ok); %Ignore failed tasks
    for i=1:numel(I)
        taskReport.essential{I(i),1}=cModel.rxns(essentialRxnMat(:,I(i)));
        taskReport.gapfill{I(i),1}=refModelNoExc.rxns(addedRxnMat(:,i));
    end
else
    taskReport=[];
end

%Fix grRules and reconstruct rxnGeneMat
[grRules,rxnGeneMat] = standardizeGrRules(model,true);
model.grRules = grRules;
model.rxnGeneMat = rxnGeneMat;
end

%This is for printing a summary of a model
function [rxnS, geneS]=printScores(model,name,hpaData,arrayData,tissue,celltype)
[a, b]=scoreModel(model,hpaData,arrayData,tissue,celltype);
rxnS=mean(a);
geneS=mean(b(~isinf(b)));
fprintf([name ':\n']);
fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']);
fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']);
fprintf(['\tMean gene score: ' num2str(geneS) '\n']);
fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']);
end

author: SysBioChalmers

repo: RAVEN


text:
function val=polygamma(k,z)
%%POLYGAMMA Evaluate the polygamma function. This function is the same as
%           the psi function except it will work with negative values of z.
%           The polygamma function is the value of the kth derivative of
%           the derivative of gamma(z) divided by gamma(z). That is, it is
%           the kth derivative of the digamma function.
%
%INPUTS: This function can be used with 1 or two inputs. If two inputs are
%        given, then it is polygamma(k,z). If only one is given, then it is
%        polygamma(z) where k is taken to be zero. The inputs are:
%        k The number of derivatives to take, k>=0. The digamma function is
%          k=0.
%        z A matrix of points at which the polygamma function is to be
%          evaluated.
%
%OUTPUTS: val The values of the polygamma function. This is the same size
%             as z.
%
%For values of z>0, this function just calls psi(k,z). For values of z<0,
%this function uses the identity of [1] to deal with negative values. For
%example, for k=0, one uses the identity:
%psi(1-x)=psi(x)+pi/tan(pi*x)
%
%REFERENCES:
%[1] Weisstein, Eric W. "Polygamma Function." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/PolygammaFunction.html
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin==1)
    z=k;
    k=0;
end

numVals=numel(z);
val=zeros(size(z));

for curVal=1:numVals
    x=z(curVal);
    if(x>0)
        val(curVal)=psi(k,x);
    else
        x=1-x;
        val(curVal)=(-1)^k*(psi(k,x)+pi^(k+1)*cotDerivVal(pi*x,k));
    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 varargout = drawCurve(varargin)
%DRAWCURVE draw a curve specified by a list of points
%
%   drawCurve(COORD);
%   packs coordinates in a single [N*2] array.
%
%   drawCurve(PX, PY);
%   specifies coordinates in separate arrays. PX and PY must be column
%   vectors with the same length.
%
%   drawCurve(..., TYPE);
%   where TYPE is either 'closed' or 'open', specifies if last point must
%   be connected to the first one ('closed') or not ('open').
%   Default is 'open'.
%
%   drawCurve(..., PARAM, VALUE);
%   specify plot options as described for plot command.
%
%   H = drawCurve(...) also return a handle to the list of line objects.
%
%   Example:
%   % Draw a curve representing an ellipse
%   t = linspace(0, 2*pi, 100)';
%   px = 10*cos(t); py = 5*sin(t);
%   drawCurve([px py], 'closed');
%   axis equal;
%
%   % The same, with different drawing options
%   drawCurve([px py], 'closed', 'lineWidth', 2, 'lineStyle', '--');
%
%   See Also:
%   polygons2d, drawPolyline
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 06/04/2004.
%

%   HISTORY
%   03/01/2007 better processing of input, and update doc (drawing
%       options and CLOSE option)
%   03/03/2010 add deprecation warning

% deprecation warning
warning('geom2d:deprecated', ...
    '''drawCurve'' is deprecated, use ''drawPolyline'' instead');

% default values
closed = false;

% If first argument is a cell array, draw each curve individually,
% and eventually returns handle of each plot.
var = varargin{1};
if iscell(var)
    h = [];
    for i=1:length(var(:))
        h = [h ; drawCurve(var{i}, varargin{2:end})]; %#ok<AGROW>
    end
    if nargout>0
        varargout{1}=h;
    end
    return;
end

% extract curve coordinate
if size(var, 2)==1
    % first argument contains x coord, second argument contains y coord
    px = var;
    if length(varargin)==1
        error('Wrong number of arguments in drawCurve');
    end
    py = varargin{2};
    varargin = varargin(3:end);
else
    % first argument contains both coordinate
    px = var(:, 1);
    py = var(:, 2);
    varargin = varargin(2:end);
end

% check if curve is closed or open
if ~isempty(varargin)
    var = varargin{1};
    if strncmpi(var, 'close', 5)
        closed = true;
        varargin = varargin(2:end);
    elseif strncmpi(var, 'open', 4)
        closed = false;
        varargin = varargin(2:end);
    end
end

% add first point at the end to close the curve
if closed
    px = [px; px(1)];
    py = [py; py(1)];
end

% plot the curve, with eventually optional parameters
h = plot(px, py, varargin{:});

% format output arguments
if nargout>0
    varargout{1}=h;
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [ objectIndex ] = GetScenario_waypointCurve(varargin)
% This scenario is designed to present a waypoint following exercise.

fprintf('[SCENARIO]\tGetting the waypoint following exercise.\n');

% DEFAULT CONFIGURATION
defaultConfig = struct('file','scenario.mat',...
                       'agents',[],...
                       'agentVelocity',0,...
                       'noiseFactor',0,...
                       'waypoints',3,...
                       'waypointRadius',0.1,...
                       'plot',0);

% Instanciate the scenario builder
SBinstance = scenarioBuilder();

% PARSE THE USER OVERRIDES USING THE SCENARIO BUILDER
[inputConfig] = SBinstance.configurationParser(defaultConfig,varargin);
% AGENT CONDITIONS
agentIndex = inputConfig.agents;
agentNumber = numel(agentIndex);                    % Declare the number of agents

%% DEFINE THE AGENT CONFIGURATION
% MOVE THROUGH THE AGENTS AND INITIALISE WITH GLOBAL PROPERTIES
fprintf('[SCENARIO]\tAssigning agent definition...\n');
for i=1:agentNumber
    agentIndex{i}.SetGLOBAL('position',[0;0;0] + inputConfig.noiseFactor*randn(3,1));
    agentIndex{i}.SetGLOBAL('velocity',[inputConfig.agentVelocity;0;0] + inputConfig.noiseFactor*randn(3,1));
    agentIndex{i}.SetGLOBAL('quaternion',[1;0;0;0]);
end

%% DEFINE THE WAYPOINT CONFIGURATION
fprintf('[SCENARIO]\tBuilding the new scenario...\n');
waypointDefiningRadius = 20;
angleOffset = -pi/2;                                % Align the first waypoint to be directly infront of the agent
waypointConfig = SBinstance.planarAngle(...
    'objects',agentNumber,...
    'radius',waypointDefiningRadius,...
    'pointA',[waypointDefiningRadius;-waypointDefiningRadius;-1],...
    'pointB',[waypointDefiningRadius;-waypointDefiningRadius;0],...
    'velocities',0,...
    'zeroAngle',angleOffset);

% MOVE THROUGH THE WAYPOINTS AND INITIALISE WITH GLOBAL PROPERTIES
fprintf('[SCENARIO]\tAssigning waypoint definitions:\n');
for index = 1:inputConfig.waypoints
    nameString = sprintf('WP-%s',agentIndex{1}.name);
    waypointIndex{index} = waypoint('radius',inputConfig.waypointRadius,'name',nameString);
    % APPLY GLOBAL STATE VARIABLES
    waypointIndex{index}.SetGLOBAL('position',waypointConfig.positions(:,index));
    waypointIndex{index}.SetGLOBAL('velocity',waypointConfig.velocities(:,index));
    waypointIndex{index}.SetGLOBAL('quaternion',waypointConfig.quaternions(:,index));
    % Assign waypoint to agent with priority
    waypointIndex{index} = waypointIndex{index}.CreateAgentAssociation(agentIndex{1},1/index);  % Create waypoint with association to agent
end

% BUILD THE COLLECTIVE OBJECT INDEX
objectIndex = horzcat(agentIndex,waypointIndex);
% PLOT THE SCENE
if inputConfig.plot
    SBinstance.plotObjectIndex(objectIndex);
end
clearvars -except objectIndex % Clean-up
end

repo: OpenMAS


text:
% Instantiate simulation from modified knowledge base.
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Author: Bonny Jain, bonny.jain@mit.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 12/20/2012
classdef Repressilator < edu.stanford.covert.cell.sim.runners.SimulationRunner
    methods
        function this = Repressilator(varargin)
            this = this@edu.stanford.covert.cell.sim.runners.SimulationRunner(varargin{:});
        end
    end

    methods (Access = protected)
        function modifyNetworkStructure(this, kb)
            import edu.stanford.covert.cell.kb.Gene;
            import edu.stanford.covert.cell.kb.ProteinMonomer;
            import edu.stanford.covert.cell.kb.ProteinComplex;
            import edu.stanford.covert.cell.kb.Stimuli;
            import edu.stanford.covert.cell.kb.TranscriptionUnit;

            %% create new genes, transcription units, etc.
            meanHL = mean([kb.mRNAGenes.halfLife]);
            expression = zeros(1, 3); %normal, cold shock, heat shock
            geneA = Gene(kb, NaN, {'GeneA'}, {'Gene A'}, {'genA'}, {''}, {'mRNA'}, 0, {''}, ...
                580500, 300, true, {''}, meanHL, expression);
            geneB = Gene(kb, NaN, {'GeneB'}, {'Gene B'}, {'genB'}, {''}, {'mRNA'}, 0, {''}, ...
                581000, 300, true, {''}, meanHL, expression);
            geneC = Gene(kb, NaN, {'GeneC'}, {'Gene C'}, {'genC'}, {''}, {'mRNA'}, 0, {''}, ...
                581500, 300, true, {''}, meanHL, expression);

            tuA = TranscriptionUnit(kb, NaN, {'TuA'}, {'Transcription unit A'},...
                -35, 6, ...
                -10, 6, ...
                -1);
            tuB = TranscriptionUnit(kb, NaN, {'TuB'}, {'Transcription unit B'},...
                -35, 6, ...
                -10, 6, ...
                -1);
            tuC = TranscriptionUnit(kb, NaN, {'TuC'}, {'Transcription unit C'},...
                -35, 6, ...
                -10, 6, ...
                -1);

            monA = ProteinMonomer(kb, NaN, {'MonomerA'}, {'Protein monomer A'},...
                {''}, {''}, {''}, ...
                10, {'dsDNA'}, {'dsDNA'}, ...
                {''}, {''}, {''}, {''}, {''}, {''}, ...
                false, {''}, {''}, 0, ...
                {''});
            monB = ProteinMonomer(kb, NaN, {'MonomerB'}, {'Protein monomer B'},...
                {''}, {''}, {''}, ...
                10, {'dsDNA'}, {'dsDNA'}, ...
                {''}, {''}, {''}, {''}, {''}, {''}, ...
                false, {''}, {''}, 0, ...
                {''});
            monC = ProteinMonomer(kb, NaN, {'MonomerC'}, {'Protein monomer C'},...
                {''}, {''}, {''}, ...
                10, {'dsDNA'}, {'dsDNA'}, ...
                {''}, {''}, {''}, {''}, {''}, {''}, ...
                false, {''}, {''}, 0, ...
                {''});

            cpxAA = ProteinComplex(kb, NaN, {'ComplexAA'}, {'Macromolecular complex AA'},...
                10, {'dsDNA'}, {'dsDNA'}, ...
                {''}, {''}, {''}, {''}, {''}, {''}, ...
                {'!StimulusA'}, {''});

            cpxBB = ProteinComplex(kb, NaN, {'ComplexBB'}, {'Macromolecular complex BB'},...
                10, {'dsDNA'}, {'dsDNA'}, ...
                {''}, {''}, {''}, {''}, {''}, {''}, ...
                {'!StimulusB'}, {''});

            cpxCC = ProteinComplex(kb, NaN, {'ComplexCC'}, {'Macromolecular complex CC'},...
                10, {'dsDNA'}, {'dsDNA'}, ...
                {''}, {''}, {''}, {''}, {''}, {''}, ...
                {''}, {''});

            stimA = Stimuli(kb, NaN, {'StimulusA'}, {'Stimulus A'});
            stimB = Stimuli(kb, NaN, {'StimulusB'}, {'Stimulus B'});

            this.createDefaultDisplacementReactions(kb, [cpxAA; cpxBB; cpxCC]);

            acgt = 'ACGT';
            genome = kb.genome;
            randSeq = ceil(rand(1, 582500-genome.sequenceLength) * 4);
            genome.sequence = [genome.sequence  acgt(randSeq)];
            genome.sequence(geneA.startCoordinate : geneA.endCoordinate) = acgt(ceil(rand(1, geneA.endCoordinate - geneA.startCoordinate + 1) * 4));
            genome.sequence(geneB.startCoordinate : geneB.endCoordinate) = acgt(ceil(rand(1, geneB.endCoordinate - geneB.startCoordinate + 1) * 4));
            genome.sequence(geneC.startCoordinate : geneC.endCoordinate) = acgt(ceil(rand(1, geneC.endCoordinate - geneC.startCoordinate + 1) * 4));
            genome.sequence(geneA.startCoordinate+(0:2)) = 'ATG';
            genome.sequence(geneB.startCoordinate+(0:2)) = 'ATG';
            genome.sequence(geneC.startCoordinate+(0:2)) = 'ATG';
            genome.sequence(geneA.endCoordinate+(-2:0)) = 'TAA';
            genome.sequence(geneB.endCoordinate+(-2:0)) = 'TAA';
            genome.sequence(geneC.endCoordinate+(-2:0)) = 'TAA';

            cytosol = kb.compartments(kb.cytosolCompartmentIndexs);
            extracellularSpace = kb.compartments(kb.extracellularCompartmentIndexs);

            %% create new relationships among new genes, transcription units, etc.
            geneA.genome = genome;
            geneA.transcriptionUnits = tuA;
            geneA.proteinMonomers = monA;
            geneA.compartment = cytosol;

            geneB.genome = genome;
            geneB.transcriptionUnits = tuB;
            geneB.proteinMonomers = monB;
            geneB.compartment = cytosol;

            geneC.genome = genome;
            geneC.transcriptionUnits = tuC;
            geneC.proteinMonomers = monC;
            geneC.compartment = cytosol;

            tuA.genome = genome;
            tuA.genes = geneA;
            tuA.geneCompartments = cytosol;
            tuA.compartment = cytosol;
            tuA.transcriptionFactorProteinComplexs = cpxBB;
            tuA.transcriptionFactorProteinComplexAffinitys = NaN; %nM
            tuA.transcriptionFactorProteinComplexActivitys = 1/10;
            tuA.transcriptionFactorProteinComplexConditions = {''};
            tuA.transcriptionFactorProteinComplexBindingSiteStartCoordinates = 580460;
            tuA.transcriptionFactorProteinComplexBindingSiteLengths = 10;
            tuA.transcriptionFactorProteinComplexBindingSiteDirections = true;
            tuA.transcriptionFactorProteinComplexCompartments = cytosol;

            tuB.genome = genome;
            tuB.genes = geneB;
            tuB.geneCompartments = cytosol;
            tuB.compartment = cytosol;
            tuB.transcriptionFactorProteinComplexs = cpxCC;
            tuB.transcriptionFactorProteinComplexAffinitys = NaN; %nM
            tuB.transcriptionFactorProteinComplexActivitys = 1/10;
            tuB.transcriptionFactorProteinComplexConditions = {''};
            tuB.transcriptionFactorProteinComplexBindingSiteStartCoordinates = 580960;
            tuB.transcriptionFactorProteinComplexBindingSiteLengths = 10;
            tuB.transcriptionFactorProteinComplexBindingSiteDirections = true;
            tuB.transcriptionFactorProteinComplexCompartments = cytosol;

            tuC.genome = genome;
            tuC.genes = geneC;
            tuC.geneCompartments = cytosol;
            tuC.compartment = cytosol;
            tuC.transcriptionFactorProteinComplexs = cpxAA;
            tuC.transcriptionFactorProteinComplexAffinitys = NaN; %nM
            tuC.transcriptionFactorProteinComplexActivitys = 1/10;
            tuC.transcriptionFactorProteinComplexConditions = {''};
            tuC.transcriptionFactorProteinComplexBindingSiteStartCoordinates = 581460;
            tuC.transcriptionFactorProteinComplexBindingSiteLengths = 10;
            tuC.transcriptionFactorProteinComplexBindingSiteDirections = true;
            tuC.transcriptionFactorProteinComplexCompartments = cytosol;

            monA.gene = geneA;
            monA.geneCompartments = cytosol;
            monA.compartment = cytosol;

            monB.gene = geneB;
            monB.geneCompartments = cytosol;
            monB.compartment = cytosol;

            monC.gene = geneC;
            monC.geneCompartments = cytosol;
            monC.compartment = cytosol;

            cpxAA.proteinMonomers = monA;
            cpxAA.proteinMonomerCompartments = cytosol;
            cpxAA.proteinMonomerCoefficients = 2;
            cpxAA.regulatedTranscriptionUnits = tuC;
            cpxAA.compartment = cytosol;
            cpxAA.complexFormationProcess = kb.processes.findobj('wholeCellModelID', 'Process_MacromolecularComplexation');
            cpxAA.stimuliRegulators = stimA;

            cpxBB.proteinMonomers = monB;
            cpxBB.proteinMonomerCompartments = cytosol;
            cpxBB.proteinMonomerCoefficients = 2;
            cpxBB.regulatedTranscriptionUnits = tuA;
            cpxBB.compartment = cytosol;
            cpxBB.complexFormationProcess = kb.processes.findobj('wholeCellModelID', 'Process_MacromolecularComplexation');
            cpxBB.stimuliRegulators = stimB;

            cpxCC.proteinMonomers = monC;
            cpxCC.proteinMonomerCompartments = cytosol;
            cpxCC.proteinMonomerCoefficients = 2;
            cpxCC.regulatedTranscriptionUnits = tuB;
            cpxCC.compartment = cytosol;
            cpxCC.complexFormationProcess = kb.processes.findobj('wholeCellModelID', 'Process_MacromolecularComplexation');

            %specificy dynamics of stimuli piece-wise
            stimA.compartments = [extracellularSpace; extracellularSpace];
            stimA.values = [true; false];
            stimA.initialTimes = [0; 10];
            stimA.finalTimes = [10; inf];
            stimA.regulatedProteinComplexs = cpxAA;

            stimB.compartments = extracellularSpace;
            stimB.values = true;
            stimB.initialTimes = 0;
            stimB.finalTimes = inf;
            stimB.regulatedProteinComplexs = cpxBB;

            genome.genes = [genome.genes; geneA; geneB; geneC];
            genome.transcriptionUnits = [genome.transcriptionUnits; tuA; tuB; tuC];

            kb.genes = genome.genes;
            kb.transcriptionUnits = genome.transcriptionUnits;
            kb.proteinMonomers = [kb.proteinMonomers; monA; monB; monC];
            kb.proteinComplexs = [kb.proteinComplexs; cpxAA; cpxBB; cpxCC];
            kb.stimulis = [kb.stimulis; stimA; stimB];

            %% class super class method
            this.modifyNetworkStructure@edu.stanford.covert.cell.sim.runners.SimulationRunner(kb);
        end

        function modifyNetworkParameters(~, sim)
            %% get handles
            g = sim.gene;
            time = sim.state('Time');
            rna = sim.state('Rna');
            trn = sim.process('Transcription');

            %% get constants
            nascentMRNAIndexs = find(any(rna.nascentRNAGeneComposition(g.mRNAIndexs, :), 1));
            [~, modTuIndexs] = ismember({'TuA'; 'TuB'; 'TuC'}, rna.wholeCellModelIDs(rna.nascentIndexs));

            %% get parameter values
            tuBindProb = trn.transcriptionUnitBindingProbabilities;
            meanBindProb = mean(tuBindProb(setdiff(nascentMRNAIndexs, modTuIndexs)));
            rnaDecayRates = rna.decayRates(rna.matureIndexs);

            %% calculate modified parameter values
            %set binding probability of new TUs
            modTuBindProb = tuBindProb;
            modTuBindProb(modTuIndexs) = meanBindProb;

            %renormalize
            modTuBindProb(nascentMRNAIndexs) = modTuBindProb(nascentMRNAIndexs) * sum(tuBindProb(nascentMRNAIndexs)) / sum(modTuBindProb(nascentMRNAIndexs));

            %update expression
            modRnaExp = (rna.nascentRNAMatureRNAComposition * modTuBindProb) ./ (log(2) / time.cellCycleLength + rnaDecayRates);
            modRnaExp = modRnaExp / sum(modRnaExp);

            %% update parameter values
            trn.transcriptionUnitBindingProbabilities = modTuBindProb;
            rna.expression(rna.matureIndexs) = modRnaExp;
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell