text:
clear all; close all; clc;
%% Info Generative Adversarial Network
%% Settings
args.maxepochs = 5; args.c_weight = 1; args.z_dim = 128;
args.batch_size = 16; args.image_size = [64,64,3];
args.lrD = 0.0002; args.lrG = 0.001; args.beta1 = 0.5;
args.beta2 = 0.999; args.cc_dim = 1; args.dc_dim = 10;
args.sample_size = 100;
%% Weights, Biases, Offsets and Scales
%% Generator
[paramsGen,stGen] = initializeGen(args);
%% Discriminator
[paramsDis,stDis] = initializeDis(args);
%% Train
out = false; epoch = 0; global_iter = 0; count = 1;
% average Gradient and average Gradient squared holders
avgG.Dis = []; avgGS.Dis = []; avgG.Gen = []; avgGS.Gen = [];

while ~out
    tic; epoch = epoch+1;
    shuffleid = randperm(100000);
    fprintf('Epoch %d\n',epoch)
    for i=1:floor(length(shuffleid)/args.batch_size)
        global_iter = global_iter+1;
        idx = (i-1)*args.batch_size+1:i*args.batch_size;
        XBatch=miniBatch(shuffleid(idx));

        [GradDis,GradGen,stDis,stGen] = ...
                dlfeval(@modelGradients,XBatch,genNoiseCDC(args),...
                paramsDis,paramsGen,args,stDis,stGen);

        % Update Discriminator network parameters
        [paramsDis,avgG.Dis,avgGS.Dis] = ...
            adamupdate(paramsDis, GradDis, ...
            avgG.Dis, avgGS.Dis, global_iter, ...
            args.lrD, args.beta1, args.beta2);

        % Update Generator network parameters
        [paramsGen,avgG.Gen,avgGS.Gen] = ...
            adamupdate(paramsGen, GradGen, ...
            avgG.Gen, avgGS.Gen, global_iter, ...
            args.lrG, args.beta1, args.beta2);

        if i==1 || rem(i,20) == 0
            % progress plot
            progressplot(args,paramsGen,stGen)
        end

        if rem(global_iter,count^2)==0 || rem(global_iter,500)==0
            h = gcf;
            % Capture the plot as an image
            frame = getframe(h);
            im = frame2im(frame);
            [imind,cm] = rgb2ind(im,256);
            % Write to the GIF File
            if count == 1
              imwrite(imind,cm,'InfoGANcelebA.gif','gif', 'Loopcount',inf);
            else
              imwrite(imind,cm,'InfoGANcelebA.gif','gif','WriteMode','append');
            end

            count = count+1;
        end

    end

    elapsedTime = toc;
    disp("Epoch "+epoch+". Time taken for epoch = "+elapsedTime + "s")

    if epoch == args.maxepochs
        out = true;
    end
end

%% Help Functions
%% miniBatch
function dlx = miniBatch(ids)
path = 'D:\44754\Documents\Data\celeba-dataset\img_align_celeba\img_align_celeba\';
filesA = dir([path '*.jpg']);
dlx = [];
for i = 1:length(ids)
    imj=imresize(imread([path filesA(ids(i)).name]),[64,64]);
    imj=preprocess(imj);
    dlx = cat(4,dlx,im2double(imj));
end
dlx = gpdl(dlx);
end
%% preprocess
function x = preprocess(x)
x = double(x)./255;
x = (x-.5)./.5;
x = reshape(x,64,64,3,[]);
end
%% Weight Initialization
function parameter = initializeGaussian(parameterSize)
parameter = randn(parameterSize, 'single') .* 0.01;
end
%% extract data
function x = gatext(x)
x = gather(extractdata(x));
end
%% gpu dl array wrapper
function dlx = gpdl(x,labels)
if nargin < 2
    labels='SSCB';
end
dlx = gpuArray(dlarray(single(x),labels));
end
%% Generator
function [dly,st] = Generator(dlx,params,st)
dly = dlx;
exstatbatch='[dly,st] = batchnormwrap(dly,params,st,%d);';
exstattran = 'dly = dltranspconv(dly,params.TCW%d,params.TCb%d,"Stride",%d,"Cropping",%d);';
exstatrelu = 'dly = relu(dly);';

for i = 1:5
    if i==1
        eval(sprintf(exstattran,i,i,1,0))
    elseif i==5
        eval(sprintf(exstattran,i,i,2,1))
        dly = tanh(dly);
    else
        eval(sprintf(exstattran,i,i,2,1))
        eval(sprintf(exstatbatch,i-1))
        eval(exstatrelu)
    end
end
end
%% Discriminator
function [dly,st] = Discriminator(dlx,params,args,st)
dly = dlx;
exstatconv ='dly = dlconv(dly,params.CNW%d,params.CNb%d,"Stride",%d,"Padding",%d);';
exstatleak = 'dly = leakyrelu(dly,.1);';
exstatbatch='[dly,st] = batchnormwrap(dly,params,st,%d);';

for i = 1:5
    if i == 1
        eval(sprintf(exstatconv,i,i,2,1))
        eval(exstatleak)
    elseif i == 5
        eval(sprintf(exstatconv,i,i,1,0))
    else
        eval(sprintf(exstatconv,i,i,2,1))
        eval(sprintf(exstatbatch,i-1))
        eval(exstatleak)
    end
end
dly = gpdl(reshape(dly,1+args.cc_dim+args.dc_dim,[]),'CB');
dly(1,:) = sigmoid(dly(1,:));
dly(args.cc_dim+2:end,:) = softmax(dly(args.cc_dim+2:end,:));

end
%% General Noise
function out = genNoiseCDC(args)
noi = randn([args.z_dim,args.batch_size]);
out = [noi;gen_cc([args.cc_dim,args.batch_size]);...
    gen_dc([args.dc_dim,args.batch_size])];
out = reshape(out,1,1,args.z_dim+args.cc_dim+args.dc_dim,[]);
out = gpuArray(dlarray(out,'SSCB'));
end
%% Generate Discrete Code
function out = gen_dc(shape)
out = zeros(shape);
random_cate = randi([1,shape(1)],shape(2),1);
for i = 1:length(random_cate)
    out(random_cate(i),i) = 1;
end
end
%% Generate Continuous Code
function out = gen_cc(shape)
out = randn(shape)*0.5;
end
%% Progress Plot
function progressplot(args,paramsGen,stGen)
fixednoise = zeros(args.z_dim,args.sample_size);
tmp = zeros(args.cc_dim,args.sample_size);
for i = 1:10
    tmp(1,(i-1)*10+1:i*10) = linspace(-2,2,10);
end
cc = tmp;

tmp = zeros(args.dc_dim,args.sample_size);
for i = 1:10
    tmp(i,(i-1)*10+1:i*10) = 1;
end
dc = tmp;

fake_data = gpdl(reshape(cat(1,fixednoise,cc,dc),1,1,...
    args.z_dim+args.cc_dim+args.dc_dim,[]));
fake_images = extractdata(Generator(fake_data,paramsGen,stGen));

fig = gcf;
if ~isempty(fig.Children)
    delete(fig.Children)
end

I = imtile(fake_images);
I = rescale(I);
imagesc(I)
title("Generated Images")

drawnow;

end
%% Report Progress
function [d_loss,g_loss] = reportprogress(x,z,paramsDis,...
    paramsGen,args,stDis,stGen)
fake_images = Generator(z,paramsGen,stGen);
d_output_real = Discriminator(x,paramsDis,args,stDis);
d_output_fake = Discriminator(fake_images,paramsDis,args,stDis);

% Loss due to true or not
d_loss_a = -mean(log(d_output_real(1,:))+log(1-d_output_fake(1,:)));
g_loss_a = -mean(log(d_output_fake(1,:)));

% cc loss
output_cc = d_output_fake(2,:);
d_loss_cc = mean((output_cc/0.5).^2);

% softmax classification loss
output_dc = d_output_fake(3:end,:);
d_loss_dc = -(mean(sum(z(args.z_dim+args.cc_dim+1:end,:).*output_dc,1))+...
    mean(sum(z(args.z_dim+args.cc_dim+1:end,:).*z(args.z_dim+args.cc_dim+1:end,:),1)));

% Discriminator Loss
d_loss = d_loss_a+args.c_weight*d_loss_cc+d_loss_dc;
% Generator Loss
g_loss = g_loss_a+args.c_weight*d_loss_cc+d_loss_dc;
end
%% Model Gradients
function [GradDis,GradGen,stDis,stGen] = modelGradients(x,z,paramsDis,...
    paramsGen,args,stDis,stGen)
[fake_images,stGen] = Generator(z,paramsGen,stGen);
d_output_real = Discriminator(x,paramsDis,args,stDis);
[d_output_fake,stDis] = Discriminator(fake_images,paramsDis,args,stDis);

% Loss due to true or not
d_loss_a = -mean(log(d_output_real(1,:))+log(1-d_output_fake(1,:)));
g_loss_a = -mean(log(d_output_fake(1,:)));

% cc loss
output_cc = d_output_fake(2,:);
d_loss_cc = mean((output_cc/0.5).^2);

% softmax classification loss
output_dc = d_output_fake(3:end,:);
z = gpdl(reshape(z,args.z_dim+args.cc_dim+args.dc_dim,[]),'CB');
d_loss_dc = -(mean(sum(z(args.z_dim+args.cc_dim+1:end,:).*output_dc,1))+...
    mean(sum(z(args.z_dim+args.cc_dim+1:end,:).*z(args.z_dim+args.cc_dim+1:end,:),1)));

% Discriminator Loss
d_loss = d_loss_a+args.c_weight*d_loss_cc+d_loss_dc;
% Generator Loss
g_loss = g_loss_a+args.c_weight*d_loss_cc+d_loss_dc;

% For each network, calculate the gradients with respect to the loss.
GradGen = dlgradient(g_loss,paramsGen,'RetainData',true);
GradDis = dlgradient(d_loss,paramsDis);

end
%% Initialize Generator params
function [paramsGen,stGen] = initializeGen(args)
paramsGen.TCW1 = dlarray(initializeGaussian([4,4,1024,...
    args.z_dim+args.cc_dim+args.dc_dim]));
paramsGen.TCb1 = dlarray(zeros(1024,1,'single'));
paramsGen.TCW2 = dlarray(initializeGaussian([4,4,512,1024]));
paramsGen.TCb2 = dlarray(zeros(512,1,'single'));
paramsGen.BNo1 = dlarray(zeros(512,1,'single'));
paramsGen.BNs1 = dlarray(ones(512,1,'single'));
paramsGen.TCW3 = dlarray(initializeGaussian([4,4,256,512]));
paramsGen.TCb3 = dlarray(zeros(256,1,'single'));
paramsGen.BNo2 = dlarray(zeros(256,1,'single'));
paramsGen.BNs2 = dlarray(ones(256,1,'single'));
paramsGen.TCW4 = dlarray(initializeGaussian([4,4,128,256]));
paramsGen.TCb4 = dlarray(zeros(128,1,'single'));
paramsGen.BNo3 = dlarray(zeros(128,1,'single'));
paramsGen.BNs3 = dlarray(ones(128,1,'single'));
paramsGen.TCW5 = dlarray(initializeGaussian([4,4,3,128]));
paramsGen.TCb5 = dlarray(zeros(3,1,'single'));

% States for Batch Norm
stGen.BN1 = []; stGen.BN2 = []; stGen.BN3 = [];
end
%% Initialize Discriminator params
function [paramsDis,stDis] = initializeDis(args)
paramsDis.CNW1 = dlarray(initializeGaussian([4,4,3,128]));
paramsDis.CNb1 = dlarray(zeros(128,1,'single'));
paramsDis.CNW2 = dlarray(initializeGaussian([4,4,128,256]));
paramsDis.CNb2 = dlarray(zeros(256,1,'single'));
paramsDis.BNo1 = dlarray(zeros(256,1,'single'));
paramsDis.BNs1 = dlarray(ones(256,1,'single'));
paramsDis.CNW3 = dlarray(initializeGaussian([4,4,256,512]));
paramsDis.CNb3 = dlarray(zeros(512,1,'single'));
paramsDis.BNo2 = dlarray(zeros(512,1,'single'));
paramsDis.BNs2 = dlarray(ones(512,1,'single'));
paramsDis.CNW4 = dlarray(initializeGaussian([4,4,512,1024]));
paramsDis.CNb4 = dlarray(zeros(1024,1,'single'));
paramsDis.BNo3 = dlarray(zeros(1024,1,'single'));
paramsDis.BNs3 = dlarray(ones(1024,1,'single'));
paramsDis.CNW5 = dlarray(initializeGaussian([4,4,1024,...
    1+args.cc_dim+args.dc_dim]));
paramsDis.CNb5 = dlarray(zeros(1+args.cc_dim+args.dc_dim,1,'single'));

% States for Batch Norm
stDis.BN1 = []; stDis.BN2 = []; stDis.BN3 = [];
end
%% batchnormwrap
function [dly,st] = batchnormwrap(dlx,params,st,num)
exstat1=sprintf('if isempty(st.BN%d),',num);
exstat2=sprintf('[dly,st.BN%d.mu,st.BN%d.sig]=batchnorm(dlx,params.BNo%d,params.BNs%d,"MeanDecay",0.8);else,',num*ones(1,4));
exstat3=sprintf('[dly,st.BN%d.mu,st.BN%d.sig]=batchnorm(dlx,params.BNo%d,params.BNs%d,st.BN%d.mu,st.BN%d.sig,"MeanDecay",0.8);end',num*ones(1,6));
eval(strcat(exstat1,exstat2,exstat3));
end

author: zcemycl

repo: Matlab-GAN

save_path: github-repos/MATLAB/zcemycl-Matlab-GAN


text:
function f = divergence(F)
%DIVERGENCE   Divergence of a CHEBFUN3V object.
%   DIVERGENCE(F) returns divergence of the CHEBFUN3V object F as a
%   CHEBFUN3. If F = U i + V j + W k, then divergence(F) = U_x + V_y + W_z.
%
% See also CHEBFUN3V/DIV.

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

% Empty check:
if ( isempty(F) )
    f = chebfun3();
    return
end

% Two components:
if ( F.nComponents == 2 )
    Fc = F.components;
    f = diff(Fc{1}, 1, 1) + diff(Fc{2}, 1, 2);

% Three components:
elseif ( F.nComponents == 3 )
    Fc = F.components;
    diff1 = diff(Fc{1}, 1, 1);
    diff2 = diff(Fc{2}, 1, 2);
    diff3 = diff(Fc{3}, 1, 3);
    %Developer Note: Instead of f = diff1 + diff2 + diff3; which calls the
    % constructor two times, we use the following trick to call it just
    % once. See CHEBFUN3/PLUS for more details:
    vscales = [vscale(diff1) + vscale(diff2), vscale(diff3)];

    m = 51; % size of sampling grid
    LVals = sample(diff1, m, m, m) + sample(diff2, m, m, m) + ...
        sample(diff3, m, m, m);
    LVscale = max(abs(LVals(:)));
    if sum(vscales) == 0
        f = chebfun3(0);
    else
        kappa = sum(vscales)/LVscale;
        pref = chebfunpref().cheb3Prefs;
        eps = pref.chebfun3eps;
        tol = eps*kappa;
        f = chebfun3(@(x,y,z) feval(diff1, x, y, z) + feval(diff2, x, y, z) + ...
           feval(diff3, x, y, z) , Fc{1}.domain, 'eps', tol);
    end
else
    error('CHEBFUN:CHEBFUN3V:divergence:notSupported', ...
        'Two or three components are needed.')

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function s = nnz( x )

%   Disciplined convex/geometric programming information for NNZ:
%      When used in CVX models, NNZ cannot deduce the numeric value of a
%      CVX variable, so it returns the number of "structurally" nonzero
%      elements; i.e., those that are not *identically* zero. This is
%      similar to the concept of structural nonzeros in sparse matrix
%      factorizations. To illustrate the distinction, consider:
%         variable x(3);
%         x(2) == 0;
%         y = [x(1);0;x(3)];
%      In this case, NNZ(x) returns 3, even though x(2) has been set to
%      zero in an equality constraint. (After all, the overall model may
%      be infeasible, in which case x(2) is arguably not zero.) However,
%      NNZ(y) returns 2, because y(2) is identically zero.

s = nnz( any( x.basis_, 1 ) );

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

author: yu-jiang

repo: radpbook

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


text:
function [stat] = ft_spike_xcorr(cfg, spike)

% FT_SPIKE_XCORR computes the cross-correlation histogram and shift predictor.
%
% Use as
%   [stat] = ft_spike_xcorr(cfg, data)
%
% The input SPIKE should be organised as the spike or the raw datatype, obtained from
% FT_SPIKE_MAKETRIALS or FT_PREPROCESSING (in that case, conversion is done
% within the function). A mex file is located in /contrib/spike/private
% which will be automatically mexed.
%
% Configurations options for xcorr general:
%   cfg.maxlag           = number in seconds, indicating the maximum lag for the
%                          cross-correlation function in sec (default = 0.1 sec).
%   cfg.debias           = 'yes' (default) or 'no'. If 'yes', we scale the
%                          cross-correlogram by M/(M-abs(lags)), where M = 2*N -1 with N
%                          the length of the data segment.
%   cfg.method           = 'xcorr' or 'shiftpredictor'. If 'shiftpredictor'
%                           we do not compute the normal cross-correlation
%                           but shuffle the subsequent trials.
%                           If two channels are independent, then the shift
%                           predictor should give the same correlogram as the raw
%                           correlogram calculated from the same trials.
%                           Typically, the shift predictor is subtracted from the
%                           correlogram.
%   cfg.outputunit       = - 'proportion' (value in each bin indicates proportion of occurence)
%                          - 'center' (values are scaled to center value which is set to 1)
%                          - 'raw' (default) unnormalized crosscorrelogram.
%   cfg.binsize          = [binsize] in sec (default = 0.001 sec).
%   cfg.channelcmb       = Mx2 cell-array with selection of channel pairs (default = {'all' 'all'}),
%                          see FT_CHANNELCOMBINATION for details
%   cfg.latency          = [begin end] in seconds, 'max' (default), 'min', 'prestim'(t<=0), or
%                          'poststim' (t>=0).%
%   cfg.vartriallen      = 'yes' (default) or 'no'.
%                          If 'yes' - accept variable trial lengths and use all available trials
%                          and the samples in every trial.
%                          If 'no'  - only select those trials that fully cover the window as
%                          specified by cfg.latency and discard those trials that do not.
%                          if cfg.method = 'yes', then cfg.vartriallen
%                          should be 'no' (otherwise, fewer coincidences
%                          will occur because of non-overlapping windows)
%   cfg.trials           = numeric selection of trials (default = 'all')
%   cfg.keeptrials       = 'yes' or 'no' (default)
%
% A peak at a negative lag for stat.xcorr(chan1,chan2,:) means that chan1 is leading
% chan2. Thus, a negative lag represents a spike in the second dimension of
% stat.xcorr before the channel in the third dimension of stat.stat.
%
% Variable trial length is controlled by the option cfg.vartriallen. If it is
% specified as cfg.vartriallen='yes', all trials are selected that have a minimum
% overlap with the latency window of cfg.maxlag. However, the shift predictor
% calculation demands that following trials have the same amount of data, otherwise,
% it does not control for rate non-stationarities. If cfg.vartriallen = 'yes', all
% trials should fall in the latency window, otherwise we do not compute the shift
% predictor.
%
% Output:
%    stat.xcorr            = nchans-by-nchans-by-(2*nlags+1) cross correlation histogram with dimord 'chan_chan_time'
%   or
%    stat.shiftpredictor   = nchans-by-nchans-by-(2*nlags+1) shift predictor with dimord 'chan_chan_time'
% and
%    stat.lags             = (2*nlags + 1) vector with lags in seconds.
%    stat.trial            = ntrials-by-nchans-by-nchans-by-(2*nlags + 1) with single trials and dimord 'rpt_chan_chan_time'
%    stat.label            = corresponding labels to channels in stat.xcorr
%    stat.cfg              = configurations used in this function

% Copyright (C) 2010-2012, Martin Vinck
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble provenance spike

% check input spike structure
spike = ft_checkdata(spike, 'datatype', 'spike', 'feedback', 'yes');

% get the default options
cfg.trials         = ft_getopt(cfg, 'trials', 'all');
cfg.latency        = ft_getopt(cfg, 'latency','maxperiod');
cfg.keeptrials     = ft_getopt(cfg, 'keeptrials', 'no');
cfg.method         = ft_getopt(cfg, 'method', 'xcorr');
cfg.channelcmb     = ft_getopt(cfg, 'channelcmb', 'all');
cfg.vartriallen    = ft_getopt(cfg, 'vartriallen', 'yes');
cfg.debias         = ft_getopt(cfg, 'debias', 'yes');
cfg.maxlag         = ft_getopt(cfg, 'maxlag', 0.01);
cfg.binsize        = ft_getopt(cfg, 'binsize', 0.001);
cfg.outputunit     = ft_getopt(cfg, 'outputunit', 'proportion');

% ensure that the options are valid
cfg = ft_checkopt(cfg, 'latency', {'char', 'ascendingdoublebivector'});
cfg = ft_checkopt(cfg, 'trials', {'char', 'doublevector', 'logical'});
cfg = ft_checkopt(cfg, 'keeptrials', 'char', {'yes', 'no'});
cfg = ft_checkopt(cfg, 'method', 'char', {'xcorr', 'shiftpredictor'});
cfg = ft_checkopt(cfg, 'channelcmb', {'char', 'cell'});
cfg = ft_checkopt(cfg, 'vartriallen', 'char', {'no', 'yes'});
cfg = ft_checkopt(cfg, 'debias', 'char', {'yes', 'no'});
cfg = ft_checkopt(cfg, 'maxlag', 'doublescalar');
cfg = ft_checkopt(cfg, 'binsize', 'doublescalar');
cfg = ft_checkopt(cfg, 'outputunit', 'char', {'proportion', 'center', 'raw'});

cfg = ft_checkconfig(cfg, 'allowed', {'latency', 'trials', 'keeptrials', 'method', 'channelcmb', 'vartriallen', 'debias', 'maxlag', 'binsize', 'outputunit'});

doShiftPredictor  = strcmp(cfg.method, 'shiftpredictor'); % shift predictor

% determine the corresponding indices of the requested channel combinations
cfg.channelcmb = ft_channelcombination(cfg.channelcmb, spike.label(:), true);
cmbindx        = zeros(size(cfg.channelcmb));
for k=1:size(cfg.channelcmb,1)
  cmbindx(k,1) = strmatch(cfg.channelcmb(k,1), spike.label, 'exact');
  cmbindx(k,2) = strmatch(cfg.channelcmb(k,2), spike.label, 'exact');
end
nCmbs 	   = size(cmbindx,1);
chansel    = unique(cmbindx(:)); % get the unique channels
nChans     = length(chansel);
if isempty(chansel), error('No channel was selected'); end

% get the number of trials or change SPIKE according to cfg.trials
nTrials = size(spike.trialtime,1);
if  strcmp(cfg.trials,'all')
  cfg.trials = 1:nTrials;
elseif islogical(cfg.trials) || all(cfg.trials==0 | cfg.trials==1)
  cfg.trials = find(cfg.trials);
end
cfg.trials = sort(cfg.trials(:));
if max(cfg.trials)>nTrials, error('maximum trial number in cfg.trials should not exceed length of DATA.trial'); end
if isempty(cfg.trials), errors('No trials were selected in cfg.trials'); end
nTrials = length(cfg.trials);

% get the latencies
begTrialLatency = spike.trialtime(cfg.trials,1);
endTrialLatency = spike.trialtime(cfg.trials,2);
trialDur = endTrialLatency - begTrialLatency;

% select the latencies
if strcmp(cfg.latency,'minperiod')
  cfg.latency = [max(begTrialLatency) min(endTrialLatency)];
elseif strcmp(cfg.latency,'maxperiod')
  cfg.latency = [min(begTrialLatency) max(endTrialLatency)];
elseif strcmp(cfg.latency,'prestim')
  cfg.latency = [min(begTrialLatency) 0];
elseif strcmp(cfg.latency,'poststim')
  cfg.latency = [0 max(endTrialLatency)];
end

% check whether the time window fits with the data
if (cfg.latency(1) < min(begTrialLatency)), cfg.latency(1) = min(begTrialLatency);
  warning('Correcting begin latency of averaging window');
end
if (cfg.latency(2) > max(endTrialLatency)), cfg.latency(2) = max(endTrialLatency);
  warning('Correcting begin latency of averaging window');
end

% get the maximum number of lags in samples
fsample = (1/cfg.binsize);
if cfg.maxlag>(cfg.latency(2)-cfg.latency(1))
  warning('Correcting cfg.maxlag since it exceeds latency period');
  cfg.maxlag = (cfg.latency(2) - cfg.latency(1));
end
nLags        = round(cfg.maxlag*fsample);
lags         = (-nLags:nLags)/fsample; % create the lag axis

% check which trials will be used based on the latency
fullDur       = trialDur>=(cfg.maxlag); % only trials which are larger than maximum lag
overlaps      = endTrialLatency>(cfg.latency(1)+cfg.maxlag) & begTrialLatency<(cfg.latency(2)-cfg.maxlag);
hasWindow     = ones(nTrials,1);
if strcmp(cfg.vartriallen,'no') % only select trials that fully cover our latency window
  startsLater    = single(begTrialLatency) > (single(cfg.latency(1)) + 0.5*cfg.binsize);
  endsEarlier    = single(endTrialLatency) < (single(cfg.latency(2)) - 0.5*cfg.binsize);
  hasWindow      = ~(startsLater | endsEarlier); % check this in all other funcs
end
trialSel           = fullDur(:) & overlaps(:) & hasWindow(:);
cfg.trials         = cfg.trials(trialSel);
begTrialLatency    = begTrialLatency(trialSel);
endTrialLatency    = endTrialLatency(trialSel);

if isempty(cfg.trials), warning('No trials were selected'); end
if length(cfg.trials)<2&&doShiftPredictor
  error('Shift predictor can only be calculated with more than 1 selected trial.');
end

% if we allow variable trial lengths
if strcmp(cfg.vartriallen,'yes') && doShiftPredictor && ~strcmp(cfg.outputunit,'proportion')
  warning('using cfg.vartriallen = "yes" and shift predictor method: please use cfg.outputunit = "proportion"');
end
nTrials   = length(cfg.trials); % only now reset nTrials

% preallocate the sum and the single trials and the shift predictor
keepTrials   = strcmp(cfg.keeptrials,'yes');
s     = zeros(nChans,nChans,2*nLags);
if keepTrials
  warning('storing single trials for cross correlation is memory expensive, please check');
  singleTrials = zeros(nTrials,nChans,nChans,2*nLags);
end

if strcmp(cfg.method,'shiftpredictor'), singleTrials(1,:,:,:) = NaN; end

if ((cfg.latency(2)-cfg.latency(1))/cfg.maxlag)<2.5
  warning('selected latency will cause highly variable xcorr at borders');
end

ft_progress('init', 'text',     'Please wait...');
for iTrial = 1:nTrials
  origTrial = cfg.trials(iTrial);
  ft_progress(iTrial/nTrials, 'Processing trial %d from %d', iTrial, nTrials);

  for iCmb = 1:nCmbs
    % take only the times that are in this trial and latency window
    indx  = cmbindx(iCmb,:);
    inTrial1 = spike.trial{indx(1)}==origTrial;
    inTrial2 = spike.trial{indx(2)}==origTrial;
    inWindow1 = spike.time{indx(1)}>=cfg.latency(1) &  spike.time{indx(1)}<=cfg.latency(2);
    inWindow2 = spike.time{indx(2)}>=cfg.latency(1) &  spike.time{indx(2)}<=cfg.latency(2);
    ts1 = sort(spike.time{indx(1)}(inTrial1(:) & inWindow1(:)));
    ts2 = sort(spike.time{indx(2)}(inTrial2(:) & inWindow2(:)));

    switch cfg.method
      case 'xcorr'
        % compute the xcorr if both are non-empty
        if ~isempty(ts1) && ~isempty(ts2)
          if indx(1)<=indx(2)
            [x]   = spike_crossx_matlab(ts1(:),ts2(:),cfg.binsize,nLags*2+1);  % removed the mex file for now, until issue on windows is resolved
          else
            [x]   = spike_crossx_matlab(ts2(:),ts1(:),cfg.binsize,nLags*2+1); % removed the mex file for now, until issue on windows is resolved
          end

          % remove the center peaks from the auto-correlogram
          if indx(1)==indx(2), x(nLags:nLags+1) = 0; end

          if strcmp(cfg.debias,'yes')
            lags = (-nLags:nLags)*cfg.binsize;
            lags = (lags(2:end)+lags(1:end-1))/2;
            T    = nanmin([endTrialLatency(iTrial);cfg.latency(2)])-nanmax([begTrialLatency(iTrial);cfg.latency(1)]);
            sc = (T./(T-abs(lags(:))));
            sc = length(sc)*sc./sum(sc);
            x  = x(:).*sc(:);
          end
          % sum the xcorr
          s(indx(1),indx(2),:) = s(indx(1),indx(2),:) + shiftdim(x(:),-2);
          s(indx(2),indx(1),:) = s(indx(2),indx(1),:) + shiftdim(flipud(x(:)),-2);

          % store individual trials if requested
          if keepTrials
            singleTrials(iTrial,indx(1),indx(2),:) = shiftdim(x(:),-3);
            singleTrials(iTrial,indx(2),indx(1),:) = shiftdim(flipud(x(:)),-3);
          end
        end
      case 'shiftpredictor'
        if iTrial>1
          % symmetric, get x21 to x12 and x22 to x11
          inTrial1_old = spike.trial{indx(1)}==cfg.trials(iTrial-1);
          ts1_old      = sort(spike.time{indx(1)}(inTrial1_old(:) & inWindow1(:)));
          inTrial2_old = spike.trial{indx(2)}==cfg.trials(iTrial-1);
          ts2_old      = sort(spike.time{indx(2)}(inTrial2_old(:) & inWindow2(:)));

          % compute both combinations
          for k = 1:2

            % take chan from this and previous channel
            if k==1
              A = ts1;
              B = ts2_old;
            else
              A = ts1_old;
              B = ts2;
            end
            if ~isempty(A) && ~isempty(B)
              if indx(1)<=indx(2)
                [x]   = spike_crossx_matlab(A(:),B(:),cfg.binsize,nLags*2+1);
              else
                [x]   = spike_crossx_matlab(B(:),A(:),cfg.binsize,nLags*2+1);
              end
              % remove the center peaks from the auto-correlogram
              if indx(1)==indx(2), x(nLags:nLags+1) = 0; end

              if strcmp(cfg.debias,'yes')
                lags = (-nLags:nLags)*cfg.binsize;
                lags = (lags(2:end)+lags(1:end-1))/2;
                T    = nanmin([endTrialLatency(iTrial);cfg.latency(2)])-nanmax([begTrialLatency(iTrial);cfg.latency(1)]);
                sc = (T./(T-abs(lags(:))));
                sc = length(sc)*sc./sum(sc);
                x  = x(:).*sc(:);
              end
              % compute the sum
              s(indx(1),indx(2),:) =  s(indx(1),indx(2),:) + shiftdim(x(:),-2);
              s(indx(2),indx(1),:) =  s(indx(2),indx(1),:) + shiftdim(flipud(x(:)),-2);
              if keepTrials
                singleTrials(iTrial,indx(1),indx(2),:) = shiftdim(x(:)/2,-3);
                singleTrials(iTrial,indx(2),indx(1),:) = shiftdim(flipud(x(:))/2,-3);
              end
            end
          end
        end
    end % symmetric shift predictor loop
  end % combinations
end % trial loop
ft_progress('close')
% multiply the shift sum by a factor so it has the same scale as raw
dofShiftPred = 2*(nTrials-1);
if doShiftPredictor, s = s*nTrials/dofShiftPred; end
switch cfg.outputunit
  case 'proportion'
    sm = repmat(nansum(s,3),[1 1 nLags*2]);
    s = s./sm;
    if keepTrials
      singleTrials       = singleTrials./repmat(shiftdim(sm,-1),[nTrials 1 1 1]);
    end
  case 'center'
    center   = repmat(s(:,:,nLags+1),[1 1 nLags*2]);
    s        = s./center;
    if keepTrials
      singleTrials       = singleTrials./repmat(shiftdim(center,-1),[nTrials 1 1 1]);
    end
end

% return the results
stat.(cfg.method)      = s;
lags = (-nLags:nLags)*cfg.binsize;
lags = (lags(2:end)+lags(1:end-1))/2;
stat.time              = lags;
stat.dimord            = 'chan_chan_time';
if keepTrials
  stat.trial   = singleTrials;
  stat.dimord  = 'trial_chan_chan_time';
end
stat.label       = spike.label(chansel);

% do the general cleanup and bookkeeping at the end of the function

ft_postamble previous   spike
ft_postamble provenance stat
ft_postamble history    stat

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [C] = spike_crossx_matlab(tX,tY,binsize,nbins)
tX = sort(tX(:));
tY = sort(tY(:));

minLag = - binsize * (nbins-1) / 2;
j = 0:nbins-1;
B = minLag + j * binsize;
tX(tX<(tY(1)+minLag) | tX>(tY(end)-minLag))   = [];
if isempty(tX)
  C = zeros(1,length(B)-1);
  return;
end
tY(tY>(tX(end)-minLag) | tY<(tX(1)+minLag)) = [];
if isempty(tY)
  C = zeros(1,length(B)-1);
  return;
end
nX = length(tX); nY = length(tY);

% compute all distances at once using a multiplication trick
if (nX*nY)<2*10^7  % allow matrix to grow to about 150 MB, should always work
  D = log(exp(-tX(:))*exp(tY(:)'));
  D = D(:);
  D(abs(D)>abs(minLag)) = [];
  [C] = histc(D,B);
  C(end) = [];
else
  % break it down in pieces such that nX*nY<2*10*7
  k   = 2;
  nXs = round(nX/k);
  while (nXs*nY)>=(2*10^7)
    k = k+1;
    nXs = round(nX/k);
  end

  % get the indices
  steps = round(nX/k);
  begs = 1:steps:steps*k;
  ends = begs+(steps-1);
  rm   = begs>nX;
  begs(rm) = [];
  ends(rm) = [];
  ends(end) = nX;
  nSteps    = length(begs);

  C = zeros(1,length(B));
  for iStep = 1:nSteps
    d = log(exp(-tX(begs(iStep):ends(iStep)))*exp(tY(:)'));
    d = d(:);
    d(abs(d)>abs(minLag)) = [];
    C = C + histc(d,B)';
  end
  C(end) = [];
end

if isempty(C)
  C = zeros(1,length(B)-1);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef IN1KMOP9 < EvoXBenchProblem
% <multi> <real> <large/none> <expensive/none>

%------------------------------- Reference --------------------------------
% Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture
% search as multiobjective optimization benchmarks: Problem formulation and
% performance assessment, IEEE Transactions on Evolutionary Computation,
% 2023.
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            config.name = 'mnv3';
            config.args.objs = 'err&params&flops&latency';
            config.args.normalized_objectives = false;
            obj.Setting@EvoXBenchProblem(config);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function Z = separate_freq(Y)
	if iscell(Y)
		Z = {};

		for m = 0:length(Y)-1
			Z{m+1} = separate_freq(Y{m+1});
		end

		return;
	end

	Z.signal = {};
	Z.meta = struct();

	r = 1;
	for k = 1:length(Y.signal)
		if ~iscell(Y.signal{k})
			% because matlab is stupid and ignores last dim if 1
			sz_orig = size(Y.signal{k});
			if numel(sz_orig) == 2
				sz_orig(3) = 1;
			end
		else
			sz_orig = [length(Y.signal{k}) 0 size(Y.signal{k}{1},2)];
		end

		for l = 1:sz_orig(1)
			if ~iscell(Y.signal{k})
				nsignal = Y.signal{k}(l,:);
			else
				nsignal = Y.signal{k}{l};
			end

			nsignal = reshape(nsignal,[numel(nsignal)/sz_orig(3) sz_orig(3)]);

			Z.signal{r} = nsignal;

			r = r+1;
		end
	end

	[temp,I] = sortrows(Y.meta.j(1:size(Y.meta.j,1),:).');

	Z.signal = Z.signal(I);

	field_names = fieldnames(Y.meta);
	for n = 1:length(field_names)
		src_value = getfield(Y.meta,field_names{n});
		Z.meta = setfield(Z.meta,field_names{n},src_value(:,I));
	end

	Z.signal = cellfun(@(x)(permute(x,[1 3 2])), Z.signal, ...
		'UniformOutput', false);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function curr_layer = F_weight2activation(curr_layer)

curr_layer.a = curr_layer.W;

end

author: singaxiong

repo: SignalGraph


text:
function n2 = dist2(x, c)
%DIST2	Calculates squared distance between two sets of points.
%
%	Description
%	D = DIST2(X, C) takes two matrices of vectors and calculates the
%	squared Euclidean distance between them.  Both matrices must be of
%	the same column dimension.  If X has M rows and N columns, and C has
%	L rows and N columns, then the result has M rows and L columns.  The
%	I, Jth entry is the  squared distance from the Ith row of X to the
%	Jth row of C.
%
%	See also
%	GMMACTIV, KMEANS, RBFFWD
%

%	Copyright (c) Ian T Nabney (1996-2001)

if nargin<2
    c = x;
end
[ndata, dimx] = size(x);
[ncentres, dimc] = size(c);
if dimx ~= dimc
	error('Data dimension does not match dimension of centres')
end

n2 = (ones(ncentres, 1) * sum((x.^2)', 1))' + ...
  ones(ndata, 1) * sum((c.^2)',1) - ...
  2.*(x*(c'));

% Rounding errors occasionally cause negative entries in n2
if any(any(n2<0))
  n2(n2<0) = 0;
end

author: BatzoglouLabSU

repo: SIMLR


text:
function overlayMapsSameSession(view);
% overlayMapsSameSession(view);
%
% Take two maps (co, amp, ph, or map) from the current
% view, and superimpose them in different color channels,
% with an interface for adjusting the thresholds of each/
% moving across slices if necessary.
%
%
% 06/04 ras: wrote it.

% dialog to get params
prompt = {...
        'map 1 field (can be co, amp, ph, or map):',...
        'map 1 clip (2 numbers, min and max--from 0 to 1 normalized):',...
        'map 1 scan:',...
        'map 1 data type #:',...
        'map 2 field (can be co, amp, ph, or map):',...
        'map 2 clip (2 numbers, min and max--from 0 to 1 normalized):',...
        'map 2 scan:',...
        'map 2 data type #:',...
        'Add Anatomical as underlay? (1 for yes, 0 for no):',...
        'Anat Clip (if 1 above):',...
    };
dlgTitle = 'Overlay Two maps from the same session....';
lineNo = 1;
defaults = {...
        viewGet(view,'displayMode'),...
        '0 1',...
        num2str(getCurScan(view)),...
        num2str(viewGet(view,'curdt')),...
        viewGet(view,'displayMode'),...
        '0 1',...
        '',...
        num2str(viewGet(view,'curdt')),...
        '1',...
        num2str(viewGet(view,'anatClip')),...
    };
AddOpts.Resize = 'on';
AddOpts.WindowStyle='normal';
AddOpts.Interpreter='tex';
answer = inputdlg(prompt,dlgTitle,lineNo,defaults,AddOpts);

% if cancel pressed, exit gracefull
if isempty(answer)
    return
end

% parse the responses
map1Name = answer{1};
map1Clip = str2num(answer{2});
map1Scan = str2num(answer{3});
map1Dt = str2num(answer{4});
map2Name = answer{5};
map2Clip = str2num(answer{6});
map2Scan = str2num(answer{7});
map2Dt = str2num(answer{8});
anatFlag = str2num(answer{9});
anatClip = str2num(answer{10});

curdt = viewGet(view,'curdt');
if map1Dt ~= curdt
    if map1Name=='map'
        fprintf('Whoops ... haven''t added the ability to get a map from another dt.\n');
        return
    end

    view = viewSet(view,'curdt',map1Dt);
    view = loadCorAnal(view);
end
field1 = viewGet(view,map1Name);
vol1 = field1{map1Scan};

curdt2 = viewGet(view,'curdt');
if map2Dt ~= curdt2
    if map2Name=='map'
        fprintf('Whoops ... haven''t added the ability to get a map from another dt.\n');
        return
    end

    view = viewSet(view,'curdt',map2Dt);
    view = loadCorAnal(view);
end
field2 = viewGet(view,map2Name);
vol2 = field2{map2Scan};

hbox = msgbox('Getting maps to overlay...');

% % apply clip
% for i = 1:size(vol1,3)
%     clip = map1Clip .* [min(vol1(:)) max(vol1(:))];
%     vol1(:,:,i) = rescale2(vol1(:,:,i),clip,clip);
%
%     clip = map2Clip .* [min(vol2(:)) max(vol2(:))];
%     vol2(:,:,i) = rescale2(vol2(:,:,i),clip,clip);
% end

name1 = sprintf('%s Scan %i',map1Name,map1Scan);
name2 = sprintf('%s Scan %i',map2Name,map2Scan);

% use a separte function to add the
% background image, do the overlaying:
overlayMaps(view,vol1,vol2,anatFlag,name1,name2);

close(hbox);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function t_runmarket(quiet)
%T_RUNMARKET  Tests for code in RUNMKT, SMARTMKT AND AUCTION.

%   MATPOWER
%   Copyright (c) 2005-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

n_tests = 20;

t_begin(n_tests, quiet);

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;

if ~have_feature('smartmarket')
    t_skip(n_tests, 'smartmarket code not available');
else
    mpc = loadcase('t_auction_case');

    mpopt = mpoption('opf.ac.solver', 'MIPS', 'out.lim.all', 1, 'out.branch', 0, 'out.sys_sum', 0, 'out.all', 0, 'verbose', 1);
    % mpopt = mpoption('out.gen', 1, 'out.branch', 0, 'out.sys_sum', 0);

    offers.P.qty = [
        12 24 24;
        12 24 24;
        12 24 24;
        12 24 24;
        12 24 24;
        12 24 24;
    ];
    offers.P.prc = [
        20 50 60;
        20 40 70;
        20 42 80;
        20 44 90;
        20 46 75;
        20 48 60;
    ];
    bids.P.qty = [
        10 10 10;
        10 10 10;
        10 10 10;
    ];
    bids.P.prc = [
        100 70 60;
%         100 64.3 20;
%         100 30.64545 0;
        100 50 20;
        100 60 50;
    ];

    offers.Q.qty = [ 60; 60; 60; 60; 60; 60; 0; 0; 0 ];
    offers.Q.prc = [ 0; 0; 0; 0; 0; 3; 0; 0; 0 ];
    bids.Q.qty = [ 15; 15; 15; 15; 15; 15; 15; 12; 7.5 ];
%     bids.Q.prc = [ 0; 0; 0; 0; 0; 0; 0; 83.9056; 0 ];
    bids.Q.prc = [ 0; 0; 0; 0; 0; 0; 0; 20; 0 ];

    t = 'marginal Q offer, marginal PQ bid, auction_type = 5';
    mkt = struct(   'auction_type', 5, ...
                    't' , [], ...
                    'u0', [], ...
                    'lim', []   );
    [r, co, cb, f, dispatch, success, et] = runmarket(mpc, offers, bids, mkt, mpopt);
    co5 = co;
    cb5 = cb;

%     [ co.P.qty co.P.prc ]
%     [ cb.P.qty cb.P.prc ]
%     [ co.Q.qty co.Q.prc ]
%     [ cb.Q.qty cb.Q.prc ]

    i2e = r.bus(:, BUS_I);
    e2i = sparse(max(i2e), 1);
    e2i(i2e) = (1:size(r.bus, 1))';
    G = find( ~isload(r.gen) );   %% real generators
    L = find(  isload(r.gen) );   %% dispatchable loads
    Gbus = e2i(r.gen(G,GEN_BUS));
    Lbus = e2i(r.gen(L,GEN_BUS));

    t_is( co.P.qty, ones(6, 1) * [12 24 0], 2, [t ' : gen P quantities'] );
    t_is( co.P.prc(1,:), 50.1578*ones(1,3), 3, [t ' : gen 1 P prices'] );
    t_is( cb.P.qty, [10 10 10; 10 0.196 0; 10 10 0], 2, [t ' : load P quantities'] );
    t_is( cb.P.prc(2,:), 56.9853*ones(1,3), 4, [t ' : load 2 P price'] );
    t_is( co.P.prc(:,1), r.bus(Gbus, LAM_P), 8, [t ' : gen P prices'] );
    t_is( cb.P.prc(:,1), r.bus(Lbus, LAM_P), 8, [t ' : load P prices'] );

    t_is( co.Q.qty, [4.2722; 11.3723; 14.1472; 22.8939; 36.7886; 12.3375; 0; 0; 0], 2, [t ' : Q offer quantities'] );
    t_is( co.Q.prc, [0;0;0;0;0;3; 0.4861; 2.5367; 1.3763], 4, [t ' : Q offer prices'] );
    t_is( cb.Q.qty, [0;0;0;0;0;0; 15; 4.0785; 5], 2, [t ' : Q bid quantities'] );
    t_is( cb.Q.prc, [0;0;0;0;0;3; 0.4861; 2.5367; 1.3763], 4, [t ' : Q bid prices'] );
    t_is( co.Q.prc, r.bus([Gbus; Lbus], LAM_Q), 8, [t ' : Q offer prices'] );
    t_is( cb.Q.prc, co.Q.prc, 8, [t ' : Q bid prices'] );

    t = 'marginal Q offer, marginal PQ bid, auction_type = 0';
    mkt.auction_type = 0;
    [r, co, cb, f, dispatch, success, et] = runmarket(mpc, offers, bids, mkt, mpopt);
    t_is( co.P.qty, co5.P.qty, 8, [t ' : gen P quantities'] );
    t_is( cb.P.qty, cb5.P.qty, 8, [t ' : load P quantities'] );
    t_is( co.P.prc, offers.P.prc, 8, [t ' : gen P prices'] );
    t_is( cb.P.prc, bids.P.prc, 8, [t ' : load P prices'] );

    t_is( co.Q.qty, co5.Q.qty, 8, [t ' : gen Q quantities'] );
    t_is( cb.Q.qty, cb5.Q.qty, 8, [t ' : load Q quantities'] );
    t_is( co.Q.prc, offers.Q.prc, 8, [t ' : gen Q prices'] );
    t_is( cb.Q.prc, bids.Q.prc, 8, [t ' : load Q prices'] );
end

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function V_corr = make_vcorr(DD,pv,nb,nl,f,Zb)
%MAKE_VCORR  Voltage Correction used in distribution power flow
%
%   V_corr = make_vcorr(DD,pv,nb,nl,f,Zb)
%
%   Calculates voltage corrections with current generators placed at PV
%   buses. Their currents are calculated with the voltage difference at PV
%   buses break points and loop impedances. The slack bus voltage is set to
%   zero. Details can be seen in
%   D. Rajicic, R. Ackovski and R. Taleski, "Voltage correction power flow,"
%   IEEE Transactions on Power Delivery, vol. 9, no. 2, pp. 1056-1062, Apr 1994.
%   https://doi.org/10.1109/61.296308
%
%   See also RADIAL_PF.

V_corr = zeros(nb,1);
I = zeros(nb,1);
I(pv) = DD;
% backward sweep
for k = nl:-1:2
    i = f(k);
    I(i) = I(i) + I(k);
end
% forward sweep
for k = 2:nl
    i = f(k);
    V_corr(k) = V_corr(i) - Zb(k) * I(k);
end


text:
function a = gproc(hyper)

%=============================================================================
% GPROC Gaussian-Process object
%=============================================================================
%
% A=GPROC(hyper) returns a GPROC object initialized with hyperparameters hyper.
% minimize log-likelihood for length runs
%
%
%  Hyperparameters, and their defaults
%
%
%   child=kernel;        -- the kernel is stored as a member called "child"
%   length=50            -- default maximum iteration of line searches.
%
%  Model
%  H			-- covariance hyperparameters
%
%  Methods:
%  d=gen(toyreg({'o=1','n=2','l=200'}))
%
%  [r,a]=train(gproc,get(d,1:100))
%  r=test(a,get(d,101:200))
%  % Note that r.X has 2 an extra output columns!
%
%   train, test
%=============================================================================
% Reference : Introduction to gaussian processes
% Author    : David J.C. MacKay
% Link      : http://www.inference.phy.cam.ac.uk/mackay/gpB.ps.gz
%=============================================================================

  % hyperparams
  a.input= [];
  a.target = [];
  a.length=50;
  a.H=[];
  a.pred_var=0;
  % model

  p=algorithm('gproc');
  a= class(a,'gproc',p);
  a.algorithm.use_signed_output=0;

  if nargin==1
    eval_hyper;
  end;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function r = ne(p,q)
% this function compares the inequality condition of x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'CellVector')&&isa(q, 'CellVector'))
    r=p;
    r.xvalue = p.xvalue~=q.xvalue;
    r.yvalue = p.yvalue~=q.yvalue;
    r.zvalue = p.zvalue~=q.zvalue;
elseif isa(p, 'CellVector')
    r=p;
    r.xvalue = p.xvalue~=q;
    r.yvalue = p.yvalue~=q;
    r.zvalue = p.zvalue~=q;
else
    r=q;
    r.xvalue = p~=q.xvalue;
    r.yvalue = p~=q.yvalue;
    r.zvalue = p~=q.zvalue;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [names, mf_settings, data, data_endo, data_endo_a, data_endo_c, data_endo_c_lags, data_exo, data_exo_a, data_exo_p, data_exo_c, data_exo_c_lags, Fperiods, Fcomp, Fcperiods, Fcenddate]= ...
    gensample_mf(startdate,enddate,VARtype,Fstartdate,Fenddate,Fendsmpl,endo,exo,frequency,lags,F,CF,pref, numendo)

% Phase 1: data loading and error checking

% first read the data from Excel
[data_m, names_m]=xlsread(pref.excelFile,'mfvar_monthly');
[data_q, names_q]=xlsread(pref.excelFile,'mfvar_quarterly');
[select, names_s] = xlsread(pref.excelFile,'mf_var_trans');

n_m = size(data_m,2);           mf_settings.Nm = n_m;
n_q = size(data_q,2);           mf_settings.Nq = n_q;
T_m = size(data_m,1);           mf_settings.T_m = T_m;
T_q = size(data_q,1);
T_b = T_q*3;                    mf_settings.T_b = T_b;    % The end of the balanced dataset (when all quarterly and monthly data exist)
ym_q = kron(data_q,ones(3,1));
% ym_q = reshape(repmat(data_q,1,3)',T_b,1);

data = [data_m [ym_q; nan(T_m - T_b,n_q)]];
data(:,select == 0) = log(data(:,select==0));
data(:,select == 1) = data(:,select==1)./100;
names =  [names_m repmat(names_m(:,end),1,n_q)];
names(1,n_m+2:end) = names_q(1,2:end);

% Check if the monthly variables
% if n_q ~= floor(n_m/3)

% now, as a preliminary step: check if there is any Nan in the data; if yes, return an error since the model won't be able to run with missing data
% a simple way to test for NaN is to check for "smaller or equal to infinity": Nan is the only number for which matlab will return 'false' when asked so
% [r,c]=size(data);
% for ii=1:r
%    for jj=1:c
%    temp=data(ii,jj);
%       if (temp<=inf)==0
%       % identify the variable and the date
%       NaNvariable=names{1,jj+1};
%       NaNdate=names{ii+1,1};
%       message=['Error: variable ' NaNvariable ' at date ' NaNdate ' (and possibly other sample entries) is identified as NaN. Please check your Excel spreadsheet: entry may be blank or non-numerical.'];
%       msgbox(message);
%       error('programme termination: data error');
%       end
%    end
% end

% identify the date strings
datestrings=names_m(2:end,1);
% identify the position of the string corresponding to the start period
startlocation=find(strcmp(datestrings,startdate));
% identify the position of the string corresponding to the end period
endlocation=find(strcmp(datestrings,enddate));

% if either the start date or the date date is not recognised, return an error message
if isempty(startlocation)
    error('bear:gensample_mf', 'Unknown start date for the sample. Please check your sample start date (remember that names are case-sensitive).')
elseif isempty(endlocation)
    error('bear:gensample_mf', 'Unknown end date for the sample. Please check your sample end date (remember that names are case-sensitive).');
end
% also, if the start date is posterior to the end date, obviously return an error
if startlocation>=endlocation==1
    error('bear:gensample_mf','Inconsistency between the start and end dates. The start date must be anterior to the end date.');
end

% identify the variable strings, endogenous and exogenous
variablestrings=[names_m(1,2:end) names_q(1,2:end)];

% identify the position of the strings corresponding to the endogenous variables
% for each variable, find the corresponding string
for ii=1:numendo
    % check first that the variable ii in endo appears in the list of variable strings
    % if not, the variable is unknown: return an error
    var=endo{ii,1};
    check=find(strcmp(variablestrings,var));
    if isempty(check)==1
        error('bear:gensample_mf', ['Endogenous variable ' var ' cannot be found on the excel data spreadsheet.']);
    end
    % if the variable is known, go on
    endolocation(ii,1)=find(strcmp(variablestrings,endo(ii,1)));
end

% identify the position of the strings corresponding to the exogenous variables
% proceed similarly to the endogenous variables, but account for the fact that exogenous may be empty
% so check first whether there are exogenous variables altogether
if isempty(exo)
    numexo=0;
else
    % if not empty, repeat what has been done with the exogenous
    numexo=size(exo,1);
    % for each variable, find the corresponding string
    for ii=1:numexo
        % check first that the variable ii in endo appears in the list of variable strings
        % if not, the variable is unknown: return an error
        var=exo{ii,1};
        check=find(strcmp(variablestrings,var));
        if isempty(check)==1
            error('bear:gensample_mf', ['Error: exogenous variable ' var ' cannot be found on the excel data spreadsheet.']);
        end
        % if the variable is known, go on
        exolocation(ii,1)=find(strcmp(variablestrings,exo(ii,1)));
    end
end

% Phase 2: creation of the data matrices data_endo and data_exo

% now create the matrix of endogenous variables for the estimation sample
% it is simply the concatenation of the vectors of each endogenous variables, over the selected sample dates
data_endo=[];
% loop over endogenous variables
for ii=1:numendo
    data_endo=[data_endo data(startlocation:endlocation,endolocation(ii,1))];
end

% Similarly, create the matrix of exogenous variables for the estimation sample
data_exo=[];
for ii=1:numexo
    data_exo=[data_exo data(startlocation:endlocation,exolocation(ii,1))];
end

% these are the transformations for the relevant variables; 0 for taking logs and 1 for dividing by 100
% They are read from the sheet mf_var_trans on line 23 in this file
mf_settings.select = select(:,endolocation(:,1)');

% Phase 3: determination of the position of the forecast start and end periods

% if both unconditional and conditional forecasts were not selected, there is no need for all the forecast-specific matrices: simply return empty matrices
if (VARtype==1 && F==0) || ((VARtype==2 || VARtype==3 || VARtype==5 || VARtype==6|| VARtype==7 ) && (F==0 && CF==0))
    data_endo_a=[];
    data_exo_a=[];
    data_exo_p=[];
    Fperiods=[];
    Fcomp=[];
    Fcperiods=[];
    data_endo_c=[];
    data_endo_c_lags=[];
    data_exo_c=[];
    data_exo_c_lags=[];
    Fcenddate=[];

    % if forecast were selected, create all the required elements
else

    % preliminary tasks
    % first, identify the date strings, and the variable strings
    datestrings=names(2:end,1);
    variablestrings=names(1,2:end);

    % identify the location of the last period in the dataset
    dataendlocation=size(datestrings,1);

    % identify the position of the start period (the easy part)
    % if the start period has been selected as the first period after the sample end, identifies it directly
    if Fendsmpl==1
        Fstartlocation=find(strcmp(datestrings,enddate))+1;
        % if the start period has not been selected as the first period after the sample end, it must be within the sample: look for it
    elseif Fendsmpl==0
        Fstartlocation=find(strcmp(datestrings,Fstartdate));
        % if the start date is not recognised, return an error message
        if isempty(Fstartlocation)==1
            error('bear:gensample_mf', 'Unknown start date for the forecasts. Select a date within a sample, or select "Start forecasts after last sample period"');
        end
    end

    % identify the position of the final forecast period (the hard part)
    % this period can be in or outside the sample, depending on the user's choice

    % if the data is yearly
    if frequency==1
        Fendlocation=str2num(Fenddate(1,1:end-1))-str2num(datestrings{1,1}(1,1:end-1))+1;

        % if the data is quarterly
    elseif frequency==2
        % first identify the year and quarter of the initial date in the whole data set (not just the sample)
        datastartyear=str2num(datestrings{1,1}(1,1:4));
        datastartquarter=str2num(datestrings{1,1}(1,6));
        % convert this date into quarters only
        datastart=datastartyear*4+datastartquarter;
        % then identify the year and quarter of the final forecast date
        forecastendyear=str2num(Fenddate(1,1:4));
        forecastendquarter=str2num(Fenddate(1,6));
        % convert this date into quarters only
        forecastend=forecastendyear*4+forecastendquarter;
        % finally, compute the number of periods that separate the two dates
        Fendlocation=forecastend-datastart+1;

        % if the data is monthly
    elseif frequency==3
        % first identify the year and month of the initial date in the whole data set (not just the sample)
        temp=datestrings{1,1};
        datastartyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,datastartmonth]=strtok(temp);
        % convert this date into months only
        datastart=datastartyear*12+str2num(datastartmonth);
        % then identify the year and month of the final forecast date
        temp=Fenddate;
        forecastendyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,forecastendmonth]=strtok(temp);
        % convert this date into months only
        forecastend=forecastendyear*12+str2num(forecastendmonth);
        Fendlocation=forecastend-datastart+1;

        % if the data is weekly
    elseif frequency==4
        % then identify the year and week corresponding to this final period
        temp=datestrings{end,1};
        dataendyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,dataendweek]=strtok(temp);
        dataendweek=str2num(dataendweek);
        % identify the year and week corresponding to the end of the forecast period
        temp=Fenddate;
        Fendyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,Fendweek]=strtok(temp);
        Fendweek=str2num(Fendweek);
        % determine whether the forecast period ends within the data set, or after the end of the data set
        if Fendyear<dataendyear
            insmpl=1;
        elseif Fendyear==dataendyear && Fendweek<=dataendweek
            insmpl=1;
        else
            insmpl=0;
        end
        % if the forecast end period lies within the dataset, simply detect its location
        if insmpl==1
            Fendlocation=find(strcmp(Fenddate,datestrings));
            % if it is outside the data set, complete until the forecast end date is reached (assuming 52 weeks per year)
        elseif insmpl==0
            % Consider two cases separately
            % first case: if the end year of the forecast is the same as the end year of the dataset
            % in this case, simply complete the missing periods
            if dataendyear==Fendyear
                complement=Fendweek-dataendweek;
                Fendlocation=dataendlocation+complement;
                % if the end year of the forecasts is posterior to the last year of the data
            elseif dataendyear<Fendyear
                % complete the first year (the one shared that ends the data set)
                complement=52-dataendweek;
                % complete the following years before the last one (if any)
                for ii=dataendyear+1:Fendyear-1
                    complement=complement+52;
                end
                % complete the final forecast year
                complement=complement+Fendweek;
                Fendlocation=dataendlocation+complement;
            end
        end

        % if the data is daily
    elseif frequency==5
        % then identify the year and day corresponding to this final period
        temp=datestrings{end,1};
        dataendyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,dataendday]=strtok(temp);
        dataendday=str2num(dataendday);
        % identify the year and day corresponding to the end of the forecast period
        temp=Fenddate;
        Fendyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,Fendday]=strtok(temp);
        Fendday=str2num(Fendday);
        % determine whether the forecast period ends within the data set, or after the end of the data set
        if Fendyear<dataendyear
            insmpl=1;
        elseif Fendyear==dataendyear && Fendday<=dataendday
            insmpl=1;
        else
            insmpl=0;
        end
        % if the forecast end period lies within the dataset, simply detect its location
        if insmpl==1
            Fendlocation=find(strcmp(Fenddate,datestrings));
            % if it is outside the data set, complete until the forecast end date is reached (assuming 261 opening days per year)
        elseif insmpl==0
            % Consider two cases separately
            % first case: if the end year of the forecast is the same as the end year of the dataset
            % in this case, simply complete the missing periods
            if dataendyear==Fendyear
                complement=Fendday-dataendday;
                Fendlocation=dataendlocation+complement;
                % if the end year of the forecasts is posterior to the last year of the data
            elseif dataendyear<Fendyear
                % complete the first year (the one shared that ends the data set)
                complement=261-dataendday;
                % complete the following years before the last one (if any)
                for ii=dataendyear+1:Fendyear-1
                    complement=complement+261;
                end
                % complete the final forecast year
                complement=complement+Fendday;
                Fendlocation=dataendlocation+complement;
            end
        end

        % finally, if the data is undated
    elseif frequency==6
        Fendlocation=str2num(Fenddate(1,1:end-1))-str2num(datestrings{1,1}(1,1:end-1))+1;
    end

    % from this, conclude the total number of forecast periods
    Fperiods=Fendlocation-Fstartlocation+1;

    if Fperiods<0
        error('bear:gensample_mf','The forecast start date needs to be prior to the forecast end date');
    end

    % Phase 4: generation of the forecast-specific matrices

    % now create the matrix of endogenous variables for the pre-forecast period
    % it is simply the concatenation of the vectors of each endogenous variables, over the selected sample dates
    data_endo_a=[];
    % loop over endogenous variables
    for ii=1:numendo
        data_endo_a=[data_endo_a data(1:Fstartlocation-1,endolocation(ii,1))];
    end
    % also, create the matrix of exogenous variables for the pre-forecast period
    data_exo_a=[];
    for ii=1:numexo
        data_exo_a=[data_exo_a data(1:Fstartlocation-1,exolocation(ii,1))];
    end

    % create the matrix of endogenous variables for the period common to actual data and forecasts (for forecast evaluation)
    % first, check that there are such common periods: it is the case if the beginning of the forecast period is anterior to the end of the dataset
    if Fstartlocation<=dataendlocation
        % return a scalar value to indicate that forecast evaluation is possible
        Fcomp=1;
        % compute the number of common periods
        % if the forecast period ends before the end of the data set, the common periods end with the end of the forecasts
        if Fendlocation<=dataendlocation
            Fcperiods=Fperiods;
            % record the end date of the common periods
            Fcenddate=Fenddate;
            % if the forecast period ends later than the data set, the common periods end at the end of the data set
        elseif Fendlocation>dataendlocation
            Fcperiods=dataendlocation-Fstartlocation+1;
            % record the end date of the common periods
            Fcenddate=datestrings{end,1};
        end

        % create a matrix of endogenous data for the common periods
        data_endo_c=[];
        for ii=1:numendo
            data_endo_c=[data_endo_c data(Fstartlocation:min(dataendlocation,Fendlocation),endolocation(ii,1))];
        end

        % create a lagged matrix of endogenous data prior to the common periods
        % the number of values is equal to "lags"; this will be used for computation of the log predictive score
        data_endo_c_lags=[];
        for ii=1:numendo
            data_endo_c_lags=[data_endo_c_lags data(Fstartlocation-lags:Fstartlocation-1,endolocation(ii,1))];
        end

        % create a matrix of exogenous data for the common periods
        data_exo_c=[];
        for ii=1:numexo
            data_exo_c=[data_exo_c data(Fstartlocation:min(dataendlocation,Fendlocation),exolocation(ii,1))];
        end

        % create a lagged matrix of exogenous data prior to the common periods
        % the number of values is equal to "lags"; this will be used for computation of the log predictive score
        data_exo_c_lags=[];
        for ii=1:numexo
            data_exo_c_lags=[data_exo_c_lags data(Fstartlocation-lags:Fstartlocation-1,exolocation(ii,1))];
        end
        % if there are no common periods, return a scalar value to indicate that forecast evaluation is not possible
    else
        Fcomp=0;
        Fcperiods=0;
        data_exo_c=[];
        data_endo_c=[];
        Fcenddate=[];
        data_endo_c_lags=[];
        data_exo_c_lags=[];
    end

    % now create the matrix data_exo_p
    % two possible cases

    % if there are no exogenous variables, simply create an empty matrix
    if isempty(exo)
        data_exo_p=[];

        % if there are exogenous variables, load from excel
    else
        % load the data from Excel
        [num txt strngs]=xlsread(pref.excelFile,'pred exo');

        % obtain the row location of the forecast start date
        [Fstartlocation,~]=find(strcmp(strngs,Fstartdate));
        % check that the start date for the forecast appears in the sheet; if not, return an error
        if isempty(Fstartlocation)
            message=['Error: a forecast application is selected for a model that uses exogenous variables. Hence, predicted exogenous values should be supplied over the forecast periods. Yet the start date for forecasts (' Fstartdate ') cannot be found on the ''pred exo'' sheet of the Excel data file. Please verify that this sheet is properly filled, and remember that dates are case-sensitive.'];
            msgbox(message);
            error('programme termination: data error');
        end
        % obtain the row location of the forecast end date
        [Fendlocation,~]=find(strcmp(strngs,Fenddate));
        % check that the end date for the forecast appears in the sheet; if not, return an error
        if isempty(Fendlocation)
            message=['Error: a forecast application is selected for a model that uses exogenous variables. Hence, predicted exogenous values should be supplied over the forecast periods. Yet the end date for forecasts (' Fenddate ') cannot be found on the ''pred exo'' sheet of the Excel data file. Please verify that this sheet is properly filled, and remember that dates are case-sensitive.'];
            msgbox(message);
            error('programme termination: data error');
        end

        % identify the strings for the exogenous variables
        % loop over exogenous
        for ii=1:numexo
            % try to find a column match for exogenous variable ii
            [~,location]=find(strcmp(strngs,exo{ii,1}));
            % if no match is found, return an error
            if isempty(location)
                message=['Error: a forecast application is selected for a model that uses exogenous variables. Hence, predicted exogenous values should be supplied over the forecast periods. Yet the exogenous variable ''' exo{ii,1} ''' cannot be found on the ''pred exo'' sheet of the Excel data file. Please verify that this sheet is properly filled, and remember that variable names are case-sensitive.'];
                msgbox(message);
                error('programme termination: data error');
                % else, record the value
            else
                pexolocation(ii,1)=location;
            end
        end

        % if everything was fine, reconstitute the matrix data_exo_p
        % initiate
        data_exo_p=[];
        % loop over exogenous variables
        for ii=1:numexo
            % initiate the predicted values for exogenous variable ii
            predexo=[];
            % loop over forecast periods
            for jj=1:Fperiods
                temp=strngs{Fstartlocation+jj-1,pexolocation(ii,1)};
                % if this entry is empty or NaN, return an error
                if (isempty(temp) || (temp<=inf)==0)
                    message=['Error: the predicted value for exogenous variable ' exo{ii,1} ' at forecast period ' strngs{Fstartlocation+jj,1} ' (and possibly other entries) is either empty or NaN. Please verify that the ''pred exo'' sheet of the Excel data file is properly filled.'];
                    msgbox(message);
                    error('programme termination: data error');
                    % if this entry is a number, record it
                else
                    predexo=[predexo;temp];
                end
            end
            % concatenate
            data_exo_p=[data_exo_p predexo];
        end

        % also, record the exogenous values on Excel
        % replace NaN entries by blanks
        strngs(cellfun(@(x) any(isnan(x)),strngs))={[]};
        % then save on Excel
        if pref.results==1
            warning off MATLAB:xlswrite:AddSheet
            xlswrite(fullfile(pref.results_path, [pref.results_sub '.xlsx']),strngs,'pred exo','A1');
            warning on MATLAB:xlswrite:AddSheet
        end
    end

end

author: european-central-bank

repo: BEAR-toolbox


text:
%% Test Newton Method 2
% Test Newton_n_sytem subroutine

% Clear memory
clear all; close all; clc;

% display format
format long

%% Examples
example = 3;
switch example
    case{1} % 2-by-2 Polynomials,
        % Define vector of functions 'f' and matrix of partial derivates 'df',
        f = inline('[x(1)^3+x(2)-1 ; x(2)^3-x(1)+1]');
        df = inline('[3*x(1)^2, 1  ; -1,  3*x(2)^2]');

        % Plot figures,
        g1 = 'x^3+y-1'; g2 = 'y^3-x+1'; xy_range = [-4,4,-4,4];
        figure(1); grid on;
        hold on;
        ez1 = ezplot(g1,xy_range);
        ez2 = ezplot(g2,xy_range);
        hold off;
        legend('g1(x,y)=x^3+y-1','g2(x,y)=y^3-x+1');
        set(ez1,'color',[0 0 1]);
        set(ez2,'color',[1 0 0]);
        title('Solving x^3+y-1 & y^3-x+1 for {x,y}')

        % Initial Guess,
        x0 = [.5;.5] ;

    case{2} % 2-by-2 Nonlinear system,
        % Define vector of functions 'f' and matrix of partial derivates 'df',
        f = inline('[10*x(1)^2+sin(x(2))-20 ; x(1)^4+5*x(2)-6]');
        df = inline('[20*x(1), cos(x(2))  ; 4*x(1)^3, 5]');

        % Plot figures,
        g1 = '10*x^2+sin(y)-20'; g2 = 'x^4+5*y-6'; xy_range = [-8,8,-8,8];
        figure(1); grid on;
        hold on;
        ez1 = ezplot(g1,xy_range);
        ez2 = ezplot(g2,xy_range);
        hold off;
        legend('g1(x,y)=10x^2+sin(y)-20','g2(x,y)=x^4+5y-6',1);
        set(ez1,'color',[0 0 1]);
        set(ez2,'color',[1 0 0]);
        title('Solving 10x^2+sin(y)=20 & x^4+5y=6 for {x,y}')

        % Initial Guess,
        x0 = [1;1];

        case{3} % 3-by-3 Polynomial System,
        % Define vector of functions 'f' and matrix of partial derivates 'df',
        f = inline('[x(1)^3-2*x(2)-2 ; x(1)^3-5*x(3)^2+7 ; x(2)*x(3)^2-1]');
        df = inline('[3*x(1)^2 -2 0;3*x(1)^2 0 -10*x(3);0 x(3)^2 2*x(2)*x(3)]');

        % Initial Guess,
        x0 = [1;1;1];
end

%% Convergence criteria
    % Define max number of iterations,
    n_max = 8;

    % Define tolerance,
    tol = 1e-6;
    tol = tol*ones(size(2,1));

%% Use Newton_n_sytem.m func,
tic;
x = x0; x_next = Newton_n_system(f,df,x); n = 1;
while and((n < n_max),(abs(x_next - x) > tol))
    x = x_next;
    x_next = Newton_n_system(f,df,x);
    n = n + 1;
end
fprintf('My solver: %1.12f \n',x) % show me the result
t1 = toc; % time for Costumb solver

%% Using Matlab solver,
tic;
x = x0; x = fsolve(f,x);
fprintf('Matlab solver: %1.12f \n',x) % show me the result
t2 = toc; % time for Matlab solver

% how fast is our Newton Solver,
fprintf('Costum Solver is: %1.1f times faster than Matlab solver \n',t2/t1)

repo: Aero-matlab


text:
% OP_KL_SHELLS: assemble the Kirchhoff-Love stiffness matrix K.
%
%   mat = op_KL_shells (spu, spv, msh, E_coeff, nu_coeff, t_coeff);
%
% INPUT:
%
%  spu:   structure representing the space of trial functions (see sp_vector/sp_evaluate_col)
%  spv:   structure representing the space of test functions  (see sp_vector/sp_evaluate_col)
%  msh:   structure containing the domain partition and the quadrature rule (see msh_cartesian/msh_evaluate_col)
%  E_coeff:  coefficients for the Young's modulus
%  nu_coeff: coefficients for the Poisson's ratio
%  t_coeff:  thickness of the shell
%
% OUTPUT:
%
%  mat:    assembled stiffness matrix
%
% Copyright (C) 2018, 2019 Pablo Antolin, Luca Coradello
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function mat = op_KL_shells (ref_sp_u, ref_sp_v, msh, E_coeff, nu_coeff, t_coeff)

   wq = msh.jacdet .* msh.quad_weights;

   [bending_stress, Kappa] = op_KL_bending_stress (ref_sp_u, ref_sp_v, msh, E_coeff .* wq, nu_coeff, t_coeff);
   [membrane_stress, Epsilon] = op_KL_membrane_stress (ref_sp_u, ref_sp_v, msh, E_coeff .* wq, nu_coeff, t_coeff);

   bending_stress = reshape(bending_stress, [], msh.nqn, ref_sp_u.nsh_max, msh.nel);
   membrane_stress = reshape(membrane_stress, [], msh.nqn, ref_sp_u.nsh_max, msh.nel);

   Kappa = reshape(Kappa, [], msh.nqn, ref_sp_v.nsh_max, msh.nel);
   Epsilon = reshape(Epsilon, [], msh.nqn, ref_sp_v.nsh_max, msh.nel);

   rows = zeros (msh.nel * ref_sp_u.nsh_max * ref_sp_v.nsh_max, 1);
   cols = zeros (msh.nel * ref_sp_u.nsh_max * ref_sp_v.nsh_max, 1);
   values = zeros (msh.nel * ref_sp_u.nsh_max * ref_sp_v.nsh_max, 1);

   ncounter = 0;

   for iel = 1:msh.nel
     if (all (msh.jacdet(:,iel)))

       ndof_u = ref_sp_u.nsh_max;
       ndof_v = ref_sp_v.nsh_max;

       bending_stress_iel = reshape(bending_stress (:, :, :, iel), [], msh.nqn, 1, ndof_u);
       membrane_stress_iel = reshape(membrane_stress (:, :, :, iel), [], msh.nqn, 1, ndof_u);

       Kappa_iel = reshape(Kappa (:, :, :, iel), [], msh.nqn, ndof_v, 1);
       Epsilon_iel = reshape(Epsilon (:, :, :, iel), [], msh.nqn, ndof_v, 1);

       tmp1 = sum (bsxfun (@times, bending_stress_iel, Kappa_iel), 1);
       tmp2 = sum (bsxfun (@times, membrane_stress_iel, Epsilon_iel), 1);

       elementary_values =   reshape (sum (tmp1, 2), ndof_v, ndof_u) + ...
                            reshape (sum (tmp2, 2), ndof_v, ndof_u);
       [rows_loc, cols_loc] = ndgrid (ref_sp_v.connectivity(:,iel), ref_sp_u.connectivity(:,iel));
       indices = rows_loc & cols_loc;
       rows(ncounter+(1:ref_sp_u.nsh(iel)*ref_sp_v.nsh(iel))) = rows_loc(indices);
       cols(ncounter+(1:ref_sp_u.nsh(iel)*ref_sp_v.nsh(iel))) = cols_loc(indices);
       values(ncounter+(1:ref_sp_u.nsh(iel)*ref_sp_v.nsh(iel))) = elementary_values(indices);
       ncounter = ncounter + ref_sp_u.nsh(iel)*ref_sp_v.nsh(iel);

%%%%%%%%%%%%%%%% old version, works when test and trial space are the same
%        % Here, the loops over the two sets of functions and the loop over the quadrature points are
%        % vectorized into a single operation.
%        ndof = ref_sp_u.nsh_max;
%        mat(ref_sp_v.connectivity(:, iel), ref_sp_u.connectivity(:, iel)) = ...
%            mat(ref_sp_v.connectivity(:, iel), ref_sp_u.connectivity(:, iel)) ...
%            + reshape(sum(sum(sum(repmat(bending_stress (:, :, :, :, iel), [1, 1, 1, 1, ndof]) ...
%                            .* permute(repmat(Kappa (:, :, :, :, iel), [1, 1, 1, 1, ndof]), [1, 2, 3, 5, 4])))) ...
%                  + sum(sum(sum(repmat(membrane_stress (:, :, :, :, iel), [1, 1, 1, 1, ndof]) ...
%                            .* permute(repmat(Epsilon (:, :, :, :, iel), [1, 1, 1, 1, ndof]), [1, 2, 3, 5, 4])))), [ndof, ndof]);

     else
       warning ('geopdes:jacdet_zero_at_quad_node',...
           'op_KL: singular map in element number %d', iel)
     end
   end
   mat = sparse (rows(1:ncounter), cols(1:ncounter), ...
                           values(1:ncounter), ref_sp_v.ndof, ref_sp_u.ndof);
end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [ x, w ] = jacobi_ek_compute ( n, alpha, beta )

%*****************************************************************************80
%
%% JACOBI_EK_COMPUTE: Elhay-Kautsky method for Gauss-Jacobi quadrature rule.
%
%  Discussion:
%
%    The integral:
%
%      integral ( -1 <= x <= 1 ) (1-x)^alpha * (1+x)^beta * f(x) dx
%
%    The quadrature rule:
%
%      sum ( 1 <= i <= n ) w(i) * f ( x(i) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2011
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Input, real ALPHA, BETA, the exponents of (1-X) and
%    (1+X) in the quadrature rule.  For simple Gauss-Legendre quadrature,
%    set ALPHA = BETA = 0.0.  -1.0 < ALPHA and -1.0 < BETA are required.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%

%
%  Define the zero-th moment.
%
  zemu = 2.0^( alpha + beta + 1.0 ) ...
    * gamma ( alpha + 1.0 ) ...
    * gamma ( beta + 1.0 ) ...
    / gamma ( 2.0 + alpha + beta );
%
%  Define the Jacobi matrix.
%
  x(1) = ( beta - alpha ) / ( 2.0 + alpha + beta );

  bj(1) = 4.0 * ( 1.0 + alpha ) * ( 1.0 + beta ) ...
    / ( ( 3.0 + alpha + beta ) * ( 2.0 + alpha + beta )^2 );

  for i = 2 : n
    abi = 2.0 * i + alpha + beta;
    x(i) = ( beta + alpha ) * ( beta - alpha ) / ( ( abi - 2.0 ) * abi );
    bj(i) = 4.0 * i * ( i + alpha ) * ( i + beta ) ...
      * ( i + alpha + beta ) ...
      / ( ( abi - 1.0 ) * ( abi + 1.0 ) * abi * abi );
  end

  bj(1:n) = sqrt ( bj(1:n) );

  w = zeros(n,1);

  w(1) = sqrt ( zemu );
  w(2:n) = 0.0;
%
%  Diagonalize the Jacobi matrix.
%
  [ x, w ] = imtqlx ( n, x, bj, w );

  w(1:n) = w(1:n).^2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [data, variable] = getRemoteVariable(variable,varargin)
%"getRemoveVariable"
%   Accesses, loads and returns data that is stored remotely.
%
%   Also returns an optional copy of the variable with remote data stored
%   in the data field, and the loaded field set to 1.
%
% JRA 10/14/04
% APA 11/23/05 - Added Local storage capability in tarball format
%
%Usage:
%   [data, variable] = getRemoteVariable(variable)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS planC
indexS = planC{end};

if isstruct(variable) & isfield(variable, 'remotePath') & isfield(variable, 'isLoaded') & isfield(variable, 'data') & isfield(variable, 'info') & strcmpi(variable.info, 'Remotely stored data.  Use getRemoteVariable to access.')'

    if variable.isLoaded
        data = variable.data;
    else
        switch upper(variable.storageType)
            case 'ZIP'
                load(fullfile(tempdir, variable.remotePath));
                data = CERRRemoteVariable;
                variable.data = data;
                variable.isLoaded = 1;
            case 'LOCAL'
                load(fullfile(variable.remotePath,variable.filename));
                data = CERRRemoteVariable;
                variable.data = data;
                variable.isLoaded = 1;
            case 'NETWORK'
                error('Local and Network storage not implemented yet.');
        end
    end
else
    error('Invalid remote data struct.')
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% l1qc_logbarrier.m
%
% Solve quadratically constrained l1 minimization:
% min ||x||_1   s.t.  ||Ax - b||_2 <= \epsilon
%
% Reformulate as the second-order cone program
% min_{x,u}  sum(u)   s.t.    x - u <= 0,
%                            -x - u <= 0,
%      1/2(||Ax-b||^2 - \epsilon^2) <= 0
% and use a log barrier algorithm.
%
% Usage:  xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)
%
% x0 - Nx1 vector, initial point.
%
% A - Either a handle to a function that takes a N vector and returns a K
%     vector , or a KxN matrix.  If A is a function handle, the algorithm
%     operates in "largescale" mode, solving the Newton systems via the
%     Conjugate Gradients algorithm.
%
% At - Handle to a function that takes a K vector and returns an N vector.
%      If A is a KxN matrix, At is ignored.
%
% b - Kx1 vector of observations.
%
% epsilon - scalar, constraint relaxation parameter
%
% lbtol - The log barrier algorithm terminates when the duality gap <= lbtol.
%         Also, the number of log barrier iterations is completely
%         determined by lbtol.
%         Default = 1e-3.
%
% mu - Factor by which to increase the barrier constant at each iteration.
%      Default = 10.
%
% cgtol - Tolerance for Conjugate Gradients; ignored if A is a matrix.
%     Default = 1e-8.
%
% cgmaxiter - Maximum number of iterations for Conjugate Gradients; ignored
%     if A is a matrix.
%     Default = 200.
%
% Written by: Justin Romberg, Caltech
% Email: jrom@acm.caltech.edu
% Created: October 2005
%

function xp = l1qc_logbarrier(x0, A, At, b, epsilon, lbtol, mu, cgtol, cgmaxiter)

largescale = isa(A,'function_handle');

if (nargin < 6), lbtol = 1e-3; end
if (nargin < 7), mu = 10; end
if (nargin < 8), cgtol = 1e-8; end
if (nargin < 9), cgmaxiter = 200; end

newtontol = lbtol;
newtonmaxiter = 50;

N = length(x0);

% starting point --- make sure that it is feasible
if (largescale)
  if (norm(A(x0)-b) > epsilon)
    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
    AAt = @(z) A(At(z));
    [w, cgres] = cgsolve(AAt, b, cgtol, cgmaxiter, 0);
    if (cgres > 1/2)
      disp('A*At is ill-conditioned: cannot find starting point');
      xp = x0;
      return;
    end
    x0 = At(w);
  end
else
  if (norm(A*x0-b) > epsilon)
    disp('Starting point infeasible; using x0 = At*inv(AAt)*y.');
    opts.POSDEF = true; opts.SYM = true;
    [w, hcond] = linsolve(A*A', b, opts);
    if (hcond < 1e-14)
      disp('A*At is ill-conditioned: cannot find starting point');
      xp = x0;
      return;
    end
    x0 = A'*w;
  end
end
x = x0;
u = (0.95)*abs(x0) + (0.10)*max(abs(x0));

% disp(sprintf('Original l1 norm = %.3f, original functional = %.3f', sum(abs(x0)), sum(u)));

% choose initial value of tau so that the duality gap after the first
% step will be about the origial norm
tau = max((2*N+1)/sum(abs(x0)), 1);

lbiter = ceil((log(2*N+1)-log(lbtol)-log(tau))/log(mu));
% disp(sprintf('Number of log barrier iterations = %d\n', lbiter));

totaliter = 0;
dispProgress('Log barrier', 0, lbiter);
for ii = 1:lbiter

  [xp, up, ntiter] = l1qc_newton(x, u, A, At, b, epsilon, tau, newtontol, newtonmaxiter, cgtol, cgmaxiter);
  totaliter = totaliter + ntiter;

%   disp(sprintf('\nLog barrier iter = %d, l1 = %.3f, functional = %8.3f, tau = %8.3e, total newton iter = %d\n', ...
%     ii, sum(abs(xp)), sum(up), tau, totaliter));
  dispProgress('Log barrier',ii/lbiter);

  x = xp;
  u = up;

  tau = mu*tau;

end
dispProgress('Log barrier', 'Close');

end

author: thomaskuestner


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% getalp.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [alp,lba,uba,ier]=getalp(alpu,alpo,gTp,pTGp);
% get minimizer alp in [alpu,alpo] for a univariate quadratic
%	q(alp)=alp*gTp+0.5*alp^2*pTGp
% lba	lower bound active
% uba	upper bound active
%
% ier	 0 (finite minimizer)
%	 1 (unbounded minimum)
%
function [alp,lba,uba,ier]=getalp(alpu,alpo,gTp,pTGp);

lba=0;
uba=0;

% determine unboundedness
ier=0;
if alpu==-inf & ( pTGp<0 | (pTGp==0 & gTp>0) ),
  ier=1; lba=1;
end;
if alpo==inf & (pTGp<0 | (pTGp==0 & gTp<0) ),
  ier=1; uba=1;
end;
if ier, alp=NaN; return; end;

% determine activity
if pTGp==0 & gTp==0,
  alp=0;
elseif pTGp<=0,
  % concave case minimal at a bound
  if alpu==-inf,     lba=0;
  elseif alpo== inf, lba=1;
  else               lba = (2*gTp+(alpu+alpo)*pTGp>0);
  end;
  uba = ~lba;
else
  alp=-gTp/pTGp;          % unconstrained optimal step
  lba = (alp <= alpu);    % lower bound active
  uba = (alp >= alpo);    % upper bound active
end;

if lba, alp=alpu; end;
if uba, alp=alpo; end;

% print?
if abs(alp)==inf, gTp,pTGp,alpu,alpo,alp,lba,uba,ier, end;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function doomsday_test02 ( )

%*****************************************************************************80
%
%% DOOMSDAY_TEST02 tests DOOMSDAY against a number of known values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 May 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'DOOMSDAY_TEST02\n' );
  fprintf ( 1, '  WEEKDAY_VALUES supplies a list of dates and weekdays.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  YYYY  MM  DD  Weekday    Weekday\n' );
  fprintf ( 1, '                Tabulated  Computed\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, y, m, d, w1 ] = weekday_values ( n_data );

    if ( n_data <= 0 )
      break;
    end
%
%  The transition from Julian to Gregorian calendars occurred in 1582
%  (for some people).  The data in "WEEKDAY_VALUES" before the transition
%  is stored in Julian format, which DOOMSDAY_GREGORIAN can't handle.
%  So let's just refuse to handle 1582 or earlier%
%
    if ( y <= 1582 )
      continue;
    end

    w2 = doomsday_gregorian ( y, m, d );

    s1 = weekday_to_name_common ( w1 );
    s2 = weekday_to_name_common ( w2 );

    fprintf ( 1, '  %4d  %2d  %2d  %10s  %10s\n', y, m, d, s1, s2 );

  end

  return
end


text:
%% This file will be used to plot the different result of DL-SINDy
% Coded By: K
% Last Updated: 2019/06/23
%%
clc;clear all; close all;

%% Add path
addpath('.\PDE_Find_Results')

%% Load data
Files=dir('.\PDE_Find_Results\*.mat');

Error=zeros(length(Files),1);
for k=1:length(Files)
   FileNames=Files(k).name;
   load(FileNames)
   Error(k,1)=minVal_Final;
   digits(4)
   Expression
end

%% Plot
close all
g0=[0;0.0001;0.001;0.01;0.1;1];
figure(1)
loglog(g0,Error,'marker','o','linewidth',2.5)
set(gca,'FontSize',18);
set(gcf,'Position',[100 100 600 400]);
set(gcf,'PaperPositionMode','auto');
grid on

author: dynamicslab

repo: SINDy-PI


text:
%% Parameter setting %%

param.nx = 64;
param.ny = 64;
param.nz = 64;

param.sx = 64; % mm
param.sy = 64; % mm
param.sz = 64; % mm

% Detector setting, according to Varian Trilogy OBI (real size)
param.su = 256;	% mm
param.sv = 200;	% mm

%The real detector panel pixel density (number of pixels)
param.nu = 128;
param.nv = 100;

% X-ray source and detector setting
param.DSD = 900;    %  Distance source to detector
param.DSO = 400;	%  X-ray source to object axis distance

% angle setting
dir = -1;   % gantry rotating direction
param.deg = 0:4.3:360;
param.deg = param.deg*dir;
param.nProj = length(param.deg);

% filter='ram-lak','cosine', 'hamming', 'hann'
param.filter='hamming'; % high pass for sintetic images

param.dx = param.sx/param.nx;
param.dy = param.sy/param.ny;
param.dz = param.sz/param.nz;
param.du = param.su/param.nu;
param.dv = param.sv/param.nv;

param.off_u = 0; param.off_v = 0; % detector rotation shift (real size)

% % % For fast CPU calculation % % %
param.xs = [-(param.nx-1)/2:1:(param.nx-1)/2]*param.dx;
param.ys = [-(param.ny-1)/2:1:(param.ny-1)/2]*param.dy;
param.zs = [-(param.nz-1)/2:1:(param.nz-1)/2]*param.dz;

param.us = (-(param.nu-1)/2:1:(param.nu-1)/2)*param.du + param.off_u;
param.vs = (-(param.nv-1)/2:1:(param.nv-1)/2)*param.dv + param.off_v;

%% projection
load Phantom64.mat % img

proj = CTprojection(img,param);

for i=1:param.nProj
    figure(1); imagesc(max(proj(:,:,i)',0)); axis off; axis equal; colormap gray; colorbar;
    title(num2str(i));
    pause(0.01);
end

%% filtered backprojection
% filter='ram-lak','shepp-logan','cosine', 'hamming', 'hann'
param.filter='hamming';
% param.filter='ram-lak';

Reconimg = CTbackprojection(proj, param);

for i=1:param.nz
    figure(2); imagesc(max(Reconimg(:,:,i),0)); axis off; axis equal; colormap gray; colorbar;
    title(num2str(i));
    pause(0.01);
end

%% MLEM
param.filter = 'none';

img = ones(param.nx, param.ny, param.nz,'single');
Norimg = CTbackprojection(ones(param.nu, param.nv, param.nProj, 'single'), param);

for iter = 1:50

    proj_ratio = proj./CTprojection(img,param);
    proj_ratio(isnan(proj_ratio)) = 0;
    proj_ratio(isinf(proj_ratio)) = 0;

    img_ratio = CTbackprojection(proj_ratio, param)./Norimg;
    img_ratio(isnan(img_ratio)) = 0;
    img_ratio(isinf(img_ratio)) = 0;

    img = img.*img_ratio;

    figure(3); imagesc(max(img(:,:,round(end/2)),0)); axis off; axis equal; colormap gray; colorbar;
    title(['Iteration - ',num2str(iter)]);
    pause(0.1);
end

author: Sable

repo: mcbench-benchmarks


text:
function doseStat = dispDoseStats(doseBinsV, volHistV, name, nameVol, planC, indexS, opt)
%Command line display of basic dose statistics
%doseBinsV is a vector of the midpoint doseBin values.
%volHistV is either a histogram of volumes or surface areas.
%LM: 14 Oct 02, JOD.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

switch lower(opt)

  case 'standarddose'
    disp('-----------------------')
    disp('')
    disp(['Structure is:  ' name])

    totalVol = sum(volHistV);
    disp(['Total volume is:  ' num2str(totalVol) ' cubic cm.'])

    disp(['Dose map name is:  ' nameVol])

    meanD = sum(doseBinsV(:).*volHistV(:))/sum(volHistV);
    disp(['Mean dose is:  ' num2str(meanD)])

    ind = max(find([volHistV~=0]));
    maxD = doseBinsV(ind);
    disp(['Max dose is:  ' num2str(maxD)])

    ind = min(find([volHistV~=0]));
    minD = doseBinsV(ind);

    disp(['Min dose is:  ' num2str(minD)])
    disp('')
    disp('-----------------------')

    doseStat.volume = totalVol;

  case 'dshdose'

    disp('-----------------------')
    disp('')
    disp(['Structure is:  ' name])

    areaV = volHistV;  %actually areas, not volumes.
    dosesV = doseBinsV;

    totalArea = sum(areaV);
    disp(['Total surface area is:  ' num2str(totalArea) ' square cm.'])

    disp(['Dose map name is:  ' nameVol])

    meanD = sum(dosesV(:).*areaV(:))/sum(areaV);
    disp(['Mean surface dose is:  ' num2str(meanD)])

    maxD = max(dosesV);
    disp(['Max dose is:  ' num2str(maxD)])

    minD = min(dosesV);

    disp(['Min dose is:  ' num2str(minD)])
    disp('')
    disp('-----------------------')

    doseStat.volume = totalArea;

end

doseStat.min    = minD;
doseStat.mean   = meanD;
doseStat.max    = maxD;


text:
function s = ymdf_to_s_hebrew ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_S_HEBREW "prints" a Hebrew YMDF date into a string.
%
%  Format:
%
%    DayNumber.Fraction MonthName YearNumber AM
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, string S, a representation of the date.
%
  sm = month_to_month_name_hebrew ( y, m );
  s = sprintf ( '%05.2f %s %d AM', d + f, sm, y );

  return
end


text:
% Application of Oleinik's entropy condition

% Copyright 2001 P. Wesseling
% This program and its subprograms may be freely used, modified and distributed
% under the GNU General Public License: http://www.gnu.org/copyleft/gpl.html

% Theory in Section 9.4 of:

% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% This program generates Fig. 9.11 in the book

		% ...............Input.........................................
a = 0.5;	% Parameter in flux function on page 371
		% ............End of input.....................................

x = 0.0:0.05:1.0; den = x.*x + 0.5*(1-x).^2;
f = (x.*x)./den;	% Buckley-Leverett flux function

x1 = 1/sqrt(3);
f1 = (x1*x1)/(x1*x1 + 0.5*(1-x1)^2);	%Point where tangent meets curve

figure(1), clf, hold on, plot(x,f,'-')
xx = [0  1];  xx = x1*xx; yy = [0  1]; yy = f1*yy; plot(xx,yy,'-')
zz = [0  f1]; xx = [x1 x1];			   plot(xx,zz,'--')


text:
% sdae_get_hidden
% Copyright (C) 2011 KyungHyun Cho, Tapani Raiko, Alexander Ilin
%
%This program is free software; you can redistribute it and/or
%modify it under the terms of the GNU General Public License
%as published by the Free Software Foundation; either version 2
%of the License, or (at your option) any later version.
%
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
%
function [h_mf] = sdae_get_hidden(x0, S, target_sparsity)

if nargin < 3
    target_sparsity = 0;
end

layers = S.structure.layers;
n_layers = length(layers);

h_mf = x0;

for l = 2:n_layers
    h_mf = bsxfun(@plus, h_mf * S.W{l-1}, S.biases{l}');

    if l < n_layers || S.bottleneck.binary
        h_mf = sigmoid(h_mf, S.hidden.use_tanh);
    end
end

if S.bottleneck.binary
    if target_sparsity > 0
        avg_acts = mean(h_mf, 1);
        diff_acts = max(avg_acts - (1 - target_sparsity), 0);
        h_mf = min(max(bsxfun(@minus, h_mf, diff_acts), 0), 1);
    end
end

author: kyunghyuncho

repo: deepmat


text:
function linplus_test535 ( )

%*****************************************************************************80
%
%% TEST535 tests R8S3_READ, R8S3_READ_SIZE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  input_file = 'r8s3_matrix.txt';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST535\n' );
  fprintf ( 1, '  For a R8S3 system,\n' );
  fprintf ( 1, '  R8S3_READ_SIZE reads the size of the matrix.\n' );
  fprintf ( 1, '  R8S3_READ reads the matrix.\n' );

  [ n, nz_num ] = r8s3_read_size ( input_file );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R8S3_READ_SIZE reports matrix size data:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N =         %d\n', n );
  fprintf ( 1, '  Matrix nonzeros NZ_NUM = %d\n', nz_num );

  [ row, col, a ] = r8s3_read ( input_file, n, nz_num );

  isym = 0;

  r8s3_print_some ( n, n, nz_num, isym, row, col, a, 1, 1, ...
    10, 10, '  Initial 10x10 block of recovered R8S3 matrix:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Deleting the matrix data file "%s"\n', input_file );

  return
end


text:
function [A,z]=jacobiancsd(fun,x)
% JACOBIANCSD    Complex Step Jacobian
% J = jacobiancsd(f,x) returns the numberical (m x n) Jacobian matrix of a
% m-vector function, f(x) at the refdrence point, x (n-vector).
% [J,z] = jacobiancsd(f,x) also returns the function value, z=f(x).
%
% Example
% f=@(x)[x(1)*x(2)^2;x(1)*x(3)^2;x(1)^2];
% x=1:3;
% [J,z]=jacobiancsd(f,x)
% J =
%     4     4     0
%     9     0     6
%     2     0     0
% z =
%     4
%     9
%     1
%
% By Yi Cao at Cranfield University, 02/01/2008
%
z=fun(x);                       % function value
n=numel(x);                     % size of independent
m=numel(z);                     % size of dependent
A=zeros(m,n);                   % allocate memory for the Jacobian matrix
h=n*eps;                        % differentiation step size
for k=1:n                       % loop for each independent variable
    x1=x;                       % reference point
    x1(k)=x1(k)+h*i;            % increment in kth independent variable
    A(:,k)=imag(fun(x1))/h;     % complex step differentiation
end


text:
function [Gu,Gs,Gn,f] = spm_csd_mtf_gu(P,M)
% Spectral desnities of innovations and noise for DCM for CSD
% FORMAT [Gu,Gs,Gn,f] = spm_csd_mtf_gu(P,M)
% FORMAT [Gu,Gs,Gn,f] = spm_csd_mtf_gu(P,f)
%
% P   - parameters
% M   - neural mass model structure (with M.Hz)
% f   - frequencies of interest (Hz)
%
% Gu  - neuronal innovations
% Gn  - channel noise (non-specific)
% Gs  - channel noise (specific)
%
% f   - frequency
%
% fluctuations and noise parameters: for n regions and c channels
%--------------------------------------------------------------------------
%    pE.a(2,n) - neuronal fluctuations        - amplitude and exponent
%    pE.b(2,c) - channel noise (non-specific) - amplitude and exponent
%    pE.c(2,c) - channel noise (specific)     - amplitude and exponent
%    pE.d(8,n) - neuronal fluctuations        - basis set coefficients
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_csd_mtf_gu.m 6856 2016-08-10 17:55:05Z karl $

% frequencies of interest
%--------------------------------------------------------------------------
try, f = M.Hz(:); catch, f = M(:); end

% Number of sources and fequencies
%--------------------------------------------------------------------------
if isfield(P,'d')
    ns = max(size(P.a,2),size(P.d,2));
else
    ns = size(P.a,2);
end
nf     = size(f,1);

% spectrum of neuronal innovations (Gu)
%==========================================================================
for i = 1:size(P.a,2)
    Gu(:,i) = exp(P.a(1,i))*f.^(-exp(P.a(2,i)));
end

if size(Gu,2) == 1, Gu = Gu*ones(1,ns); end

% add structured innovations - a discrete cosine set of order length(P.d)
%--------------------------------------------------------------------------
if isfield(P,'d')
    nd = size(P.d,1);
    X  = spm_dctmtx(nf,nd + 1);
    Mu = exp(X(:,2:end)*P.d);
else
    Mu = ones(nf,1);
end

if size(Mu,2) == 1, Mu = Mu*ones(1,ns); end

% return unless channel noise is required
%--------------------------------------------------------------------------
Gu     = Gu.*Mu;
if nargout < 2, return, end

% spectrum of channel noise (non-specific)
%==========================================================================
Gn  = exp(P.b(1) - 2)*f.^(-exp(P.b(2)));

% and spectrum of channel noise (specific: with the same exponent)
%--------------------------------------------------------------------------
for i = 1:size(P.c,2)
    Gs(:,i) = exp(P.c(1,i) - 2)*f.^(-exp(P.c(2,1)));
end

author: spm


text:
%MAIN_fsm

clc; clear;

%%%% Derive equations of motion %%%%
% EoM_flight();
% EoM_slide();
% EoM_hinge();

%%%% Interesting Parameters %%%%
MoI = 0.333;  %Dimensionless moment of inertia
    % 0   => Point mass at center of rod
    % 1/3 => Slender rod (uniform distribution)
    % 1   => Point mass at both ends of the rod
CoeffFriction = 0.2;
%%%% ~~~~~~~~~~~~~~~~~~~~~~ %%%%

P.g = 9.81;
P.L = 1;  %Distance between END OF STICK AND CENTER OF MASS
P.FullLength = 2*P.L;   %Length of the entire stick
P.m = 1;
P.u = CoeffFriction;   %Coefficient of friction
P.I = MoI*P.m*P.L^2;  %Theoretical Maximum - two point masses
setup.P = P;

setup.IC.th = 1e-2;
setup.IC.x = 0;
setup.IC.y = 0;
setup.IC.dth = topple_angularRate(setup.IC.th,P);
setup.IC.dx = 0;
setup.IC.dy = 0;

setup.Tspan = [0,10];  %only used for timeout of the simulation
setup.tol = 1e-12;   %Accuracy of the intergation method
setup.dataFreq = 750;   %How much data to return?
setup.solver = @ode45; %Integration method
setup.odeMaxStep = 1e-1;  %Biggest allowable step

D = runSimulation(setup);

plotData(D);

for i=1:length(D.phase)
    disp([D.phase{i} '  ->  ' D.code{i}])
end

timeRate = 0.5;   % slow motion < 1 < fast forward
animation(D,timeRate);

author: MatthewPeterKelly

repo: dscTutorials