text:
function  [c1, lmin1, lmax1, as, pm, sig8] = cmbaccur(van, omt, bmark, one_percent, lmax, stfname)
% cmbaccur.m :
% this program calculates the radiation and matter spectra to an accuracy of 1% or 0.1% for the open, flat and
% closed Friedman-Robertson-Walker universe.
% The radiation spectra include the temperature- polarisation and cross-polarisation spectra.
% We follow the method of Uros Seljak amd Matias Zaldarriaga in 'A line-of-sight integration approach
% to Cosmic microwave background anisotropies', ApJ, nr 469, 437-444 (1996), see also :
% Ma and Bertschinger, ApJ nr 455, 7-25 (1995) and Zaldarriaga, Seljak, Bertschinger
% ApJ nr 494, 491 (1998) and Zaldarriaga and Seljak, ApJ nr129, 431 (2000) and W Hu e.a, 'A complete treatment
% of CMB anisotropies in a FRW universe', Phys Rev D57, 3290-3301 (1998), for reionisation see N Sugiyama,
% ApJ nr 419 L1-L4 (1993) and ApjS nr 100, 281 (1995)
%
% input : the "vanilla-set" of Cosmological parameters, omt, the curvature of the universe,
%   bmark : a benchmark (file) to compare with (see remarks), one_percent: the accuracy
%   is 1% resp 0.1% for ~one_percent==0 resp one_percent==0, the default is 1%, lmax, the maximum
%   multipoint number, the default lmax value is lmax=1500 (lmax must be <= 1500. otherwise a special
%   la-vector has to be defined), name of the startfile for the ultra-spherical function.
% no inputs : if no inputs are given the functions takes the default values (see remarks),
% if 5th argument (lmax) is missing, the default value lmax=1500 is taken
%
% output :the normalisation constant (asc), the positions of the peaks (minima and maxima)
% in the temperature anisotropy spectrum (lmin & lmax), as: the temperature anisotropy spectrum, as.ctt
% the E-polarisation spectrum, as.cee and the cross-correlation spectrum, as.cte, pm: the matter spectrum

% D Vangheluwe 29 jun 2005 finished 4 july 2005
% remark 1: the vanilla-set is (As, omega_lambda, omega_d, omega_b, ns, tau), which are : the amplitude correction
% at l=10, dark energy density, dark matter density, baryon density, scalar spectral index, optical depth.
% The vanilla-set is converted internally to the parameters (c1, om, omt, omb, h, nprimtilt, zri), which are:
% the amplitude correction at l=10, dark matter omega, total omega, baryon omega, Hubble factor, scalar spectral index,
% reionisation redshift, see below for the default values. For example: omt = omega_lambda + omega_b + omega_d
% remark 1 this is a compilation of the programs cmbacc1.m and ambacc2.m for resp open (K<=0) and closed space (K>0)
% remark 2 the discrete values of la for closed space are chosen different from the ones for open space : reason that
%  for closed space the peaks in the spectrum shift to lower la and we need more points there : therefore the
%  step change in la is chosen at la =100 resp la = 250 for open/flat resp closed space.
% remark 3 :the lamda content of the universe is automatically : omt - omr -om, for example with omt = 1.1,
%  om = 0.27, h = 0.72 and omr = omp/(1 - fv)= 8.01e-5, we have oml = 0.83, we have chosen here not to
%  specify the lamda energy content of primordial space. For the Vanilla parameters see ambaccur(.)
% remark 4: the different path taken for closed resp open/flat space can be seen at the if (K>0) statements.
%  for K>0 we take a different set of k1 values to solve the Bolzamn equations and a different set of k2
%  (also 'beta' in the literature) values for the calculation of the anistropy spectrum
% remark 6: the connections with the input of cmbfast and our input is a bit difficult : we give the connection
%  for resp omega of the vacuum (omegav), cold dark matter (omegac) and baryons (omegab) :
%  omegav = omt - om, omegac = om - omb, omegab = omb (om is the omega for all matter, baryons+dark matter and omt
%  is the total omega)
% remark 7: try the vanilla-lite model from Tegmark et al. astro-ph/0310723, 15jan 2004:
%   van = [0.89, 0.72, 0.12, 0.024, 1, 0.17] and omt =1, it is the best fit to WMAP2003 +SDSS assuming ns==1 and tensor field zero.
%   or the old cdm model, a universe without dark energy : van = [1, 0.954, 0.0, 0.0115, 1, 0] and omt =1.
% remark 5 (4 aug 2005): if the spectrum needs only to be calculated for a multipole range smaller than the
%  default lmax0=1500, specify this value <1500 as the last argument, the accuracy is
%  kept in line with this choice (several changes are necessary such as : smaller nkst, extended k1max etc.)
%  The calculation time is proportional to sqrt(lmax), not with lmax!

global GL_dt_table   GL_pp_dtau;
global  GL_cmb_c  GL_cmb_h0   GL_cmb_t0   GL_cmb_T0   GL_cmb_rv   GL_cmb_fv   GL_cmb_kg1   GL_cmb_ka1   ...
  GL_cmb_yp   GL_cmb_ncr   GL_cmb_cr   GL_cmb_pcm   GL_cmb_dha;
% give the global physical constants a value
cmbglobl;

% tolerance
tol = 1e-6;
trace = [];

%cosmological parameters : (As, omega_lambda, omega_d, omega_b, ns, zri) and omt for spacial curvature
if nargin == 0  %take the default vanilla-set (which gives the best fit to wmap2003)
% the default set is : c10=1, om=0.27, omb=0.046, h=0.72, ns=0.99, zri= 17.28
    van = [1, 0.73, 0.1161216, 0.0238464, 0.99, 0.17295];
    omt = 1;
    lmax = 1500;
% the default accuracy is : one_percent == 1, for an of accuracy 1%, the accuracy = 0.1% only for one_percent == 0.
    one_percent = 1
% benchmark calculations to compare : available are cmbfast for h=0.72: 'h72', 'h72i' for reionisation with zri = 17.3
%  'h72o1' for h = 0.72 open universe with omt=1.1, 'h72c1' for h = 0.72 closed universe with omt=0.9,
% 'cdm', for the CDM model (flat space) (see for more info cmbdata3.m)
    bmark = 'h72i';
else
   if nargin == 1, error('the benchmark must be defined'), return; end
   if nargin == 2, error('the accuracy must be defined'), return; end
   if nargin == 4, lmax = 1500; end
end
if  ~(one_percent == 0), one_percent = 1; end;

% the hubble factor h
%h = 0.71
%h = 0.72
if (omt - van(2) > 0)
   h = sqrt((van(3) + van(4))/(omt - van(2)))
else
   error('omt < van(2) not allowed');
   return
end

% baryon content of the universe
%omb = 0.0444
%omb = 0.046
omb = van(4)/h^2

% matter content of the universe
%om = 0.27
%om = 0.27
om = (van(3) + van(4))/h^2

% redshift for reionisation
%zri = 0
%zri = 17.28   % tauri = 0.170
%zri = 0
zri = 92*(0.03 * h * van(6)/van(4))^(2/3) * om^(1/3)

% primordial tilt : the exponent of the power spectrum of matter
%nprimtilt = 0.99;
%nprimtilt = 0.99;
nprimtilt = van(5)

% normalisation of the anisotropy spectrum
%c10 = 1.00;
c10 = van(1);

% the cmb temperature
%cmbtemp = 2.725;
cmbtemp = GL_cmb_t0;
% velocity of light in m/s
%c = 2.998e8;
c = GL_cmb_c;
% the hubble constant at present in h Mpc^-1, see my notes p71
%h0 = 1e5/c;
h0 = GL_cmb_h0;
% the ratio of neutrino to total radiation density fv, assuming three massless neutrinos
%rv = (21/8) * (4/11)^(4/3);
%fv = rv/(1 + rv);
fv = GL_cmb_fv;
% amplitude of the primordial density perturbation at the horizon, see Dodelson (8.76)
%dha = 4.6e-5;
dha = GL_cmb_dha;
dha2 = dha^2;
% the ratio of the radiation density and the critical density, see Dodelson (omega_gamma)
%kg1 = 2.47e-5;
kg1 = GL_cmb_kg1;
%ka1 = kg1/(1 - fv);
ka1 = GL_cmb_ka1;

aeq = ka1/(om * h^2);
keq = h0 * h * sqrt(om*2/aeq);
% omega of radiation and neutrinos separately
omp = kg1/h^2;
omr = omp/(1 - fv);
omn = omp * fv/(1 - fv);
% the constant K for curvature energy content from the Friedman equation
K = (omt - 1) * (h0 * h)^2;

if  nargin < 6
  if K > 0, stfname = 'stfile1';
  else, stfname = 'stfile';
  end
end

tol_ode23 = 1e-5;
if one_percent == 1
  tol_ode23 = 2e-4;  % default 2e-4
end

% find the recombination curve and the inverse optical depth, dtau in Mpc^-1,
% tau and the visibility function (vsb1) :
dt = recomb(om, omb, omt, h, zri);
%sdata = load('cmbxe');
%dt = sdata.dt;
GL_dt_table = dt.dt_table;
%tauri = dt.tauri;
tauri = van(6);

zdec = dt.zdec;
sizetable = size(GL_dt_table, 2)
at1 = GL_dt_table(1,:);
dtau1 = GL_dt_table(3,:);
%tau1 = GL_dt_table(4,:);
%vsb1 = GL_dt_table(5,:);
% make a spline of dtau
GL_pp_dtau = spline(GL_dt_table(1,:), GL_dt_table(3,:));

% calculate the current conformal time :ctc0 and ctime at reionisation and at decoupling
ctc0 = quadl(@conftim1, 0, 1, tol, trace, om, omr, omt, h);
ctad = quadl(@conftim1, 0, 1/(zdec + 1), tol, trace, om, omr, omt, h);
ctri = quadl(@conftim1, 0, 1/(zri + 1), tol, trace, om, omr, omt, h);

% if the spectrum has become isotropic for K > 0 (all photons are coming from one point) break off:
if (K > 0)
    ct4sym = 0.5 * pi/sqrt(K);
    if  (2*ct4sym < ctc0),  message('the spectrum is isotropic, return'); return; end
    if sqrt(K) * ctc0 > 2, message('curvature K * ctc0 > 2'); end
end

% define the range of l (angle number)
%lmax = 1500;
lmin = 2;

% the wavenumber kw in Mpc^1 (range : 0.2 < kw < 2e-4)
% define the k-range
if one_percent == 1
%  nkst = 44
  nkst = 30 + 10*(lmax - 900)/600;   % adapt the k1step for lmax~1500
else
  nkst = 88
end

% find the maximum multipole number in the startfile for K>0 (lmax0 =default 1500)
stv1 = load(stfname);
lmax0 = stv1.ust.la(end);
if (K > 0)
    k1min = 3*sqrt(K);
%    k1int = 4e-3;
    k1int = 8e-3;
% if K*ctc0 > 2, it is necessary to increase k1max : important change for K>0!!
% and extend the range of k1 in case lmax < maximum of the startfile (lmax0 is default 1500)
    k1max = max(2*sqrt(lmax*lmax0)/ctc0, 1.5*sqrt(lmax*lmax0)*sqrt(K));
    k1step = (k1int - k1min)/8;
    k1step2 = (k1max - k1int)/nkst;
    k1 = [[k1min : k1step : (k1int - k1step)], [k1int : k1step2 : k1max]];
else
    k1min = 1e-5;
    k1int = 8e-3;
% extend the range of k1 in case lmax < maximum of the startfile (lmax0 is default 1500)
    k1max = 2*sqrt(lmax*lmax0)/ctc0;
    k1step = 0.0012;
    k1step2 = (k1max - k1int)/nkst;
    k1 = [k1min, 4e-5, 8e-5, 2e-4, 4e-4, 6e-4, 8e-4 : k1step : k1int, k1int+k1step2 : k1step2 : k1max];
end  %define k1 range
lk1 = size(k1, 2)

% take smaller steps in the same range of the conformal time now ct2:
step0 = ctc0/7000;
if one_percent == 1
  ct2step = step0;
else
  ct2step = step0/2;
end
ct2min = 0.5;
ct2max = ctc0;
ct2 = ct2min : ct2step : ct2max;
lct2 = size(ct2, 2);

%555555555555555555
% solve the Bolzman equations :
[table_srt, table_srp, deltam, a2] = ...
        bolzmslv(ct2, k1, tol_ode23, om, omb, omt, h, ctc0, ctad, zri, ctri, K, one_percent);

% find the conformal time at recombination :
vsb2 = spline(GL_dt_table(1,:), GL_dt_table(5,:), a2);
[vsbm irm] = max(vsb2);
etar = ct2(irm);
% calculate the thickness of the recombination surface in conformal time
ids = find(vsb2 >= 0.5* vsbm);
% half width ct2 points are:
hw1 = ct2(ids(1));
hw2 = ct2(ids(end));
trs = hw2 - hw1;
% number of points over the half width of the visibility function :
nphw = length(ids)
% find an accurate value for the decoupling time (recombination)
ct3 = hw1:0.01:hw2;
vsb3 = spline(ct2, vsb2, ct3);
[vsbm irm] = max(vsb3);
ctad1 = ct3(irm);

% define a k2 range:
if (K > 0)
% define a k-range with more points ->k2  :where kb = 1/3..
    logstep = 0;
    nn2min = fix(k1min/sqrt(K));
    nn2max = fix(k1max/sqrt(K));
    nksteps = 1.5* sqrt(lmax * lmax0);
    nn2step = max(1, round((nn2max - nn2min)/nksteps));
    nn2step = min(2, nn2step);   % limit to nn2step==2
    nn2 = nn2min : nn2step : nn2max;
    k2 = nn2 * sqrt(K);
else

    logstep = 1;
    k2min = k1min;
    k2max = k1max;
    if one_percent == 1
       ek2step = 0.002;
    else
       ek2step = 0.0015;
    end
    ek2 = log(k2min)/log(10) : ek2step : log(k2max)/log(10);
    k2 = 10 .^ ek2;
end % define k2 range
lk2 = size(k2, 2);
ik2 = floor(lk2/2);
ik4 = floor(lk2/4);

% calculate the anisotropy- and polariation spectrum for l > 2 in 2 steps: step1 :l<100 and step2 :l>100
% for l > 100 :select the interval in the conformal time range ct2 where the source is none zero
% for l < 100 :the source extends over a rather large cf time range of 200-ctc0 and we found that
% the ctime integration range must be extended in order to achieve 0.1% accuracy for l < 100
% With reionisation (zri > 0) for 1% accuracy, integration should be extended over the total ct time range.

%##anisotropy spectrum, step1 : l < 100  (2 and lct2-1 to stay within interpolation limits)
ict4min = 2;
%ict4max = lct2 - 1;
ict4max = lct2;
if one_percent == 1
   ict4step = 2;
else
   ict4step = 1;
end
ict4 = ict4min : ict4step : ict4max;
ct4 = ct2(ict4);
lct4 = size(ict4, 2);

% take as much splines ready for k2 interpolations as there are cftime points in ct4 :
% for each time we have one spline (the splines are in a column vector)
% this steps needs a lot of memory : if short of memory use srcpint.m with pp4, pp5 as argument : this
% is the case for one_percent ~= 1
pp4 = spline(k1, table_srt(:, ict4)');
pp5 = spline(k1, table_srp(:, ict4)');
if one_percent == 1
   src4 = [ppval(pp4, k2(1:ik2)), ppval(pp4, k2(ik2+1:end))];
   src5 = [ppval(pp5, k2(1:ik2)), ppval(pp5, k2(ik2+1:end))];
end

% calculate the anisotropy spectrum for l <= 100
if (one_percent == 1)  &  (K > 0 )
   lstep = 20;
else
   lstep = 10;
end
%la1 = 2;
if (K > 0)
   la1 = [2,6,10, 20:lstep:220];
   lla1 = size(la1,2);
   if one_percent == 1
      [cct1, ta1, te1] = srcintf1(la1, nn2, ct4, src4, src5, ctc0, ctad, nprimtilt, lmax, K, stfname);
   else
      [cct1, ta1, ta2] = srcint1(la1, nn2, ct4, pp4, pp5, ctc0, ctad, nprimtilt, lmax, K, stfname);
   end

else
   la1 = [2,6,10,20:lstep:90];
   lla1 = size(la1,2);
   if one_percent == 1
      [cct1, ta1, te1] = srcintf(la1, k2, ct4, src4, src5, ctc0, nprimtilt, lmax, K, stfname, logstep);
   else
      [cct1, ta1, ta2] = srcint(la1, k2, ct4, pp4, pp5, ctc0, nprimtilt, lmax, K, stfname, logstep);
   end

end

% change the stepsize in the k2 vector
if (K > 0)
   nksteps = 0.75* sqrt(lmax * lmax0);
% take care that in case the peak in the source region is splitup by the symmetry point, we take all
% modes (even and uneven nn2(i)) in the summation for the anistropy spectrum and nn2min==3
   if  ( (ctc0 - ct4sym) > 0 &  (ctc0 - ct4sym) < (ctad + 100) )
      nn2step = 1;
   else
      nn2step = max(1, round((nn2max - nn2min)/nksteps));
      nn2step = min(2, nn2step);   % limit to nn2step==2
   end
   nn2 = nn2min : nn2step : nn2max;
   k2 = nn2 * sqrt(K);
else

   clear('k2')
   logstep = 0;
   if one_percent == 1
     nksteps = 1.5*sqrt(lmax *lmax0);
   else
     nksteps = 2*sqrt(lmax *lmax0);
   end
   k2step = (k2max - k2min)/nksteps;
   k2 = k2min : k2step : k2max;
end % define k2 range
lk2 = size(k2, 2);
ik2 = floor(lk2/2);
ik4 = floor(lk2/4);

%##anisotropy spectrum, step 2 :l > 100, note that the cftime extends to the full range, ict4max= lt2-1
if one_percent == 1
   ict4step = 2;
else
   ict4step = 1;
end
ict4 = ict4min : ict4step : ict4max;
ct4 = ct2(ict4);
lct4 = size(ict4, 2);

% take the splines for all k2 interpolations within the conformal time interval ict4
% for each time we have one spline (the splines are in a column vector)
% this steps needs a lot of memory : if short of memory use srcpint.m with pp4, pp5 as argument
pp4 = spline(k1, table_srt(:, ict4)');
pp5 = spline(k1, table_srp(:, ict4)');
if one_percent == 1
   src4 = [ppval(pp4, k2(1:ik2)), ppval(pp4, k2(ik2+1:end))];
   src5 = [ppval(pp5, k2(1:ik2)), ppval(pp5, k2(ik2+1:end))];
end

% calculate the anisotropy spectrum for l > 100
if one_percent == 1
   lstep = 50;
else
   lstep = 25;
end
%la2 = 500
if (K > 0)
   la2 = 250:lstep:lmax;
   lla2 = size(la2,2);
   if  one_percent == 1
      [cct2, ta1] = srcintf1(la2, nn2, ct4, src4, src5, ctc0, ctad, nprimtilt, lmax, K,stfname);
   else
      [cct2, ta1] = srcint1(la2, nn2, ct4, pp4, pp5, ctc0, ctad, nprimtilt, lmax, K,stfname);
   end
else
   la2 = 100:lstep:lmax;
   lla2 = size(la2,2);
   if one_percent == 1
      [cct2, ta1] = srcintf(la2, k2, ct4, src4, src5, ctc0, nprimtilt, lmax, K, stfname, logstep);
   else
      [cct2, ta1] = srcint(la2, k2, ct4, pp4, pp5, ctc0, nprimtilt, lmax, K, stfname, logstep);
   end
end

% collect the results and go to the final step
la = [la1, la2];
ctth = [cct1.ctt, cct2.ctt];
ceeh = [cct1.cee, cct2.cee];
cteh = [cct1.cte, cct2.cte];

% calculate the amplitude of the anistropy spectrum in terms of the cmb temperature cmbtemp as
% cmbtemp^2*l*(l+1)*C1/(2*pi) in units of Kelvin^2 (including the effect of primordial tilt of the
% power spectrum)
% calculate the present (a==1) grow factor D1 by integration (see Dodelson 7.77)
D1 = 2.5 * om * quadl(@cmd3, 0, 1.0, tol, trace, om, omt);
% calculate the proportionality factor of the temperature spectrum, see Sachs and Wolfe and
% Dodelson (8.74) :
cswl1 = (50/9) * cmbtemp^2* (h*h0)^(1-nprimtilt) * (om/D1)^2;

% calculate the power spectrum of matter psm in units (Mpc/h)^3 as a function of kh = k1/h in h/Mpc.
kh = k1/h;
% use deltam(k) = delta(k)/psi, see Dodelson formulas (6.100) and (7.9) and delta is the mass overdensity.
psm1 =(50 *pi^2/9) .* (h./k1) .^3 .* (k1/(h*h0)) .^ (nprimtilt-1) .* (deltam .^2) * (om/D1)^2;
deltam(1)

% find psm at k1=0.05/Mpc and normalise the spectra with the density fluctuation at horizon crossing (dha)
% the c10 normalisation for the temperature spectrum is taken at l=10 :
% the best fit value for wmap2003 is 727 microKelvin^2, which gives C10 = 5.59e-12 and sqrt(C10) = 2.36e-6.
% the density fluctuation at the horizon for the best fit (dha) is in the global file of constants (cmbglobl.m)
% we use it to normalise both temperature and matter spectra.
kfid = 0.05
% normalise the temperature and matter spectra:
psm = psm1 * dha2;
cswl = cswl1 * dha2;

% calculate sigma8, the extended rms mass overdensity, sampled with a sphere of radius of radius 8 Mpc/h
% extend the matter power spectrum from k1 into small scales with k3 > 0.2 Mpc^-1,
k3step = 0.01;
k3max = 1.0;
k3 = (k1max+k3step) : k3step : k3max;
%for the extension use the asymptotic expansion of delta, see Dodelson formula (7.9) with the transfer function of (7.69)
% as psm is in units (Mpc/h)^3, we have to include the Hubble factor^3 in the Bolzman part of psm
psm2 = psm/h^3;
psm3 = psm2(end) * (k1(end) ./ k3) .^(4-nprimtilt) .* (log(k3/(8*keq))/log(k1(end)/(8*keq))) .^2;
k4 = [k1, k3];
psm4 = [psm2, psm3];
sig8h = (0.5/pi^2) * quadl(@sigma_8, k4(1), k4(end), tol, trace, k4, psm4, h);
sig8 = sqrt(sig8h);
% plot both psm from Bolzman calculation and the extension from formula (7.69)
figure(1)
loglog(k4, psm4)
%semilogy(k4, psm4)
title(['full range mass power spectrum: omt=', num2str(omt), ' h=', num2str(h), ',om=', num2str(om), ...
   ',omb=', num2str(omb), ',ns=', num2str(nprimtilt), ',tauri=', num2str(tauri)])
ylabel('P(k) in (Mpc/h)^3')
xlabel('log of the wavenumber in h/Mpc')

% compare with wmap results and find the normalisation constant c10
fn = 'cmbansp';
pm.k = kh;
pm.psm = psm;
as.la = la;
as.ctt = ctth * cswl;
as.cee = ceeh * cswl;
as.cte = cteh * cswl;
% c1==1 gives a model with a perfect fit

[c1, lmin1, lmax1] = cmbdata1(as, pm, h, omt, om, omb, nprimtilt, tauri, bmark);
% sigma* has to be corrected with sqrt of the amplitude correction ,sqrt(c1) :standard case sig8 =0.915
sig8 = sig8 * sqrt(c1);

% save the anisotropy spectra (as/pm) and matter spectrum (pm) and make plots
%save cmbansp as pm
% bring the cosmological parameters and other data that have been used in the routine to the output:
fid = 1;
fprintf(fid, '\n');
fprintf(fid, 'the redshift at decoupling (maximum of the visibility function) = %g\n', dt.zdec);
fprintf(fid, 'the thickness of the decoupling surface (half width of vsb in redshift) = %g\n', dt.tds);
fprintf(fid, 'the age of the universe at decoupling = %g kyear\n', dt.tdec);
fprintf(fid, 'the age of the universe at reionisation = %g kyear\n', dt.tri);
fprintf(fid, 'the current age of the universe = %g Gyear\n', dt.trc);
fprintf(fid, '\n');
fprintf(fid, 'the conformal distance at decoupling = %g Mpc\n', dt.etadec);
fprintf(fid, 'the thickness of the decoupling surface (half width of vsb in redshift) = %g Mpc\n', trs);
fprintf(fid, 'the number of time integration points over the thickness of the decoupling surface = %g\n', 2*nphw);
fprintf(fid, '\n');
fprintf(fid, 'the cosmological parameter values used in the routine are :\n');
fprintf(fid, 'the density fluctuation at the horizon = %g\n', dha);
fprintf(fid, 'the amplitude correction (at l=10) = %g\n', c1);
fprintf(fid, 'total omega (omt) = %g\n', omt);
fprintf(fid, 'dark matter omega = %g\n', om);
fprintf(fid, 'baryon omega = %g\n', omb);
fprintf(fid, 'Hubble factor = %g\n', h);
fprintf(fid, 'spectral density index = %g\n', nprimtilt);
fprintf(fid, 'optical depth at reionisation = %g\n', tauri);
fprintf(fid, 'sigma8, rms mass overdensity sampled with a sphere of 8 Mpc/h = %g\n', sig8);
fprintf(fid, ['in the figures we compare the results with ', bmark, '\n'])
fprintf(fid, 'the minima in the temperature spectrum are at: ')
fprintf(fid, '%g,  ', lmin1)
fprintf(fid, '\n')
fprintf(fid, 'the maxima in the temperature spectrum are at: ')
fprintf(fid, '%g,  ', lmax1)
fprintf(fid, '\n')


text:
%SETCOST Reset classification cost matrix of dataset
%
%   A = SETCOST(A,COST,LABLIST)
%
% The classification cost matrix of the dataset A is reset to COST.
% COST should have size [C,C+n], n >= 0, if C is the number of classes.
% COST(I,J) are the costs of classifying an object of class I
% as class J. Columns C+j generate an alternative reject classes and
% may be omitted, yielding a size of [C,C].
% An empty cost matrix, COST = [] (default) is interpreted as
% COST = ONES(C) - EYE(C) (identical costs of misclassification).
%
% In LABLIST the corresponding class labels may be supplied.
% LABLIST may have only class names of the existing classes in A.
% Reset class names first by SETLABLIST if necessary.
%
% Alternatively, for classification matrices, LABLIST may refer to
% the class names stored in the feature labels. This should be used
% with care as it may disturb the existing lableling of A.
%
% If LABLIST is not given, the order defined by the existing LABLIST
% for A (determined by [NLAB,LABLIST] = renumlab(LABELS)) is used.

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [min] = us2min(us)
% Convert time from microseconds to minutes.
% Chad Greene 2012
min = us*1.666666666667e-8;

author: Sable

repo: mcbench-benchmarks


text:
function [x,fval,exitflag,info] = opti_sedumi(At,b,c,K,opts)
%OPTI_SEDUMI Solve a SDP using SEDUMI
%
%   x = opti_sedumi(At,b,c,K) solves a LP or SDP in SeDuMi format.
%
%   x = opti_sedumi(At,b,c,K,opts) uses opts to pass optiset options to the
%   solver.
%
%   [x,fval,exitflag,info] = opti_sedumi(...) returns the objective value at
%   the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR SEDUMI (GPL)

%   Copyright (C) 2013 Jonathan Currie (IPL)

if(nargin < 5), opts = optiset; end
if(nargin < 4), K = []; end
if(nargin < 3), error('This function expects at least 3 arguments!'); end

%Set SEDUMI options
pars = opts.solverOpts;
pars.fid = min(dispLevel(opts.display),1);
pars.maxiter = opts.maxiter;
pars.eps = opts.tolrfun;

t = tic;
%Run Solver
[X,x,sinfo] = sedumi(At,b,c,K,pars);

%Assign primal fval
fval = -b'*x;

%Assign Outputs
info.Iterations = sinfo.iter;
info.Time = toc(t);
info.Algorithm = 'SeDuMi: Self-Dual-Minimization SDP Solver';

%Assign Exit Flag
if(sinfo.iter >= opts.maxiter)
    info.Status = 'Exceeded Maximum Iterations';
    exitflag = 0;
elseif(sinfo.numerr > 0)
    info.Status = 'Numerical Errors';
    exitflag = -2;
elseif(sinfo.pinf==1 || sinfo.dinf==1)
    info.Status = 'Infeasible';
    exitflag = -1;
elseif(sinfo.pinf==0 && sinfo.dinf==0)
    info.Status = 'Optimal';
    exitflag = 1;
end
info.DualObjective = c'*X;
info.X = X;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function dtiWriteFibersPdb(fg, xformToAcpc, pdbFname)
%
% dtiWriteFibersPdb([fg=uigetfile], [xformToAcpc=eye(4)], [[pdbFname=uiputfile])
%
% Writes the fiber group to disk in the specified file using the PDB format
% recognized by CINCH.
% * Note that .pdb extension will be added to pdbFname.
%
% Example usage:
% dtiWriteFibersPdb(fg,dt.xformToAcpc,fullfile(fiberDir,'Lfibers'));
%
% Note that CINCH doesn't use the xformToAcpc, so you can safely skip it.
%
%
% HISTORY:
% 2008.07.02 LMP: wrote it.
% 2011.01.27 LMP: Now just uses mtrExport fibers.
%
%
disp([mfilename ' is now obsolete -- wrapping your inputs for mtrExportFibers...']);

if(~exist('fg','var')||isempty(fg))
    [fg,fgName] = dtiReadFibers;
    [p,f,e] = fileparts(fgName);
    fgName = fullfile(p,[f '.pdb']);
else
    fgName = fullfile(pwd,'fibers.pdb');
end

if(~exist('xformToAcpc','var')||isempty(xformToAcpc))
    xformToAcpc = [];
end

if(~exist('pdbFname','var')||isempty(pdbFname))
    [f,p] = uiputfile(fgName,'Save pdb file as...');
    if(isequal(f,0)), disp('user canceled.'); return; end
    pdbFname = fullfile(p,[f '.pdb']);
end
mtrExportFibers(fg, pdbFname);

return

% if(~exist('fg','var')||isempty(fg))
%     [fg,fgName] = dtiReadFibers;
%     [p,f,e] = fileparts(fgName);
%     fgName = fullfile(p,[f '.pdb']);
% else
%     fgName = fullfile(pwd,'fibers.pdb');
% end
%
% if(~exist('xformToAcpc','var')||isempty(xformToAcpc))
%     xform = eye(4);
% else
%     % The xform that CINCH/dtiQuery wants assumes 1mm voxels:
%     [t,r,s,k] = affineDecompose(xformToAcpc);
%     xform = affineBuild(t,r,sign(s),k);
% end
%
% if(~exist('pdbFname','var')||isempty(pdbFname))
%     [f,p] = uiputfile(fgName,'Save pdb file as...');
%     if(isequal(f,0)), disp('user canceled.'); return; end
%     pdbFname = fullfile(p,f);
% end
%
% % add the .pdb file extension
% if(~strcmpi(pdbFname(end-3:end),'.pdb'))
%     pdbFname = [pdbFname,'.pdb'];
% end
%
% dtiWritePDBHeader(xform, pdbFname);
% offsets = dtiAppendPathwaysToPDB(fg, pdbFname);
% dtiAppendFileOffsetsToPDB(offsets, pdbFname);
%
% return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% function vectorized_bspline_coeff
% ------
%
% Input
%  vi [n x m]
%  vs: [n x m]
%
%  see Kristin Branson's "A Practical Review of Uniform B-splines"
%
% Output
%   C [n x 1]: the coefficients
function C = vectorized_bspline_coeff(vi,vs)

    assert(isequal(size(vi),size(vs)));

    % Go through conditions
    C = zeros(size(vi));

    sel1 = vs >= vi & vs < vi+1;
    C(sel1) = (1/6)*(vs(sel1)-vi(sel1)).^3;

    sel2 = vs >= vi+1 & vs < vi+2;
    C(sel2) = (1/6)*(-3*(vs(sel2)-vi(sel2)-1).^3 + 3*(vs(sel2)-vi(sel2)-1).^2 + 3*(vs(sel2)-vi(sel2)-1)+1);

    sel3 = vs >= vi+2 & vs < vi+3;
    C(sel3) = (1/6)*(3*(vs(sel3)-vi(sel3)-2).^3 - 6*(vs(sel3)-vi(sel3)-2).^2 + 4);

    sel4 = vs >= vi+3 & vs < vi+4;
    C(sel4) = (1/6)*(1-(vs(sel4)-vi(sel4)-3)).^3;

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Torr Vision Group (TVG)
%  University of Oxford - UK
%
%  Qizhu Li <liqizhu@robots.ox.ac.uk>
%  August 2018
% ------------------------------------------------------------------------
% This file is part of the weakly-supervised training method presented in:
%    Qizhu Li*, Anurag Arnab*, Philip H.S. Torr,
%    "Weakly- and Semi-Supervised Panoptic Segmentation,"
%    European Conference on Computer Vision (ECCV) 2018.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
% This function removes thing-class predictions outside their bounding
% boxes.
%  INPUT:
%  - pred_label : prediction
%  - bbox_masks : object bounding box masks produced by make_bbox_masks.m
%  - thing_classes : list of train ids for thing classes
%  - ignore_label  : label to ignore, normally 255
%
%  OUTPUT:
%  - modified_label: processed prediction
%
%  DEMO:
%  - See scripts/clean_label.m
% ------------------------------------------------------------------------

function modified_label = apply_bbox_prior(pred_label, bbox_masks, thing_classes, ignore_label)
modified_label = pred_label;
for j = reshape(thing_classes, 1, [])
    if isempty(bbox_masks{j+1})
        modified_label(modified_label == j) = ignore_label;
    else
        prediction_bbox_prior_mask = modified_label == j;
        gt_bbox_mask = bbox_masks{j+1};
        ignore_mask = and(prediction_bbox_prior_mask, ~gt_bbox_mask);
        modified_label(ignore_mask) = ignore_label;
    end
end
end

author: qizhuli

repo: Weakly-Supervised-Panoptic-Segmentation


text:
function disp( x, prefix )
if nargin < 2,
    prefix = '';
end
str = [ prefix, 'cvx: ', cvx_class( x, true, true, true ) ];
if numel( x ) > 1,
    str = [ str, ' ', type( x ) ];
elseif any( strcmp( str(end-5:end), { 'nstant', ': zero' } ) ),
    str = [ prefix, 'cvx: ', num2str(cvx_constant(x)) ];
end
disp( str );
dual = cvx_getdual( x );
if ~isempty( dual ) && isstruct( dual ),
    disp( [ prefix, '   tied to dual variable: ', dual.subs ] );
end

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

author: yu-jiang

repo: radpbook

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


text:
%% Zero Padding Function
% Input data is the sequence s, Target length N
% output:
% zs :zero padded signal
% nozs = number of zeros
% [nozs zs] = zpp(s, N)
% 22/04/2011
function [nozs zs]=zpp(s, N)
L=length(s);
sz=zeros(N,1);
nozs=N-L;
sz(1:L/2)=s(1:L/2);

sz(end-L/2+1:end)=s(end-L/2+1:end);
zs=sz;


text:
function matrix = rmAverageTime(matrix,nrep);
% rmAverageTime - average non-unique epochs in time (1st) dimension
%
%  out = rmAverageTime(in,nrep);
%
% 2006/03 SOD: wrote it.

% sanity check (<1 no averaging needed)
if nrep <= 1,
  return;
else,
  matrixin = matrix;
end;

% get total size input
sz  = size(matrixin);

% get total size output
len      = sz(1)./nrep;
szout    = sz;
szout(1) = len;

% initiate matrixout
matrix   = matrixin(1:len,:);

% repeat (add) process
start = len;
for n=1:nrep-1,
    matrix = matrix + matrixin(start+1:start+len,:);
    start     = start + len;
end;

% mean
matrix = matrix ./ nrep;

% reshape if necesary
matrix = reshape(matrix,szout);

return


text:
classdef noutput < dml.method
%NOUTPUT wrapper class to make methods handle multiple outputs.
%
%   DESCRIPTION
%
%
%   REFERENCE
%
%
%   EXAMPLE
%
%
%   DEVELOPER
%   Marcel van Gerven (m.vangerven@donders.ru.nl)

  properties

    method  % the method which should be replicated

  end

  methods

    function obj = noutput(varargin)

     obj@dml.method(varargin{:});

     if ~iscell(obj.method)
       obj.verbose = obj.method.verbose;
     else
       obj.verbose = obj.method{1}.verbose;
     end

    end

    function obj = train(obj,X,Y)

      ny = size(Y,2);

      if ~iscell(obj.method)
        obj.method = repmat({obj.method},[1 ny]);
      end

      for c=1:ny

        if obj.method{c}.verbose
          fprintf('training output %d of %d\n',c,ny);
        end

        obj.method{c} = obj.method{c}.train(X,Y(:,c));
      end

    end

    function Y = test(obj,X)

      ny = length(obj.method);

      Y = [];
      for c=1:ny

        if obj.method{c}.verbose
          fprintf('testing output %d of %d\n',c,ny);
        end

        Y = [Y obj.method{c}.test(X)];

      end

    end

    function m = model(obj)

      ny = length(obj.method);

      m = cell(1,ny);
      for c=1:ny

        m{c} = obj.method{c}.model;

      end

    end

  end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Automatic flight control of F8 aircraft using model predictive control
% and various models

% MPC execution time for sine3 models
% 83.5243, 188.2009, 2.1348e+03
% 79.5731, 179.1829, 2.1387e+03

clear all, close all, clc

SystemModel = 'F8';

figpath = ['../FIGURES/',SystemModel,'/'];mkdir(figpath)
datapath = ['../DATA/',SystemModel,'/'];mkdir(datapath)
addpath('../utils');

%% Load Models
Nvar = 3;
InputSignalTypeModel = 'sine3';             % Use actuation to collect training data
InputSignalTypeModel_Validation = 'sine2';  % Use actuation to collect validation data to test for generalization

% DelayDMDc
ModelTypeDMDc = 'DMDc'; % Model: DelayDMDc , DMDc
load(fullfile(datapath,['EX_',SystemModel,'_SI_',ModelTypeDMDc,'_',InputSignalTypeModel,'.mat']))
Models.DelayDMDc = Model;

% SINDYc
load(fullfile(datapath,['EX_',SystemModel,'_SI_SINDYc','_',InputSignalTypeModel,'.mat']))
Models.SINDYc = Model;

% NARX
NARX_SUBSTRACT_MEAN = 0;
load(fullfile(datapath,['EX_',SystemModel,'_SI_NARX','_',InputSignalTypeModel,'.mat']))
Models.NARX = Model;

%% Get training data
ENSEMBLE_DATA = 0;
ONLY_TRAINING_LENGTH = 1;
InputSignalType = InputSignalTypeModel;
Ndelay = Models.DelayDMDc.Ndelay;
getTrainingData

%% Get validation data: Model comparison
InputSignalType = InputSignalTypeModel_Validation;
Ndelay = Models.DelayDMDc.Ndelay;
tspanV =[tspan(end):dt:150];

% Forcing
if strcmp(InputSignalTypeModel_Validation,'sine2')
    forcing = @(x,t) [(.05* (sin(0.7*t).*sin(.1*t).*sin(.2*t).*sin(.05*t)) )];
elseif strcmp(InputSignalTypeModel_Validation,'sine3')
    forcing = @(x,t) (.5*sin(5*t).*sin(.5*t)+0.1).^3;
end

% Reference
[tA,xA] = ode45(@(t,x)F8Sys(t,x,forcing(x,t)),tspanV,x(end,1:3),options);   % true model
u_valid = forcing(0,tA)';

% DelayDMDc / DMDc
if Ndelay == 1
    x0      = [x(end,1:Nvar)];
    Hunew   = [u(end),u_valid(1:end-1)];
    [xB,tB] = lsim(Models.DelayDMDc.sys,Hunew,tspanV,[x0-[Models.DelayDMDc.xmean]]');
elseif Ndelay > 1
    x0      = [x(end-Ndelay+1,1:Nvar),x(end,1:Nvar)];
    Hunew   = [ u(end-Ndelay+1:end),u_valid(1:end-Ndelay);
        u(end),u_valid(1:end-1)];
    [xB,tB] = lsim(Models.DelayDMDc.sys,Hunew,tspanV,[x0-[Models.DelayDMDc.xmean]]');
    xB = xB(:,Nvar+1:2*Nvar); xB = xB + repmat(xmean,[size(xB,1) 1]);
end
xB = xB + repmat(Models.DelayDMDc.xmean,[length(tB) 1]);

% SINDYc
[tC,xC]=ode45(@(t,x)sparseGalerkinControl(t,x,forcing(x,t),Models.SINDYc.Xi(:,1:Nvar),Models.SINDYc.polyorder,Models.SINDYc.usesine),tspanV,x(end,1:Nvar),options);  % approximate

% NARX
Udummy = [u(end),u_valid(1:end)];
Hdummy = zeros(Nvar,size(Udummy,2));
if NARX_SUBSTRACT_MEAN == 1
    NARX_xmean = Models.NARX.xmean;
else
    NARX_xmean = zeros(size(Models.NARX.xmean'));
end

Hdummy(:,1) = x(end,1:Nvar)'-NARX_xmean';
[Us,Ui,Si] = preparets(Models.NARX.net,con2seq(Udummy),{},con2seq(Hdummy));
xD = Models.NARX.net(Us,Ui,Si); % Predict on validation data
xD = cell2mat(xD)'; xD = [x0;xD(1:end-1,:)];
xD(2:end,:) = xD(2:end,:) + repmat(NARX_xmean,[size(xD,1)-1 1]);

%% Show results
clear ph
h = figure;
subplot(3,1,1), box on, hold on
plot([tB(1) tB(1)],[-100 100],'--k')
plot(t,x(:,1),'Color',[.4 .4 .4],'LineWidth',1.5);
plot(tA,xA(:,1),'k','LineWidth',1.5);
plot(tB,xB(:,1),'r-','LineWidth',1.5);
plot(tC,xC(:,1),'g--','LineWidth',1.5);
plot(tA(1:end),xD(:,1),'-.','Color',[0.7,0.7,1],'LineWidth',1.5);
grid on
ylim([-0.4 0.3])
ylabel('x_1','FontSize',13)
set(gca,'FontSize',13)

subplot(3,1,2), box on, hold on
plot([tB(1) tB(1)],[-100 100],'--k')
plot(t,x(:,2),'Color',[.4 .4 .4],'LineWidth',1.5);
plot(tA,xA(:,2),'k','LineWidth',1.5);
plot(tB,xB(:,2),'r-','LineWidth',1.5);
plot(tC,xC(:,2),'g--','LineWidth',1.5);
plot(tA(1:end),xD(:,2),'-.','Color',[0.7,0.7,1],'LineWidth',1.5);
ylim([-4 0.1])
set(gca,'FontSize',13)
ylabel('x_2','FontSize',13)

subplot(3,1,3), box on, hold on
plot([tB(1) tB(1)],[-100 100],'--k')
ph(1) = plot(t,x(:,3),'Color',[.4 .4 .4],'LineWidth',1.5);
ph(2) = plot(tA,xA(:,3),'k','LineWidth',1.5);
ph(3) = plot(tB,xB(:,3),'r-','LineWidth',1.5);
ph(4) = plot(tC,xC(:,3),'g--','LineWidth',1.5);
ph(5) = plot(tA(1:end),xD(:,3),'-.','Color',[0.7,0.7,1],'LineWidth',1.5);
l1=legend(ph,'Training','Validation',ModelTypeDMDc,'SINDYc','NARX');
set(l1,'Location','NorthWest')
grid on
ylim([-0.9 0.5])
ylabel('x_3','FontSize',13)
set(gca,'FontSize',13)
xlabel('Time','FontSize',13)

set(h,'Units','Inches');
set(gcf,'Position',[1 1 6. 5.5])
pos = get(h,'Position');
set(h,'PaperPositionMode','Auto','PaperSize',[pos(3), pos(4)])
print(h,'-painters','-depsc2',  '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_Comparison_Validation_',InputSignalTypeModel,'_',InputSignalType,'.eps'],'-r0');

%% Actuation signal
t_valid = tspanV;
clear ph
figure,box on, hold on,
ccolors = get(gca,'colororder');
plot([tA(1),tA(1)],[-15 260],':','Color',[0.4,0.4,0.4],'LineWidth',1.5)
plot(t,u,'-k','LineWidth',1);
plot(t_valid,u_valid,'-k','LineWidth',1);
grid off
ylim([min([u u_valid])+0.05*min([u u_valid]) max([u u_valid])+0.05*max([u u_valid])])
xlim([0 t_valid(end)])
xlabel('Time')
ylabel('Input')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_Comparison_Validation_',InputSignalTypeModel,'_',InputSignalType,'_Actuation.eps']);

%% True SINDYc model parameters
Xi0 = zeros(size(Models.SINDYc.Xi));
Xi0([2,4,5,6,8,10,16,18,19,25,35],1) = [-0.877,1,-0.215,0.47,-0.088,-0.019,3.846,-1,0.28,0.47,0.63];
Xi0([4],2) = 1;
Xi0([2,4,5,6,16,19,25,35],3) = [-4.208,-0.396,-20.967,-0.47,-3.564,6.265,46,61.4];

%% MPC
% Parameters MPC
options = optimoptions('fmincon','Algorithm','sqp','Display','none', ...
    'MaxIterations',100);
Duration = 6;                   % Run for 'Duration' time units
Ton = 0;                        % Time units when control is turned on
getMPCparams

x0n=[0.1 0 0]'; %xA(end,:)';    % Initial condition
Tcontrol = tA(end);             % Time offset to combine training, prediction and control phase

% Parameters Models
ModelCollection = {ModelTypeDMDc, 'SINDYc', 'NARX'}; % 51.8924, 46.0083, 1.4662e+03 execution times
Nmodels = length(ModelCollection);

clear Results
Results(1:Nmodels) = struct('x',[], 'u', [], 't', [], 'xref', [], 'J', [], 'elapsed_time', []);

Ts = Models.SINDYc.dt;

%%
for iM = 1:Nmodels
    select_model = ModelCollection{iM}; % DelayDMDc; DMDc; NARX ; SINDYc

    % Prepare variables
    Nt = (Duration/Ts)+1;
    uopt0    = 0;
    xhat     = x0n;
    uopt     = uopt0.*ones(Nu,1);
    xHistory = zeros(Nvar,Nt); xHistory(:,1) = xhat;
    uHistory = zeros(1,Nt); uHistory(1)   = uopt(1);
    tHistory = zeros(1,Nt); tHistory(1)   = Tcontrol;
    rHistory = zeros(1,Nt);

    % Parameters Model
    switch select_model
        case 'DelayDMDc'
            pest.dt     = Models.DelayDMDc.dt;
            pest.sys    = Models.DelayDMDc.sys;
            pest.xmean  = xmean';
            pest.udelay = zeros(1,N);
            pest.xdelay = [x(end-Ndelay+1,1:2)]';
            pest.Nxlim  = size(x,1);
        case 'DMDc'
            pest.dt = Models.DelayDMDc.dt;
            pest.sys = Models.DelayDMDc.sys;
            pest.xmean = Models.DelayDMDc.xmean';
            pest.Nxlim = size(x,1);
        case 'SINDYc'
            pest.ahat = Models.SINDYc.Xi(:,1:Nvar); % Replace with Xi0 to test true model
            pest.polyorder = Models.SINDYc.polyorder;
            pest.usesine = Models.SINDYc.usesine;
            pest.dt = Models.SINDYc.dt;
        case 'NARX'
            pest.net = Models.NARX.net;
            pest.Ndelay = 1;
            pest.udelay = zeros(1,Ndelay);
            pest.xdelay = zeros(2,length(pest.udelay));
            pest.xdelay(:,1:Ndelay) = [x(end-Ndelay+1:end,1:2)]';
            pest.Nxlim = size(x,1);
    end

    % Start simulation
    fprintf('Simulation started.  It might take a while...\n')
    tic
    for ct = 1:(Duration/Ts)

        % Set references over prediction horizon
        tref = (ct:ct+N-1).*Ts;
        xref = [xrefFUN(tref); zeros(1,N); zeros(1,N)];

        % NMPC with full-state feedback
        COSTFUN = @(u) ObjectiveFCN_models(u,xhat,N,Nu,xref,uHistory(:,ct),pest,diag(Q),R,Ru,select_model);
        CONSFUN = @(u) ConstraintFCN_models(u,uHistory(:,ct),xhat,N,LBo,UBo,LBdu,UBdu,pest,select_model);
        uopt = fmincon(COSTFUN,uopt,[],[],[],[],LB,UB,CONSFUN,options);

        % Run without constraints
%         uopt = fmincon(COSTFUN,uopt,[],[],[],[],[],[],[],options);

        % Integrate system
        xhat = rk4u(@F8Sys,xhat,uopt(1),Ts/10,10,[],0);
        xHistory(:,ct+1) = xhat;
        uHistory(:,ct+1) = uopt(1);
        tHistory(:,ct+1) = ct*Ts+Tcontrol;
        rHistory(:,ct+1) = xref(1,2);

    end
    tElapsed = toc
    fprintf('Simulation finished!\n')

    % Collect results
    Results(iM).eval_time = tElapsed;
    Results(iM).xref = rHistory;
    Results(iM).x = xHistory;
    Results(iM).u = uHistory;
    Results(iM).t = tHistory;
    Results(iM).J = evalObjectiveFCN(Results(iM).u,Results(iM).x,Results(iM).xref,diag(Q),R,Ru);
    Results(iM).elapsed_time = tElapsed;

    VIZ_SI_Validation_MPC

    if iM==Nmodels % Plot ensemble training for NN / NARX
       VIZ_SI_Validation_MPC_ensemble(ct,t_valid,u_valid,xHistory,tHistory,uHistory,Results,t,tA,xA,xB,xC,xD,select_model,SystemModel,figpath,N,InputSignalTypeModel,Nmodels,ModelCollection)
    end

end

% MPC Execution times using N=10: 120.9128, 96.8967, 2.1755e+03

%% Unforced
[tU,xU] = ode45(@(t,x)F8Sys(t,x,0,[]),tHistory,xA(end,1:3),options);   % true model

%% Show cost
clear ph
figure, hold on, box on
ph(1) = plot(Results(1).t,cumsum(Results(1).J),'-r','LineWidth',2);
ph(3) = plot(Results(1).t,cumsum(Results(3).J),'-.','Color',[0.7,0.7,1],'LineWidth',2);
ph(2) = plot(Results(1).t,cumsum(Results(2).J),'--g','LineWidth',2);

xlim([Results(1).t(1) Results(1).t(end)])
ylabel('Cost','FontSize',14)
xlabel('Time','FontSize',14)
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')

print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_Comparison_Validation_',InputSignalTypeModel,'_',InputSignalType,'_Cost','_N_',num2str(N),'.eps']);

%% Show control stage separate
clear ph
figure, hold on, box on
ph(4) = plot(Results(1).t(2:end),(Results(1).xref(2:end)),'-k','LineWidth',2);
ph(1) = plot(Results(1).t,(Results(1).x(1,:)),'-r','LineWidth',2);
ph(3) = plot(Results(1).t,(Results(3).x(1,:)),'-.','Color',[0.7,0.7,1],'LineWidth',2);
ph(2) = plot(Results(1).t,(Results(2).x(1,:)),'--g','LineWidth',2);

l1=legend(ph,ModelCollection,'Ref');
set(l1,'Location','NorthEast')
xlim([Results(1).t(1) Results(1).t(end)])
ylim([-0.25 .25])
ylabel('AoA','FontSize',14)
xlabel('Time','FontSize',14)
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')

print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_Comparison_',InputSignalTypeModel,'_',InputSignalType,'_AoA','_N_',num2str(N),'.eps']);

clear ph
figure, hold on, box on
ph(1) = plot(Results(1).t,(Results(1).u(:)),'-r','LineWidth',2);
ph(3) = plot(Results(1).t,(Results(3).u(:)),'-.','Color',[0.7,0.7,1],'LineWidth',2);
ph(2) = plot(Results(1).t,(Results(2).u(:)),'--g','LineWidth',2);

xlim([Results(1).t(1) Results(1).t(end)])
ylabel('Input','FontSize',14)
xlabel('Time','FontSize',14)
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto')

print('-depsc2', '-loose','-cmyk', [figpath,'EX_',SystemModel,'_SI_Comparison_',InputSignalTypeModel,'_',InputSignalType,'_Input','_N_',num2str(N),'.eps']);

%% Save results
save(fullfile(datapath,['MPC_',SystemModel,'_Results_AllModels_',InputSignalTypeModel,'_',InputSignalType,'_N_',num2str(N),'.mat']),'Results')
save(fullfile(datapath,['MPC_',SystemModel,'_Results_AllModels_',InputSignalTypeModel,'_',InputSignalType,'_All','_N_',num2str(N),'.mat']))

author: eurika-kaiser

repo: SINDY-MPC


text:
function inplane = flat2ipAllROIs(flat,inplane,volume)
%
% inplane = flat2ipAllROIs(flat,inplane,[volume])
%
% Author: AAB/BW
%   Transform ROIs from flat view to inplane, going through the volume
%   representation.
%
% Examples:
%    INPLANE{1} = flat2ipAllROIs(FLAT{4},INPLANE{1});
%

if ieNotDefined('volume'), volume = initHiddenGray; end

[volume,volRoiIndexList] = flat2volAllROIs(flat,volume);
inplane = vol2ipAllROIs(volume,inplane);

return;


text:
function bd = blk_diag(A,n)
%BLK_DIAG Make or extract a sparse block diagonal matrix
% function bd = blk_diag(A,n);
% If A is not sparse, then
% returns a sparse block diagonal "bd", diagonalized from the
% elements in "A".
% "A" is ma x na, comprising bdn=(na/"n") blocks of submatrices.
% Each submatrix is ma x "n", and these submatrices are
% placed down the diagonal of the matrix.
%
% If A is already sparse, then the operation is reversed, yielding a block
% row matrix, where each set of n columns corresponds to a block element
% from the block diagonal.
%
% Routine uses NO for-loops for speed considerations.

% Copyright (c) 1993-1995, The Regents of the University of California.
% This software was produced under a U.S. Government contract
% (W-7405-ENG-36) by Los Alamos National Laboratory, which is operated
% by the University of California for the U.S. Department of Energy,
% and was funded in part by NIH grant R01-MH53213 through the University
% of Southern California to Los Alamos National Laboratory,
% and was funded in part by NIH grant R01-EY08610 to Los Alamos
% National Laboratory.
% The U.S. Government is licensed to use, reproduce, and distribute this
% software.  Permission is granted to the public to copy and use this
% software without charge, provided that this Notice and any statement
% of authorship are reproduced on all copies.  Neither the Government
% nor the University makes any warranty, express or implied, or assumes
% any liability or responsibility for the use of this software.
%
% Author: John C. Mosher, Ph.D.
% Los Alamos National Laboratory
% Group ESA-MT, MS J580
% Los Alamos, NM 87545
% email: mosher@LANL.Gov

% July 29, 1993 Author
% September 28, 1993 JCM Conversion to sparse
% July 27, 1995 JCM inverse block diagonal added

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

if(~issparse(A)),		% then make block sparse
    [ma,na] = size(A);
    bdn = na/n; 			% number of submatrices

    if(bdn - fix(bdn)),
        error('Width of matrix must be even multiple of n');
    end

    if(0)
        i = [1:ma]';
        i = i(:,ones(1,n));
        i = i(:); 			% row indices first submatrix

        ml = length(i); 		% ma*n

        % ndx = [0:(bdn-1)]*ma; 	% row offsets per submatrix
        ndx = [0:ma:(ma*(bdn-1))]; 	% row offsets per submatrix

        i = i(:,ones(1,bdn)) + ndx(ones(ml,1),:);
    else
        tmp = reshape([1:(ma*bdn)]',ma,bdn);
        i = zeros(ma*n,bdn);
        for iblock = 1:n,
            i((iblock-1)*ma+[1:ma],:) = tmp;
        end
    end

    i = i(:); 			% row indices foreach sparse bd

    j = [1:na];
    j = j(ones(ma,1),:);
    j = j(:); 			% column indices foreach sparse bd

    bd = sparse(i,j,A(:));

else 				% already is sparse, unblock it

    [mA,na] = size(A);		% matrix always has na columns
    % how many entries in the first column?
    bdn = na/n;			% number of blocks
    ma = mA/bdn;			% rows in first block

    % blocks may themselves contain zero entries.  Build indexing as above
    if(0)
        i = [1:ma]';
        i = i(:,ones(1,n));
        i = i(:); 			% row indices first submatrix

        ml = length(i); 		% ma*n

        % ndx = [0:(bdn-1)]*ma; 	% row offsets per submatrix
        ndx = [0:ma:(ma*(bdn-1))]; 	% row offsets per submatrix

        i = i(:,ones(1,bdn)) + ndx(ones(ml,1),:);
    else
        tmp = reshape([1:(ma*bdn)]',ma,bdn);
        i = zeros(ma*n,bdn);
        for iblock = 1:n,
            i((iblock-1)*ma+[1:ma],:) = tmp;
        end
    end

    i = i(:); 			% row indices foreach sparse bd

    if(0)
        j = [1:na];
        j = j(ones(ma,1),:);
        j = j(:); 			% column indices foreach sparse bd

        % so now we have the complete two dimensional indexing. Convert to
        % one dimensional

        i = i + (j-1)*mA;
    else
        j = [0:mA:(mA*(na-1))];
        j = j(ones(ma,1),:);
        j = j(:);

        i = i + j;
    end

    bd = full(A(i)); 	% column vector
    bd = reshape(bd,ma,na);	% full matrix
end

author: brainstorm-tools


text:
function y = trapmf(x, params)
%TRAPMF Trapezoidal membership function.
%   TRAPMF(X, PARAMS) returns a matrix which is the trapezoidal
%   membership function evaluated at X. PARAMS = [A B C D] is a 4-element
%   vector that determines the break points of this membership function.
%   We require that A <= B and C <= D. If B >= C, this membership
%   function becomes a triangular membership function that could have
%   a height less than unity. (See the example below.)
%
%   For example:
%
%       x = (0:0.1:10)';
%       y1 = trapmf(x, [2 3 7 9]);
%       y2 = trapmf(x, [3 4 6 8]);
%       y3 = trapmf(x, [4 5 5 7]);
%       y4 = trapmf(x, [5 6 4 6]);
%       plot(x, [y1 y2 y3 y4]);
%       set(gcf, 'name', 'trapmf', 'numbertitle', 'off');
%
%   See also DSIGMF, EVALMF, GAUSS2MF, GAUSSMF, GBELLMF, MF2MF, PIMF, PSIGMF,
%   SIGMF, SMF, TRIMF, ZMF.

%   Roger Jang, 6-28-93, 10-5-93, 4-14-94.
%   Copyright 1994-2002 The MathWorks, Inc.
%   $Revision: 1.22 $  $Date: 2002/04/14 22:21:13 $

if nargin ~= 2
    error('Two arguments are required by the trapezoidal MF.');
elseif length(params) < 4
    error('The trapezoidal MF needs at least four parameters.');
end

a = params(1); b = params(2); c = params(3); d = params(4);

if a > b,
    error('Illegal parameter condition: a > b');
elseif c > d,
    error('Illegal parameter condition: c > d');
end

y1 = zeros(size(x));
y2 = zeros(size(x));

% Compute y1
index = find(x >= b);
if ~isempty(index),
    y1(index) = ones(size(index));
end
index = find(x < a);
if ~isempty(index),
    y1(index) = zeros(size(index));
end
index = find(a <= x & x < b);
if ~isempty(index) & a ~= b,
    y1(index) = (x(index)-a)/(b-a);
end

% Compute y2
index = find(x <= c);
if ~isempty(index),
    y2(index) = ones(size(index));
end
index = find(x > d);
if ~isempty(index),
    y2(index) = zeros(size(index));
end
index = find(c < x & x <= d);
if ~isempty(index) & c ~= d,
    y2(index) = (d-x(index))/(d-c);
end

% Compute y
y = min(y1, y2);

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [imgThr, imgEb] = WardComputeThreshold(img, wardPercentile, wardTolerance)
%
%
%       [imgThr, imgEb] = WardComputeThreshold(img, wardPercentile, wardTolerance)
%
%       This function computes the Ward's MTB.
%
%       Input:
%           -img: an input image
%           -wardPercentile: a value for thresholding the image. This is
%           typically set to 0.5.
%           -wardTolerance: a tolerance threshold for classifying pixels
%           falling around edges.
%
%       Output:
%           -imgThr: Ward's threshold image. This image is set to 1 if the
%           the pixel value is greater or equal to the median value.
%           -imgEb: a tolerance mask of pixels around edges of imgThr.
%
%     Copyright (C) 2012  Francesco Banterle
%
%     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/>.
%

if(~exist('wardPercentile', 'var'))
    wardPercentile = 0.5;
end

if(~exist('wardTolerance', 'var'))
    wardTolerance = 4 / 256;
end

if(size(img, 3) == 1)
    grey = img;
else
    grey = (54 * img(:,:,1) + 183 * img(:,:,2) +  19 * img(:,:,3)) / 256;
end

medVal = MaxQuart(grey, wardPercentile);

imgThr = zeros(size(grey));
imgThr(grey > medVal) = 1.0;

A = medVal - wardTolerance;
B = medVal + wardTolerance;
imgEb = ones(size(grey));
imgEb((grey >= A) & (grey <= B)) = 0.0;

end

author: banterle


text:
function cout=filterbank(cin,h,g,order,basis)

%tstoolbox/@core/filterbank
%   Syntax:
%     * filterbank(cin,H,G,ORDER,BASIS)
%
%   Input Arguments:
%     * H - lowpass filter
%     * G - highpass filter
%     * ORDER - indicates the type of tree:
%          + 0 - band sorting according to the filter bank
%          + 1 - band sorting according to the frequency decomposition
%     * BASIS - desired subband decomposition
%
%   calculates the Wavelet Packet Transform of cin. It can be obtained
%   using a selection algorithm function. It may be switched from one
%   format to another using CHFORMAT. The different bands are sorted
%   according to ORDER and BASIS. If BASIS is omitted, the output is a
%   matrix with the coefficients obtained from all the wavelet packet
%   basis in the library. Each column in the matrix represents the outputs
%   for a level in the tree. The first column is the original signal. If
%   the length of X is not a power of 2, the columns are zero padded to
%   fit the different lengths. Run the script 'BASIS' for help on the
%   basis format.
%   See also: IWPK, CHFORMAT, PRUNEADD, PRUNENON, GROWADD, GROWNON.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

%       Auxiliary input argument for recursion:
%       flag: if set to 1, indicates that no basis is specified
%          and the output will be a matrix
%       flag2: if set to 1 the function was called with the detail signal,
%          then the coefficients must be rearranged

cout = core(wpk(data(cin), h,g,order,basis));

function w = wpk(x, h,g,order,basis,flag,flag2)

if ((nargin==4)|(nargin==5))
        if nargin==4
                L=floor(log2(length(x)));       % initialization when
                basis=ones(2^L,1)*(L-1);        % no basis was specified
                flag=1;
                flag2=0;
        else                                    % basis was specified
                if (~basis)
                        w=x;
                        return
                end
                basis=basis-1;
                flag=0;
                flag2=0;
        end
        if size(x,2)~=1                 % shapes the signal as column
                x=x(:);                 % vector if necessary
                fil=1;
        else
                fil=0;
        end
end

wx=wt(x,h,g,1);                         % perform one analysis level
                                        % into the analysis tree

N=length(wx);                           % separate approximation and detail
a=wx(1:N/2);                            % at the analysis output
d=wx(N/2+1:N);

if (flag2*order)                        % if the recursion is going to be
        s1=d;s2=a;                      % performed on the detail signal,
else                                    % and it is frequency sorted, the
        s1=a;s2=d;                      % outputs must be swapped
end

if all(basis==0)                           % two ending nodes achieved
        w=[s1;s2];
        if (nargin==5)&(fil), w=w'; end
        return
end;

tope=1;                                 % finds the point where the
suma=2^(-basis(1));                     % basis vector must be divided
i=1;                                    % in order to call the recursion
while (suma<tope)
        i=i+1;
        suma=suma+2^(-basis(i));
end

if all(basis(1:i)~=0)                                   % a level with an ending node
        wa=wpk(s1,h,g,order,(basis(1:i)-1),flag,0);     % but the other node continues
else wa=s1;
end

if all(basis(i+1:length(basis))~=0)
        wd=wpk(s2,h,g,order,(basis(i+1:length(basis))-1),flag,1);   % complementary case
else wd=s2;
end

if flag                                 % output is a matrix
        lwd=length(wd);
        lwx=length(wx);
        nz=2*lwd-lwx;                   % find the number of zeros
                                        % to be inserted for padding

        w=[ [s1;zeros(nz/2,1);s2;zeros(nz/2,1)] , [wa;wd] ];
else
        w=[wa;wd];
end

if nargin==4
        N=size(w,1);
        w=[ [x;zeros(N-length(x),1)] w];
elseif (nargin==5)&(fil), w=w'; end

function y=wt(x,h,g,k,del1,del2)

% WT   Discrete Wavelet Transform.
%
%      WT(X,H,G,K) calculates the wavelet transform of vector X.
%      If X is a matrix (2D), WT will calculate the one dimensional
%      wavelet transform of each row vector. The second argument H
%      is the lowpass filter and the third argument G the highpass
%      filter.
%
%      The output vector contains the coefficients of the DWT ordered
%      from the low pass residue at scale K to the coefficients
%      at the lowest scale, as the following example ilustrates:
%
%      Output vector (k=3):
%
%      [------|------|------------|------------------------]
%         |       |        |                |
%         |       |        |                `-> 1st scale coefficients
%         |       |        `-----------> 2nd scale coefficients
%         |       `--------------------> 3rd scale coefficients
%         `----------------> Low pass residue  at 3rd scale
%
%
%      If X is a matrix, the result will be another matrix with
%      the same number of rows, holding each one its respective
%      transformation.
%
%      WT (X,H,G,K,DEL1,DEL2) calculates the wavelet transform of
%      vector X, but also allows the users to change the alignment
%      of the outputs with respect to the input signal. This effect
%      is achieved by setting to DEL1 and DEL2 the delays of H and
%      G respectively. The default values of DEL1 and DEL2 are
%      calculated using the function WTCENTER.

% -----------------------------------
%    CHECK PARAMETERS AND OPTIONS
% -----------------------------------

h=h(:)';        % Arrange the filters so that they are row vectors.
g=g(:)';

if length(x)<2^k
        disp('The scale is too high. The maximum for the signal is:')
        floor(log2(length(x)))
        return
end

[liy,lix]=size(x);

if lix==1               % And arrange the input vector to a row if
        x=x';           % it's not a matrix.
        trasp=1;        % (and take note of it)
        [liy,lix]=size(x);
else
        trasp=0;
end

%--------------------------
%    DELAY CALCULATION
%--------------------------

% Calculate delays as the C.O.E. of the filters
dlp=wtcenter(h);
dhp=wtcenter(g);

if rem(dhp-dlp,2)~=0            % difference between them.
        dhp=dhp+1;              % must be even
end;

if nargin==6,                   % Other experimental filter delays
        dlp=del1;               % can be forced from the arguments
        dhp=del2;
end;

%------------------------------
%    WRAPPAROUND CALCULATION
%------------------------------
llp=length(h);                  % Length of the lowpass filter
lhp=length(g);                  % Length of the highpass filter.

L=max([lhp,llp,dlp,dhp]);       % The number of samples for the
                                % wrapparound. Thus, we should need to
                                % move along any L samples to get the
                                % output wavelet vector phase equal to
                                % original input phase.

%------------------------------
%     START THE ALGORITHM
%------------------------------

for it=1:liy,           % For every row of the input matrix...
                        % (this makes one wavelet transform
                        % for each of the rows of the input matrix)
        tm=[];
        t=x(it,:);                      % Copy the vector to transform.

        for i=1:k                       % For every scale (iteration)...
                lx=length(t);
                if rem(lx,2)~=0         % Check that the number of samples
                        t=[t,0];        % will be even (because of decimation).
                        lx=lx+1;
                end
                tp=t;                   % Build wrapparound. The input signal
                pl=length(tp);          % can be smaller than L, so it can
                while L>pl              % be necessary to repeat it several
                        tp=[tp,t];      % times
                        pl=length(tp);
                end

                t=[tp(pl-L+1:pl),t,tp(1:L)];    % Add the wrapparound.

                yl=conv(t,h);           % Then do lowpass filtering ...
                yh=conv(t,g);           % ... and highpass filtering.

                yl=yl((dlp+1+L):2:(dlp+L+lx));    % Decimate the outputs
                yh=yh((dhp+1+L):2:(dhp+L+lx));    % and leave out wrapparound

                tm=[yh,tm];             % Put the resulting wavelet step
                                        % on its place into the wavelet
                                        % vector...
                t=yl;                   % ... and set the next iteration.
        end

        y(it,:)=[t,tm];                 % Wavelet vector (1 row vector)

end                             % End of the "rows" loop.

%------------------------------
%    END OF THE ALGORITHM
%------------------------------

if trasp==1                     % If the input data was a column vector
        y=y';                   % then transpose it.
end

function d=wtcenter(x,op);

%  WTCENTER Calculates the delay of filters for alignment.
%
%           WTCENTER (X) calculates the integer aproximation
%           of delay for filter X using the method set with
%           the WTMETHOD function, for alignment operations
%           in Wavelet transforms.
%
%           For a non integer value, use the CENTER function.
%
%           See also: WTMETHOD, CENTER, WT
%

global WTCENTERMETHOD

if size(WTCENTERMETHOD)==[0,0]
        WTCENTERMETHOD=0;
end

if WTCENTERMETHOD>3 | WTCENTERMETHOD<0
        WTCENTERMETHOD=0
end

d=floor(center(x,WTCENTERMETHOD));

% (Another long function !!!)

function d=center(x,op);

%  CENTER  Delay calculation for Wavelet transform alignment.
%
%          CENTER (X, OP) calculates the delay for filter in X
%          according to the alignment method indicated in OP.
%          This delay is used by Wavelet transform functions.
%          The value of OP can be:
%              0 : First Absolute Maxima Location
%              1 : Zero delay in analysis (Full for synthesis).
%              2 : Mass center (sum(m*d)/sum(m))
%              3 : Energy center (sum(m^2 *d)/sum(m^2))
%
%          If no output argument is given, then the vector X will
%          be plotted in the current figure, and a color line will be
%          marking the result.(red: OP=0; green: OP=1; cyan: OP=2;
%          blue: OP=4)

lx=length(x);
l=1:lx;

if op==1
        d=0;
else
        if op==2
                xx=abs(x(:)');
                L=l;
        end
        if op==3
                xx=x(:)'.^2;
                L=l;
        end
        if op==0
                [mx,d]=max(abs(x));
        else

                d=sum(xx.*L)/sum(xx);
        end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%% v2struct
% v2struct Pack/Unpack Variables to/from a scalar structure.
function varargout = v2struct(varargin)

%% Description
%    v2struct has dual functionality in packing & unpacking variables into structures and
%    vice versa, according to the syntax and inputs.
%
%    Function features:
%      * Pack variables to structure with enhanced field naming
%      * Pack and update variables in existing structure
%      * Unpack variables from structure with enhanced variable naming
%      * Unpack only specific fields in a structure to variables
%      * Unpack without over writing existing variables in work space
%
%    In addition to the obvious usage, this function could by highly useful for example in
%    working with a function with multiple inputs. Packing variables before the call to
%    the function, and unpacking it in the beginning of the function will make the function
%    call shorter, more readable, and you would not have to worry about arguments order any
%    more. Moreover you could leave the function as it is and you could pass same inputs to
%    multiple functions, each of which will use its designated arguments placed in the
%    structure.
%
%% Syntax
%    Pack
%      S = v2struct
%      S = v2struct(x,y,z,...)
%      S = v2struct(fieldNames)
%      S = v2struct(A,B,C,..., fieldNames)
%      S = v2struct(x,..., nameOfStruct2Update, fieldNames)
%      v2struct
%      v2struct(x,y,z,...)
%      v2struct(fieldNames)
%      v2struct(A,B,C,..., fieldNames)
%      v2struct(x,..., nameOfStruct2Update, fieldNames)
%
%    Unpack
%      v2struct(S)
%      [a,b,c,...] = v2struct(S)
%      v2struct(S,fieldNames)
%      [a,b,c,...] = v2struct(S,fieldNames)
%
%% Inputs & Outputs
%    Pack - inputs
%      x,y,z,...           - any variable to pack. can be replaced by fieldNames below.
%      nameOfStruct2Update - optional, name of structure to update if desired.
%      fieldNames          - optional, cell array of strings, which must include a cell
%                            with the string 'fieldNames' and must be the last input.
%    Pack - outputs
%      S - the packed structure. If there is no output argument then a structure named
%          Sv2struct would be created in the caller workspace.
%
%    Unpack - inputs
%      S          - name of structure to be unpacked.
%      fieldNames - optional, cell array of strings, which must include a cell with the
%                   string 'fieldNames' and must be the last input.
%    Unpack - outputs
%      a,b,c,... - variables upacked from the structure.
%                  if there are no output arguments then variables would be created in
%                  the caller workspace with naming according to name of inputs.
%
%% Examples
%  % see 'Usage example' section below for convenient presentation of these examples.
%
%    % NOTE: whenever using filedNames cell array please note the following
%    %  1. fieldNames cell array must include a cell with the string 'fieldNames'
%    %  2. fieldNames cell array input must be the last input.
%
%  % Pack
%      x = zeros(3); x2 = ones(3); y = 'Testing123'; z = cell(2,3);
%      fieldNames1 = {'fieldNames', 'x', 'y', 'z'};
%      fieldNames2 = {'fieldNames', 'a', 'b', 'c'};
%      fieldNames3 = {'fieldNames', 'x'};
%      nameOfStruct2Update = 'S';
%
%       % The four examples below return structure S with same values however the
%       % structure's field names are defined differently in every syntax.
%      % Example 1.
%      % structure field names defined by variables names.
%       S = v2struct(x,y,z)
%      % Example 2.
%      % structure field names defined according to the cell array fieldNames.
%       % NOTE: variables with the names in fieldNames1 must exist in the caller workspace.
%       S = v2struct(fieldNames1)
%      % Example 3.
%      % same as #1. but arguments are passed explicitly
%       S = v2struct(zeros(3), 'Testing123', cell(2,3), fieldNames1)
%      % Example 4.
%      % field names defined by content of fieldNames2 while
%      % the values are set according to the passed arguments. In this case the structure
%      % S returned would be: S.a=x, S.b=y, S.c=z
%       S = v2struct(x,y,z, fieldNames2)
%
%      % Example 5.
%      % update structure S. The fields that would be updated are according to content
%      % of fieldNames3. Note that you must pass a variable with the name
%      % 'nameOfStruct2Update' placed before 'fieldNames3'. This variable should contain
%      % the name of the structure you want to update as a string. Also note that if you
%      % set an output structure name which is different than the one stated in
%      % nameOfStruct2Update a new structure would be created and the structure that was
%      % meant to be updated would not get updated.
%       S.oldField = 'field to be saved for future use'
%       S = v2struct(x2, nameOfStruct2Update, fieldNames3)
%
%      % Example 6.
%      % pack all variables in caller workspace. Call without input arguments.
%        S = v2struct
%
%      % The following examples return the same results as the examples above but the
%      % structure would be returned with the default name 'Sv2struct'. Be cautious as
%      % this might lead to overriding of arguments.
%      % Example 7.
%       v2struct(x,y,z)
%      % Example 8.
%       v2struct(fieldNames1)
%      % Example 9.
%       v2struct(zeros(3), 'Testing123', cell(2,3), fieldNames1)
%      % Example 10.
%       v2struct(x,y,z, fieldNames2)
%      % Example 11.
%       S.oldField = 'field to be saved for future use'
%       v2struct(x2, nameOfStruct2Update, fieldNames3)
%      % Example 12.
%       v2struct
%
%  % Unpack
%      clear S x x2 y z fieldNames1 fieldNames2 fieldNames3 nameOfStruct2Update
%      S.x = zeros(3); S.y = 'Testing123'; S.z = cell(2,3);
%      fieldNames3 = {'fieldNames','x','z'};
%
%      % Example 1.
%      % This example creates or overwrites variables x, y, z in the caller with the
%      % contents of the corresponding named fields.
%       v2struct(S)
%
%      % Example 2.
%      % This example assigns the contents of the fields of the scalar structure
%      % S to the variables a,b,c rather than overwriting variables in the caller. If
%      % there are fewer output variables than there are fields in S, the remaining fields
%      % are not extracted.
%       [a,b,c] = v2struct(S)
%
%      % Example 3.
%      % This example creates or overwrites variables x and z in the caller with the
%      % contents of the corresponding named fields.
%       v2struct(S, fieldNames3)
%
%      % Example 4.
%      % This example assigns the contents of the fields 'x' and 'z' defined by
%      % fieldNames3 of the scalar structure S to the variables a and b rather than
%      % overwriting variables in the caller. If there are fewer output variables than
%      % there are fields in S, the remaining fields are not extracted.
%       [a,b] = v2struct(S, fieldNames3)
%
%       % This example unpacks variables 'y' and 'z' only without overwriting variable 'x'.
%       % NOTE the addition of the field named 'avoidOverWrite' to the structure to be
%       % unpacked. This is mandatory in order to make this functionality work. The
%       % contents of this field can be anything, it does not matter.
%      S.avoidOverWrite = '';
%      x = 'do not overwrite me';
%      v2struct(S)
%
%% Usage example (includes sub-functions)
%    1. run attached v2structDemo1.m file for on screen presentation of examples.
%    2. run attached v2structDemo2.m file and read comments in file for a suggestion of
%       how to use v2struct in managing input to other functions with improved usability.
%
%% Revision history
%    2011-05-19, Adi N., Creation
%    2011-05-29, Adi N., Update structure added, some documentation and demo function changes
%    2011-06-02, Adi N., Fixed updating structure functionality
%    2011-06-05, Adi N., Added functionality: avoid overwritring existing variables, added
%                        unpacking examples to demo1 .m file.
%    2011-06-30, Adi N., fieldNames usage corrected, now must include a specific string to
%                        be triggered. Documentation enhanced. Code tweaked.
%    2011-07-14, Adi N., Fixed bug in packing with variables only
%    2011-08-14, Adi N., Clarified warning and error when packing/unpacking with
%                        fieldNames.
%    2011-09-12, Adi N., Added easy packing of all variables in caller workspace (thanks
%                        to Vesa Lehtinen for the suggestion), fixed bug in warning
%                        handling in packing case, edited comments.
%
%    Inspired by the function: mmv2truct - D.C. Hanselman, University of Maine, Orono, ME
%    04469 4/28/99, 9/29/99, renamed 10/19/99 Mastering MATLAB 5, Prentice Hall,
%    ISBN 0-13-858366-8

% parse input for field names
if isempty(varargin)
   gotCellArrayOfStrings = false;
   toUnpackRegular = false;
   toUnpackFieldNames = false;
   gotFieldNames = false;
else
   gotCellArrayOfStrings = iscellstr(varargin{end});

   toUnpackRegular = (nargin == 1) && isstruct(varargin{1});
   if toUnpackRegular
      fieldNames = fieldnames(varargin{1})';
      nFields = length(fieldNames);
   end

   gotFieldNames = gotCellArrayOfStrings & any(strcmpi(varargin{end},'fieldNames'));
   if gotFieldNames
      fieldNamesRaw = varargin{end};
      % indices of cells with actual field names, leaving out the index to 'fieldNames' cell.
      indFieldNames = ~strcmpi(fieldNamesRaw,'fieldNames');
      fieldNames = fieldNamesRaw(indFieldNames);
      nFields = length(fieldNames);
   end
   toUnpackFieldNames = (nargin == 2) && isstruct(varargin{1}) && gotFieldNames;
end

% Unpack
if toUnpackRegular || toUnpackFieldNames

   struct = varargin{1};
   assert(isequal(length(struct),1) , 'Single input nust be a scalar structure.');
   CallerWS = evalin('caller','whos'); % arguments in caller work space

   % update fieldNames according to 'avoidOverWrite' flag field.
   if isfield(struct,'avoidOverWrite')
      indFieldNames = ~ismember(fieldNames,{CallerWS(:).name,'avoidOverWrite'});
      fieldNames = fieldNames(indFieldNames);
      nFields = length(fieldNames);
   end

   if toUnpackRegular % Unpack with regular fields order
      if nargout == 0 % assign in caller
         for iField = 1:nFields
            assignin('caller',fieldNames{iField},struct.(fieldNames{iField}));
         end
      else % dump into variables
         for iField = 1:nargout
            varargout{iField} = struct.(fieldNames{iField});
         end
      end

   elseif toUnpackFieldNames % Unpack with fields according to fieldNames
      if nargout == 0 % assign in caller, by comparing fields to fieldNames
         for iField = 1:nFields
            assignin('caller',fieldNames{iField},struct.(fieldNames{iField}));
         end
      else % dump into variables
         assert( isequal(nFields, nargout) , ['Number of output arguments',...
            ' does not match number of field names in cell array']);
         for iField = 1:nFields
            varargout{iField} = struct.(fieldNames{iField});
         end
      end
   end

% Pack
else
   % build cell array of input names
   CallerWS = evalin('caller','whos');
   inputNames = cell(1,nargin);
   for iArgin = 1:nargin
      inputNames{iArgin} = inputname(iArgin);
   end
   nInputs = length(inputNames);

      % look for 'nameOfStruct2Update' variable and get the structure name
   if ~any(strcmpi(inputNames,'nameOfStruct2Update')) % no nameOfStruct2Update
      nameStructArgFound = false;
      validVarargin = varargin;
   else % nameOfStruct2Update found
      nameStructArgFound = true;
      nameStructArgLoc = strcmp(inputNames,'nameOfStruct2Update');
      nameOfStruct2Update = varargin{nameStructArgLoc};
      % valid varargin with just the inputs to pack and fieldNames if exists
      validVarargin = varargin(~strcmpi(inputNames,'nameOfStruct2Update'));
      % valid inputNames with just the inputs name to pack and fieldNames if exists
      inputNames = inputNames(~strcmpi(inputNames,'nameOfStruct2Update'));
      nInputs = length(inputNames);
      % copy structure from caller workspace to enable its updating
      if ismember(nameOfStruct2Update,{CallerWS(:).name}) % verify existance
         S = evalin('caller',nameOfStruct2Update);
      else
         error(['Bad input. Structure named ''',nameOfStruct2Update,...
            ''' was not found in workspace'])
      end
   end

   % when there is no input or the input is only variables and perhaps
   % also nameOfStruct2Update
   if ~gotFieldNames
      % no input, pack all of variables in caller workspace
      if isequal(nInputs, 0)
         for iVar = 1:length(CallerWS)
            S.(CallerWS(iVar).name) = evalin('caller',CallerWS(iVar).name);
         end
         % got input, check input names and pack
      else
         for iInput = 1:nInputs
            if gotCellArrayOfStrings % called with a cell array of strings
               errMsg = sprintf(['Bad input in cell array of strings.'...
                           '\nIf you want to pack (or unpack) using this cell array as'...
                           ' designated names'...
                           '\nof the structure''s fields, add a cell with the string'...
                           ' ''fieldNames'' to it.']);
            else
               errMsg = sprintf(['Bad input in argument no. ', int2str(iArgin),...
                                 ' - explicit argument.\n'...
                          'Explicit arguments can only be called along with a matching'...
                          '\n''fieldNames'' cell array of strings.']);
            end
            assert( ~isempty(inputNames{iInput}), errMsg);
            S.(inputNames{iInput}) = validVarargin{iInput};
         end

         % issue warning for possible wrong usage when packing with an input of cell array of
         % strings as the last input without it containing the string 'fieldNames'.
         if gotCellArrayOfStrings
            name = inputNames{end};
            % input contains structure and a cell array of strings
            if (nargin == 2) && isstruct(varargin{1})
               msgStr = [inputNames{1},''' and ''',inputNames{2},''' were'];
               % input contains any arguments with an implicit cell array of strings
            else
               msgStr = [name, ''' was'];
            end
            warnMsg = ['V2STRUCT - ''%s packed in the structure.'...
               '\nTo avoid this warning do not put ''%s'' as last v2struct input.'...
               '\nIf you want to pack (or unpack) using ''%s'' as designated names'...
               ' of the'...
               '\nstructure''s fields, add a cell with the string ''fieldNames'' to'...
               ' ''%s''.'];
            fprintf('\n')
            warning('MATLAB:V2STRUCT:cellArrayOfStringNotFieldNames',warnMsg,msgStr,...
                     name,name,name)
         end
      end
   % fieldNames cell array exists in input
   elseif gotFieldNames
      nVarToPack = length(varargin)-1-double(nameStructArgFound);
      if nVarToPack == 0 % no variables to pack
         for iField = 1:nFields
            S.(fieldNames{iField}) = evalin('caller',fieldNames{iField});
         end

         % else - variables to pack exist
         % check for correct number of fields vs. variables to pack
      elseif ~isequal(nFields,nVarToPack)
         error(['Bad input. Number of strings in fieldNames does not match',...
                'number of input arguments for packing.'])
      else
         for iField = 1:nFields
            S.(fieldNames{iField}) = validVarargin{iField};
         end
      end

   end % if ~gotFieldNames

if nargout == 0
   assignin( 'caller', 'Sv2struct',S );
else
   varargout{1} = S;
end

end % if nargin

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
[x,y,z]=peaks(7); figure(1),  mesh(x,y,z)

[xi,yi]=meshgrid(-3:0.2:3,-3:0.2:3);

z1=interp2(x,y,z,xi,yi,'nearest');

z2=interp2(x,y,z,xi,yi,'linear');

z3=interp2(x,y,z,xi,yi,'spline');

z4=interp2(x,y,z,xi,yi,'cubic');

figure(2),  subplot(2,2,1)

 mesh(xi,yi,z1)

title('nearest')

subplot(2,2,2)

mesh(xi,yi,z2)

 title('linear')

subplot(2,2,3)

 mesh(xi,yi,z3)

title('spiine')

subplot(2,2,4)

mesh(xi,yi,z4)

 title('cubic')

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [timelock, cfg] = freq2timelock(cfg, freq)

% FREQ2TIMELOCK  transform the frequency data into something
% on which the timelocked source reconstruction methods can
% perform their trick.
%
% This function also performs frequency and channel selection, using
%   cfg.frequency
%   cfg.channel
%
% After source reconstruction, you should use TIMELOCK2FREQ.

% Copyright (C) 2005, Robert Oostenveld
%
% 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$

if isfield(freq, 'fourierspctrm')
  fprintf('constructing real/imag data representation from single trial fourier representation\n');
  % select the complex amplitude at the frequency of interest
  cdim = dimnum(freq.dimord, 'chan');  % should be 2
  fdim = dimnum(freq.dimord, 'freq');  % should be 3
  for i=1:numel(cfg.frequency)
    fbin(i) = nearest(freq.freq, cfg.frequency(i));
  end
  cfg.frequency = freq.freq(fbin);
  if cdim==2 && fdim==3
    % other dimords are not supported, since they do not occur
    spctrm = dimindex(freq.fourierspctrm, fdim, {fbin});
    spctrm = permute(spctrm, [2 1 3]);
  end
  % select the desired channels in the data
  cfg.channel = ft_channelselection(cfg.channel, freq.label);
  [dum, chansel] = match_str(cfg.channel, freq.label);
  spctrm = spctrm(chansel,:);
  % concatenate the real and imaginary part
  avg = [real(spctrm) imag(spctrm)];
elseif isfield(freq, 'crsspctrm')
  fprintf('constructing real/imag data representation from csd matrix\n');
  % hmmm... I have no idea whether this is correct
  cfg.channel = ft_channelselection(cfg.channel, freq.label);
  % this subfunction also takes care of the channel selection
  [Cf, Cr, Pr, Ntrials, dum] = prepare_freq_matrices(cfg, freq);
  cfg.frequency = dum.frequency;
  cfg.channel   = dum.channel;
  if length(size(Cf))==3
    % average the cross-spectrum over trials
    Cf = squeeze(mean(Cf,1));
  end
  % reconstruct something that spans the same space as the fft of the data, hmmm...
  [u, s, v] = svd(Cf);
  spctrm = u * sqrt(s);
  % concatenate the real and imaginary part
  avg = [real(spctrm) imag(spctrm)];
else
  ft_error('cannot convert this representation of frequency domain data');
end

timelock        = [];
timelock.avg    = avg;
timelock.label  = cfg.channel;
timelock.time   = 1:size(timelock.avg,2);
if isfield(freq, 'cfg'), timelock.cfg = freq.cfg; end
timelock.dimord = 'chan_time';

if isfield(freq, 'grad')
  timelock.grad = freq.grad;
end

if isfield(freq, 'elec')
  timelock.elec = freq.elec;
end


text:
function inplane = vol2ipAllROIs(volume,inplane)
%
% inplane = vol2ipAllROIs(volume,inplane)
%
% Calls vol2ipROI with all ROIs.
%
% rmk, 1/15/99
% djh, 2/2001, replaced globals with local variables

% if ~volume.selectedROI
%   myErrorDlg('Must have a selected ROI in the Volume window before it can be transformed to the Inplane window.');
% end

for r=1:length(volume.ROIs)
  volROI = volume.ROIs(r);
  ipROI = vol2ipROI(volROI,volume,inplane);

    % if no coords exist in the inplanes, offer user the option of not xforming
    % this ROI
    if isempty(ipROI.coords)
        q = 'This ROI lies outside the inplane range. Xform empty ROI?';
        resp = questdlg(q, ipROI.name);
        if isequal(resp, 'No')
            continue
        elseif ~isequal(resp, 'Yes'); % cancel or close dialog
            fprintf('%s: user aborted. \n', mfilename);
            return
        end
    end

    % add but don't select (we'll select the last one at the end)
    inplane = addROI(inplane, ipROI, 0);
end

% select last one (so only 1 redraw at most):
inplane = selectROI(inplane, length(inplane.ROIs));

return


text:
% Create sample waypoint data for trajectory generation
% NOTE: Modify this script with your own rigid body tree and
%       trajectory reference points
%       (We recommend saving a copy of this file)
%
% Copyright 2019 The MathWorks, Inc.

%% Common parameters
% Rigid Body Tree information
load gen3
load gen3positions
eeName = 'Gripper';
numJoints = numel(gen3.homeConfiguration);
ikInitGuess = gen3.homeConfiguration;

% Maximum number of waypoints (for Simulink)
maxWaypoints = 20;

% Positions (X Y Z)
waypoints = toolPositionHome' + ...
            [0 0 0.2 ; -0.1 0.2 0.4 ; -0.2 0 0.1 ; -0.1 -0.2 0.4 ; 0 0 0.2]';

% Euler Angles (Z Y X) relative to the home orientation
orientations = [0     0    0;
                pi/8  0    0;
                0    pi/2  0;
               -pi/8  0    0;
                0     0    0]';

% Array of waypoint times
waypointTimes = 0:4:16;

% Trajectory sample time
ts = 0.2;
trajTimes = 0:ts:waypointTimes(end);

%% Additional parameters

% Boundary conditions (for polynomial trajectories)
% Velocity (cubic and quintic)
waypointVels = 0.1 *[ 0  1  0;
                     -1  0  0;
                      0 -1  0;
                      1  0  0;
                      0  1  0]';

% Acceleration (quintic only)
waypointAccels = zeros(size(waypointVels));

% Acceleration times (trapezoidal only)
waypointAccelTimes = diff(waypointTimes)/4;

author: mathworks-robotics

repo: trajectory-planning-robot-manipulators


text:
function [c,r] = interactive()
%INTERACTIVE Illustrates inputs from keyboard and mouse.
%   [C,R] = INTERACTIVE() accepts from the keyboard an image filename of
%   the form 'filename.ext' (with single or double quotes). Upon typing
%   Enter (Return), the function displays the image, and activates
%   crosshairs superimposed on the image. The mouse is used to control
%   the location of the crosshairs. The coordinates [C,R] of the center
%   of the crosshairs are captured by a left click from the mouse. Data
%   acquisition is terminated by typing Enter (Return) on the keyboard.
%   The output, [C,R], consists of two column vectors of size K-by-1,
%   containing the column and corresponding row coordinates,
%   respectively. K is the number of samples taken. This function
%   assumes that the image file is in the MATLAB path.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Get image.
disp(' '); % Blank separator line.
imageName = input('Type image file name, then Enter (Return): ');

% Read the image from disk.
f = imread(imageName);

% Display instructions.
disp(' '); % Blank separator line.
disp('Image and crosshairs will be displayed next. Use the mouse')
disp('to control the location of the crosshairs, and left-click to')
disp('sample the coordinates.')
disp(' '); % Blank separator line.

% Message to user:
input('Type Enter (Return) to start AND to stop when finished:');
% Wait for user to type Enter (Return).

% Display the image and activate the crosshairs.
figure, imshow(f)
[c,r] = ginput;
close % Close the image window.

author: dipum

repo: dipum-toolbox

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


text:
function GUI = sessionGUI_editMenu(GUI);
% Attach an Edit Menu to the session GUI figure, including
% callbacks for all editing-related operations.
%
% GUI = sessionGUI_editMenu(GUI);
%
%
% ras, 07/06.
GUI.menus.edit = uimenu('Label', 'Edit', 'Separator', 'on');

% add option to edit mrSESSION
cb = '[mrSESSION ok] = EditSession(mrSESSION,0); if ok, saveSession(1); end; ';
uimenu(GUI.menus.edit, 'Label', 'Session', 'Separator', 'off', 'Callback', cb);

% add option to edit alignment for this session
% Edit alignment(using mrRx):
uimenu(GUI.menus.edit, 'Label', 'Edit/View Alignment...', 'Separator', 'off', ...
      'Callback', 'rxAlign; ');

% Attach submenus
submenu_dataType(GUI.menus.edit);
submenu_eventRelated(GUI.menus.edit);
submenu_retinotopy(GUI.menus.edit);
submenu_segmentation(GUI.menus.edit);

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

% /---------------------------------------------------------------------/ %
function h = submenu_dataType(parent);
% attach a submenu for Study-related operations.
h = uimenu(parent, 'Label', 'Data Type', 'Separator', 'off');

cb = ['dt = GUI.settings.dataType; ' ...
      '[dataTYPES(dt) ok] = EditDataType(dataTYPES(dt)); ' ...
      'if ok, saveSession(1); end; clear dt; '];
uimenu(h, 'Label', 'Edit Data Type Parameters', 'Separator', 'off', ...
          'Callback', cb);

cb = 'removeDataType(dataTYPES(GUI.settings.dataType).name);';
uimenu(h, 'Label', 'Delete Current Data Type', 'Separator', 'on', ...
          'Callback', cb);

cb = 'duplicateDataType(INPLANE{1});';
uimenu(h, 'Label', 'Duplicate Current Data Type', ...
        'Separator', 'off', 'Callback', cb);

cb = 'groupScans(INPLANE{1});';
uimenu(h, 'Label', 'Group Scans Into Data Type', ...
        'Separator', 'off', 'Callback', cb);

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

% /---------------------------------------------------------------------/ %
function h = submenu_eventRelated(parent);
% attach a submenu for editing event-related settings.
h = uimenu(parent, 'Label', 'Event-Related', 'Separator', 'off');

uimenu(h, 'Label', 'Group Scans Into Scan Group', 'Separator', 'off', ...
          'Callback', 'er_groupScans(INPLANE{1}, guiGet(''scans''), 2); ');

uimenu(h, 'Label', 'Assign Scan Group to Selected Scan', 'Separator', 'off', ...
          'Callback', 'er_groupScans(INPLANE{1}, guiGet(''scans''), 1); ');

uimenu(h, 'Label', 'Assign .par files to Scans', 'Separator', 'off', ...
          'Callback', 'er_assignParfilesToScans(INPLANE{1}, guiGet(''scans'')); ');

uimenu(h, 'Label', 'Show parfiles/scan group', 'Separator', 'off', ...
          'Callback', 'er_displayParfiles(INPLANE{1}); ');

cb = ['params = er_getParams(INPLANE{1}); ' ...
      'params = er_editParams(params); ' ...
	  'er_setParams(INPLANE{1}, params); '];
uimenu(h, 'Label', 'Edit Event-Related Parameters (scan group)', ...
          'Separator', 'off', 'Callback', cb);

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

% /---------------------------------------------------------------------/ %
function h = submenu_retinotopy(parent);
% attach a submenu for operations related to simple mapping of
% traveling-wave analysis results to retinotopic parameters.
% (This is not the full Retinotopy Model Serge created; just parameters
% that will convert corAnal phase to estimated physical units.)
% TO DO: add functionality to the analysis menu to convert corAnal files
% into polar angle and eccentricity parameter maps.
h = uimenu(parent, 'Label', 'Retinotopy', 'Separator', 'off');

uimenu(h, 'Label', 'Set Retinotopy Params (selected scans)', ...
          'Callback', 'retinoSetParams(INPLANE{1}, guiGet(''scans'')); ');

cb = 'retinoSetParams(INPLANE{1}, [], guiGet(''scans''), ''none''); ';
uimenu(h, 'Label', 'Un-set Params (selected scans)', 'Callback', cb);

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

% /---------------------------------------------------------------------/ %
function h = submenu_segmentation(parent);
% attach a submenu for segmentation-related operations.
h = uimenu(parent, 'Label', 'Segmentation', 'Separator', 'off');

uimenu(h, 'Label', 'Install / Reinstall Segmentation (cur session)', ...
          'Separator', 'off', 'Callback', 'installSegmentation; ');

uimenu(h, 'Label', 'Segmentation Info', ...
          'Separator', 'off', 'Callback', 'segmentInfo(VOLUME{1}); ');

uimenu(h, 'Label', 'Load Segmentation into mrViewer', ...
          'Separator', 'on', 'Callback', 'sessionGUI_loadSegmentation; ');

return


text:
function [hpe] = Btuph2hpe(Btuph)
% Convert power from British thermal units per hour to electrical horsepower.
% Chad A. Greene 2012
hpe = Btuph*0.000392857;


text:
function pop_wp(host, port)
% POP_WP(HOST, PORT)
%
% Pops current waypoint in Flightgear over socket connection
% (HOST, PORT).
% Run Flightgear on HOST (e.g. '127.0.0.1' = localhost) with
% parameter "--props = PORT".

import java.net.Socket
import java.io.*

fg_socket = Socket(host, port);
out = PrintWriter(fg_socket.getOutputStream(), true);

out.println(sprintf('set /autopilot/route-manager/input @pop\r\n'));

out.close();
fg_socket.close();

end


text:
function [ id, index ] = cvx_id( x )
id = x.id_;
index = x.index_;

% 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.