function IRSerrors()
% calculates statistics of spectrum prediction errors
% run after spectrum calculi for molecular model

global InfTable 
tic
% define foldername
modelfolder = getmodel;
fullname = strcat(modelfolder,'Bstruct.mat');
load(fullname);
fullname = strcat(modelfolder,'Umiu.mat');
load(fullname);

% read model files and initial parameters
InfTable = zeros(1,9);
filename = strcat(modelfolder,'InfTable.txt');
fid = fopen(filename,'r');
InfTable = fscanf(fid,'%i',9);
fclose(fid);
mfreqs = InfTable(9);
matrb(modelfolder);

filename = strcat(modelfolder,'params_std.txt');
load(filename);

% bootstrap loop
nB = 100; % B suffix is for Bootstrap
freqsB = zeros(nB,mfreqs);
absintsB = zeros(nB,mfreqs);
for iB = 1:nB
    params_set(modelfolder,params_std)
    MatrixU
    diag2
    intir
    filename = strcat(modelfolder,'freqs.txt');
    fid = fopen(filename,'rt');
    freqsB(iB,:) = fscanf(fid,'%f',mfreqs);
    fclose(fid);
    filename = strcat(modelfolder,'IRSabsint.txt');
    load(filename);
    absintsB(iB,:) = IRSabsint(:,1)';
    clear IRSabsint
end % iB
% end bootstrap loop

IRSstatistics = zeros(4,mfreqs);
IRSstatistics(1,:) = mean(freqsB);
IRSstatistics(2,:) = std(freqsB);
IRSstatistics(3,:) = mean(absintsB);
IRSstatistics(4,:) = std(absintsB);

filename = strcat(modelfolder,'IRSstatistics.txt');
fid = fopen(filename,'w');
fprintf(fid,' %7.2f %5.2f %5.2f %5.2f \n',IRSstatistics);
fclose(fid);

filename = strcat(modelfolder,'FreqsInts.mat');
save(filename,'freqsB','absintsB');
toc

% subfunctions ---------------------
function modelfolder = getmodel()
% gets model name in LevInfinite

filename = 'C:\LevInfinite\Config\Molecules.ini';
fid = fopen(filename,'rt');
molfolder = fgetl(fid);
fclose(fid);

filename = 'C:\LevInfinite\Config\ToCalculate.ini';
fid = fopen(filename,'rt');
modelname = fgetl(fid);
fclose(fid);

modelfolder = strcat(molfolder,modelname);
modelfolder = strcat(modelfolder,'\');

% --------------------------------------------------------------------
function params_set(pname,params_std)
% Transfer UqType to UqSparse.txt

global InfTable Bonds AtomMasses
global Alphas Gammas Ros His UqType  Qtypes Bx By Bz

t = clock;
randn('seed',t(6));
% prepare new parameters random set
nTypes = length(UqType);
new_params = zeros(nTypes,4);
for i = 1:nTypes
    u = UqType(i).u;
    param_std = params_std(i,1);
    if param_std > 0,u = u + rnd_limited(param_std,1.0);end
    new_params(i,1) = u;
    miu = UqType(i).miu;
    param_std = params_std(i,2);
    if param_std > 0,miu = miu + rnd_limited(param_std,1.0);end
    new_params(i,2) = miu;
    dm1dq2 = UqType(i).dm1dq2;
    param_std = params_std(i,3);
    if param_std > 0,dm1dq2 = dm1dq2 + rnd_limited(param_std,1.0);end
    new_params(i,3) = dm1dq2;
    dm2dq1 = UqType(i).dm2dq1;
    param_std = params_std(i,4);
    if param_std > 0,dm2dq1 = dm2dq1 + rnd_limited(param_std,1.0);end
    new_params(i,4) = dm2dq1;
end % i

% new parameters output for bootstrap
fullname = strcat(pname,'Uqsparse.txt');
Ufid = fopen(fullname,'w');

nBonds = InfTable(2);
nAlphas = InfTable(3);
nGammas = InfTable(4);
nRos = InfTable(5);
n = InfTable(7);

for i = 1:n
    iu = i;
    Qtyp1 = Qtypes(iu);
    switch Qtyp1
        case 'q'
            Atoms1 = [Bonds(iu,1),Bonds(iu,2)];
        case 'a'
            ind1 = iu - nBonds;  
            Atoms1 = [Alphas(ind1,1),Alphas(ind1,2),Alphas(ind1,3)];
        case 'g'
            ind1 = iu - (nBonds+nAlphas);  
            Atoms1 = [Gammas(ind1,1),Gammas(ind1,2),Gammas(ind1,3)];
        case 'r'
            ind1 = iu - (nBonds+nAlphas+nGammas);  
            Atoms1 = [Ros(ind1,1),Ros(ind1,2),Ros(ind1,3),Ros(ind1,4)];
        case 'h'
            ind1 = iu - (nBonds+nAlphas+nGammas+nRos);  
            Atoms1 = [His(ind1,1),His(ind1,2),His(ind1,3), ...
                    His(ind1,7),His(ind1,8),His(ind1,9)];
    end %  switch Qtyp1

    for j = i:n
        ju = j;
        Qtyp2 = Qtypes(ju);
        switch Qtyp2
            case 'q'
                Atoms2 = [Bonds(ju,1),Bonds(ju,2)];
            case 'a'
                ind2 = ju - nBonds;  
                Atoms2 = [Alphas(ind2,1),Alphas(ind2,2),Alphas(ind2,3)];
            case 'g'
                ind2 = ju - (nBonds+nAlphas);  
                Atoms2 = [Gammas(ind2,1),Gammas(ind2,2),Gammas(ind2,3)];
            case 'r'
                ind2 = ju - (nBonds+nAlphas+nGammas);  
                Atoms2 = [Ros(ind2,1),Ros(ind2,2),Ros(ind2,3),Ros(ind2,4)];
            case 'h'
                ind2 = ju - (nBonds+nAlphas+nGammas+nRos);  
                Atoms2 = [His(ind2,1),His(ind2,2),His(ind2,3), ...
                        His(ind2,7),His(ind2,8),His(ind2,9)];
        end %  switch Qtyp2
        [M1,rc1] = mass_centrum(Atoms1);
        [M2,rc2] = mass_centrum(Atoms2);
        sample.QType1 = Qtyp1;
        sample.QType2 = Qtyp2;
        sample.M = M1 + M2;
        sample.R12 = norm(rc2 - rc1);
        Atoms = [Atoms1,Atoms2];
        tenz = tenzor(Atoms);
        sample.t1 = tenz(1);
        sample.t2 = tenz(2);
        sample.t3 = tenz(3);
        sample.diag = 0;
        if iu == ju
            sample.diag = 1;
        end % if
        i1 = Atoms1(2);
        i2 = Atoms2(2);
        sample.csn = [Bx(i,i1), By(i,i1), Bz(i,i1)]* ...
                     [Bx(j,i2), By(j,i2), Bz(j,i2)]';
        sample.csn = abs(sample.csn);         
        % check iType coincide
        [coincide, uTypeID] = check_type(sample);
        if coincide
            u = new_params(uTypeID,1);
            if abs(u) > 0.0001
                fprintf(Ufid,'%4i %4i %9.5f \n',iu, ju, u);
            end % if
        end % if coincide
    end % j
end % i
fclose(Ufid);

Miu = zeros(1,nBonds);
fullname = strcat(pname,'dmdqSparse.txt');
dMdqfid = fopen(fullname,'w');
Qtyp1 = 'q';
for i = 1:nBonds
    iu = i;
    Atoms1 = [Bonds(iu,1),Bonds(iu,2)];
    % is the bond symmetrical?
    a1 = Atoms1(1);
    a2 = Atoms1(2);
    symmetrical = 0;
    if AtomMasses(a1) == AtomMasses(a2),symmetrical = 1;end
    if symmetrical
        % try envirements of the two atoms
        M1 = 0.0;
        M2 = 0.0;
        for k = 1:nBonds
            if Bonds(k,1) == a1,M1 = M1 + AtomMasses(Bonds(k,2));end
            if Bonds(k,2) == a1,M1 = M1 + AtomMasses(Bonds(k,1));end
            if Bonds(k,1) == a2,M2 = M2 + AtomMasses(Bonds(k,2));end
            if Bonds(k,2) == a2,M2 = M2 + AtomMasses(Bonds(k,1));end
        end % k
        if M1 ~= M2,symmetrical = 0;end 
    end % if symmetrical
%     if symmetrical, continue, end
    for j = i:n
        ju = j;
        Qtyp2 = Qtypes(ju);
        switch Qtyp2
            case 'q'
                Atoms2 = [Bonds(ju,1),Bonds(ju,2)];
            case 'a'
                ind2 = ju - nBonds;  
                Atoms2 = [Alphas(ind2,1),Alphas(ind2,2),Alphas(ind2,3)];
            case 'g'
                ind2 = ju - (nBonds+nAlphas);  
                Atoms2 = [Gammas(ind2,1),Gammas(ind2,2),Gammas(ind2,3)];
            case 'r'
                ind2 = ju - (nBonds+nAlphas+nGammas);  
                Atoms2 = [Ros(ind2,1),Ros(ind2,2),Ros(ind2,3),Ros(ind2,4)];
            case 'h'
                ind2 = ju - (nBonds+nAlphas+nGammas+nRos);  
                Atoms2 = [His(ind2,1),His(ind2,2),His(ind2,3), ...
                        His(ind2,7),His(ind2,8),His(ind2,9)];
        end %  switch Qtyp2
        [M1,rc1] = mass_centrum(Atoms1);
        [M2,rc2] = mass_centrum(Atoms2);
        sample.QType1 = Qtyp1;
        sample.QType2 = Qtyp2;
        sample.M = M1 + M2;
        sample.R12 = norm(rc2 - rc1);
        Atoms = [Atoms1,Atoms2];
        tenz = tenzor(Atoms);
        sample.t1 = tenz(1);
        sample.t2 = tenz(2);
        sample.t3 = tenz(3);
        sample.diag = 0;
        if iu == ju
            sample.diag = 1;
        end % if
        i1 = Atoms1(2);
        i2 = Atoms2(2);
        sample.csn = [Bx(i,i1), By(i,i1), Bz(i,i1)]* ...
                     [Bx(j,i2), By(j,i2), Bz(j,i2)]';
        sample.csn = abs(sample.csn);         
        % check iType coincide
        [coincide, uTypeID] = check_type(sample);
        if coincide
            miu = new_params(uTypeID,2);
            dm1dq2 = new_params(uTypeID,3);
            dm2dq1 = new_params(uTypeID,4);
            if symmetrical
                miu = 0.0;
                dm1dq2 = 0.0;
            end % if symmetrical
            if sample.diag
                Miu(i) = miu;
                if abs(dm1dq2) > 0.0001, ...
                    fprintf(dMdqfid,'%4i %4i %9.5f \n',iu, ju, dm1dq2); end
            end % diag
            if (~sample.diag)&&(Qtyp2=='q')
                if (abs(dm1dq2) > 0.0001)&&(abs(dm2dq1) < 0.000001)
                  fprintf(dMdqfid,'%4i %4i %9.5f \n',iu, ju, dm1dq2);
                  fprintf(dMdqfid,'%4i %4i %9.5f \n',ju, iu, dm1dq2);
                end % if
                if abs(dm2dq1) > 0.0001
                  fprintf(dMdqfid,'%4i %4i %9.5f \n',ju, iu, dm2dq1); 
                end % if
            end % Qtyp2=='q'
            if (~sample.diag)&&(Qtyp2~='q')&&(abs(dm1dq2) > 0.0001)
                fprintf(dMdqfid,'%4i %4i %9.5f \n',iu, ju, dm1dq2);
            end % Qtyp2~='q'
        end % if coincide
    end % j = i:n    
end % i = 1:nBonds
fclose(dMdqfid);

fullname = strcat(pname,'Miu.txt');
Miufid = fopen(fullname,'w');
fprintf(Miufid,'%9.5f \n',Miu);
fclose(Miufid);

%-----------------------------------------
function matrb(pname)
% MatrB is for checking near interactions between vibr coords with
% common atoms

global  InfTable Bx By Bz % AtomMasses
nAtoms = InfTable(1);
nBonds = InfTable(2);
nAlphas = InfTable(3);
nGammas = InfTable(4);
nRos = InfTable(5);
nHis = InfTable(6);
n = InfTable(7);
Bx = zeros(n,nAtoms);
By = zeros(n,nAtoms);
Bz = zeros(n,nAtoms);

fullname = strcat(pname,'BqStructures.txt');
load(fullname);

% fullname = strcat(pname,'AtomMasses.txt');
% load(fullname);

ib = 0;
for i = 1:nBonds
    ib = ib + 1;
    ia = BqStructures(i,1);
    Bx(ib,ia) = BqStructures(i,3);
    By(ib,ia) = BqStructures(i,4);
    Bz(ib,ia) = BqStructures(i,5);
    ia = BqStructures(i,2);
    Bx(ib,ia) = BqStructures(i,6);
    By(ib,ia) = BqStructures(i,7);
    Bz(ib,ia) = BqStructures(i,8);
end % i
clear BqStructures
if nAlphas > 0
    fname = 'BaStructures.txt';
    fullname = strcat(pname,fname);
    load(fullname);
	for i = 1:nAlphas
        ib = ib + 1;
        ia = BaStructures(i,1);
        Bx(ib,ia) = BaStructures(i,4);
        By(ib,ia) = BaStructures(i,5);
        Bz(ib,ia) = BaStructures(i,6);
        ia = BaStructures(i,2);
        Bx(ib,ia) = BaStructures(i,7);
        By(ib,ia) = BaStructures(i,8);
        Bz(ib,ia) = BaStructures(i,9);
        ia = BaStructures(i,3);
        Bx(ib,ia) = BaStructures(i,10);
        By(ib,ia) = BaStructures(i,11);
        Bz(ib,ia) = BaStructures(i,12);
	end % i = 1:nAlphas
    clear BaStructures
end % if nAlphas > 0
if nGammas > 0
    fname = 'BgStructures.txt';
    fullname = strcat(pname,fname);
    load(fullname);

    for i = 1:nGammas
        ib = ib + 1;
        ia = BgStructures(i,1);
        Bx(ib,ia) = BgStructures(i,4);
        By(ib,ia) = BgStructures(i,5);
        Bz(ib,ia) = BgStructures(i,6);
        ia = BgStructures(i,2);
        Bx(ib,ia) = BgStructures(i,7);
        By(ib,ia) = BgStructures(i,8);
        Bz(ib,ia) = BgStructures(i,9);
        ia = BgStructures(i,3);
        Bx(ib,ia) = BgStructures(i,10);
        By(ib,ia) = BgStructures(i,11);
        Bz(ib,ia) = BgStructures(i,12);
	end % i = 1:nGammas
    clear BgStructures
end % if nGammas > 0
if nRos > 0
    fname = 'BrStructures.txt';
    fullname = strcat(pname,fname);
    load(fullname);
 
    for i = 1:nRos
        ib = ib + 1;
        ia = BrStructures(i,1);
        Bx(ib,ia) = BrStructures(i,5);
        By(ib,ia) = BrStructures(i,6);
        Bz(ib,ia) = BrStructures(i,7);
        ia = BrStructures(i,2);
        Bx(ib,ia) = BrStructures(i,8);
        By(ib,ia) = BrStructures(i,9);
        Bz(ib,ia) = BrStructures(i,10);
        ia = BrStructures(i,3);
        Bx(ib,ia) = BrStructures(i,11);
        By(ib,ia) = BrStructures(i,12);
        Bz(ib,ia) = BrStructures(i,13);
        ia = BrStructures(i,4);
        Bx(ib,ia) = BrStructures(i,14);
        By(ib,ia) = BrStructures(i,15);
        Bz(ib,ia) = BrStructures(i,16);
    end % i = 1:nRos
    clear BrStructures
end % if nRos > 0
if nHis > 0
    fname = 'BhStructures.txt';
    fullname = strcat(pname,fname);
    load(fullname);
	for i = 1:nHis
        ib = ib + 1;
        ia = BhStructures(i,1);
        Bx(ib,ia) = BhStructures(i,7);
        By(ib,ia) = BhStructures(i,8);
        Bz(ib,ia) = BhStructures(i,9);
        ia = BhStructures(i,2);
        Bx(ib,ia) = BhStructures(i,10);
        By(ib,ia) = BhStructures(i,11);
        Bz(ib,ia) = BhStructures(i,12);
        ia = BhStructures(i,3);
        Bx(ib,ia) = BhStructures(i,13);
        By(ib,ia) = BhStructures(i,14);
        Bz(ib,ia) = BhStructures(i,15);
        ia = BhStructures(i,4);
        Bx(ib,ia) = BhStructures(i,16);
        By(ib,ia) = BhStructures(i,17);
        Bz(ib,ia) = BhStructures(i,18);
        ia = BhStructures(i,5);
        Bx(ib,ia) = BhStructures(i,19);
        By(ib,ia) = BhStructures(i,20);
        Bz(ib,ia) = BhStructures(i,21);
        ia = BhStructures(i,6);
        Bx(ib,ia) = BhStructures(i,22);
        By(ib,ia) = BhStructures(i,23);
        Bz(ib,ia) = BhStructures(i,24);
    end % i = 1:nHis
    clear BhStructures
end % if nHis > 0

%-----------------------------------------
function [M,rc] = mass_centrum(Atoms)
% Find total M and rc for a vibrational coordinate

global AtomCoords AtomMasses
M = sum(AtomMasses(Atoms));
rc = [];
for i = Atoms
%    Masses = [Masses; AtomMasses(i)];
    rc = [rc; AtomMasses(i)*AtomCoords(i,:)];
end % i
rc = sum(rc)/M;

%-----------------------------------------
function tenz = tenzor(Atoms)
% Calculate inertia tenzor for Atoms in interacting vibr coordinates
% Then turn the tenzor to its main axes

global AtomCoords AtomMasses
n = length(Atoms);
Coords = AtomCoords(Atoms,:);
Masses = AtomMasses(Atoms);
M = sum(Masses);
rc = [0 0 0];
for i = 1:n
    for j = 1:3
        rc(j) = rc(j) + Masses(i)*Coords(i,j);
    end  % j
end % i
rc = rc/M;
for i = 1:n
    Coords(i,:) = Coords(i,:) - rc;
end % i
tenz0 = zeros(3);
for i = 1:n
    tenz0(1,1) = tenz0(1,1) + Masses(i)*(Coords(i,2)^2 + Coords(i,3)^2);
    tenz0(1,2) = tenz0(1,2) - Masses(i)*Coords(i,1)*Coords(i,2);
    tenz0(1,3) = tenz0(1,3) - Masses(i)*Coords(i,1)*Coords(i,3);
    tenz0(2,2) = tenz0(2,2) + Masses(i)*(Coords(i,1)^2 + Coords(i,3)^2);
    tenz0(2,3) = tenz0(2,3) - Masses(i)*Coords(i,2)*Coords(i,3);
    tenz0(3,3) = tenz0(3,3) + Masses(i)*(Coords(i,1)^2 + Coords(i,2)^2);
end % i
tenz0(2,1) = tenz0(1,2);
tenz0(3,1) = tenz0(1,3); 
tenz0(3,2) = tenz0(2,3); 
tenz = eig(tenz0);

%-----------------------------------------
function [coincide, uTypeID] = check_type(sample)
% Compare sample structure with existing UqType
global UqType 
nType = length(UqType);
coincide = 0;
uTypeID = 0;
sample_row = [sample.M,sample.R12,sample.t1,sample.t2,sample.t3, ...
              sample.diag, sample.csn];
for i = 1:nType
    if sample.QType1 ~= UqType(i).QType1 || ...
       sample.QType2 ~= UqType(i).QType2
       continue
    end % if
    UqType_row = [UqType(i).M,UqType(i).R12,UqType(i).t1,UqType(i).t2, ...
        UqType(i).t3,UqType(i).diag, UqType(i).csn]; 
    if norm(sample_row - UqType_row) > 0.01          
       continue
    else
        coincide = 1;
        uTypeID = i;
        return;
    end % if
end % i

%-----------------------------------------
function y = rnd_limited(param_std,slimit)
% Randn limited

while 1
y = param_std*randn;
    if abs(y) < param_std*slimit
        break
    end % if    
end % while
