%% Melodic Organization and Sequential Ordering of Galant Schemata
% Implications for Eighteenth-Century European Musicianship
% Gilad Rabinovitch and Aaron Carter-Ényì
% Music Theory Online 30.1.08
% Data analysis script for MATLAB R2023a

%% Data Loading
dat=readtable('GR_ACE_2024_Data.csv','Delimiter',...
    ',','Format','%s%f%s%s%f%f%f%f%f%f'); % load CSV data into table format

[n,~]=size(dat); % assign 'n' to be the number of rows in the table...
% which is the total number of schema in the corpora

%the following removes extra characters (double quotes)
%added to the specific schema column during import (fixes an import bug)
for i=1:n % for loop where 'i' is the row out of 'n' rows
    temp=char(dat{i,4}); %specific schema text
    ind=strfind(temp,'"'); %find locations ('ind') of double quotes
    temp(ind)=[]; %remove text at 'ind' locations
    dat{i,4}={temp}; %replace with corrected text
end

%the following creates a data structure p, and arrays of generic schema,
% specific schema, and transitions between schema
gen{1}=char(dat{1,3}); %initialize cell array for generic schema type
spc{1}=char(dat{1,4}); %initialize cell array for specific schema type
temp=[gen{1} ' (' spc{1} ')']; % text of the first schema
schema{1}=temp; % initiatize the cell array 'schema' with first schema
trans=[]; % an array of schema transitions (both generic and specific labels)
p(dat{1,2}).gen{1,1}=char(dat{1,3}); %p is a data structure
p(dat{1,2}).spc{1,1}=char(dat{1,4});
p(dat{1,2}).dur(1,1)=dat{1,8}; % add duration of each schema
p(dat{1,2}).stg(1,1)=dat{1,5}; % add stage categorization
for i=2:n
    gen{i}=char(dat{i,3});
    spc{i}=char(dat{i,4});
    schema{i}=[gen{i} ' (' spc{i} ')'];
    if dat{i,2}==dat{i-1,2} % if it is part of the same piece
        p(dat{i,2}).gen{end+1}=char(dat{i,3}); 
        p(dat{i,2}).spc{end+1}=char(dat{i,4});
        p(dat{i,2}).dur(end+1)=dat{i,8};
        p(dat{i,2}).stg(end+1)=dat{i,5};
        % add adjacent schema within the same piece to transitions arrays
        if isempty(trans)==1
            trans{1,1}=schema{end-1};
            trans{1,2}=schema{end};
            transgen{1,1}=gen{end-1};
            transgen{1,2}=gen{end};
            transspc{1,1}=spc{end-1};
            transspc{1,2}=spc{end};
        else
            trans{end+1,1}=schema{end-1};
            trans{end,2}=schema{end};
            transgen{end+1,1}=gen{end-1};
            transgen{end,2}=gen{end};
            transspc{end+1,1}=spc{end-1};
            transspc{end,2}=spc{end};
        end
    else % if it is a new piece, start a new substructure p(2), p(3) ...
        p(dat{i,2}).gen{1}=char(dat{i,3});
        p(dat{i,2}).spc{1}=char(dat{i,4});
        p(dat{i,2}).dur(1)=dat{i,8};
        p(dat{i,2}).stg(1)=dat{i,5};
    end
end

%% Lexicons
uschema=unique(schema); % sorted lexicon of complete schema
ugen=unique(gen); % sorted lexicon of generic schema
uspc=unique(spc); % sorted lexicon of specific schema

% Frequencies for lexicons
[ugen,~,IC]=unique(gen);
genfreq=zeros(length(ugen),1); % frequency of generic schema
for i=1:numel(IC)
    genfreq(IC(i))=genfreq(IC(i))+1;
end

[uspc,~,IC]=unique(spc);
spcfreq=zeros(length(uspc),1); % frequency of specific schema
for i=1:numel(IC)
    spcfreq(IC(i))=spcfreq(IC(i))+1;
end

[uschema,~,IC]=unique(schema);
schemafreq=zeros(length(uschema),1); % frequency of complete schema
for i=1:numel(IC)
    schemafreq(IC(i))=schemafreq(IC(i))+1;
end

%% Confusion Matrices
cm=zeros(length(ugen)); %initialize confusion matrix (cm)
for k=1:length(transgen)
    a=transgen{k,1};
    b=transgen{k,2};
    for i=1:numel(ugen)
        if strcmp(a,ugen{i})==1
            break
        end
    end
    for j=1:numel(ugen)
        if strcmp(b,ugen{j})==1
            break
        end
    end
    cm(i,j)=cm(i,j)+1;
end

% Confusion Matrix by percent ('a' goes to 'b' out of total 'a')
for i=1:length(cm)
    linesum=sum(cm(i,:));
    for j=1:length(cm)
        cmp(i,j)=cm(i,j)/linesum;
    end
end

%% Most common transitions
[n,y]=size(cm);
a=[];
b=[];
abfreq=[];
abperc=[];
for i=1:n
    for j=1:y
        if cm(i,j)>5
            if isempty(a)==1
                a{1,1}=ugen{i};
                b{1,1}=ugen{j};
                abfreq(1,1)=cm(i,j);
                abperc(1,1)=cmp(i,j);
            else
                a{end+1,1}=ugen{i};
                b{end+1,1}=ugen{j};
                abfreq(end+1,1)=cm(i,j);
                abperc(end+1,1)=cmp(i,j);
            end
        end
    end
end

%% Openers and Closers
for i=1:numel(p)
    o{i,1}=p(i).gen{1};
    c{i,1}=p(i).gen{end};
end

[uop,~,IC]=unique(o);
opfreq=zeros(length(uop),1); % frequency of opening schema
for i=1:numel(IC)
    opfreq(IC(i))=opfreq(IC(i))+1;
end

[ucl,IA,IC]=unique(c);
clfreq=zeros(length(ucl),1); % frequency of closing schema
for i=1:numel(IC)
    clfreq(IC(i))=clfreq(IC(i))+1;
end

%% Scale Degree Transitions
ots=transspc;
for i=1:length(transspc)
    for k=1:2
        transspc{i,k}=[];
        for j=1:length(ots{i,k})
            temp=char(ots{i,k});
            if isnan(str2double(temp(j)))
                if strcmp(temp(j),'#')==1
                    transspc{i,k}(end)=transspc{i,k}(end)+.4;
                else
                    transspc{i,k}(end)=transspc{i,k}(end)-.4;
                end
            else
                transspc{i,k}(end+1)=str2double(temp(j));
            end
        end
    end
end

for i=1:length(transspc)
    transSD(i,1)=transspc{i,1}(end);
    transSD(i,2)=transspc{i,2}(1);
end
SDlist1=unique(transSD(:,1));
SDlist2=unique(transSD(:,2));
cmSD=zeros(numel(SDlist1),numel(SDlist2)); %initialization of ...
% confusion matrix of Scale Degree transitions between Schema
for i=1:length(transspc)
    for j=1:length(SDlist1)
        if transSD(i,1)==SDlist1(j)
            for k=1:length(SDlist2)
                if transSD(i,2)==SDlist2(k)
                    cmSD(j,k)=cmSD(j,k)+1;
                end
            end
        end
    end
end

% all scale-degree transitions
for i=1:numel(p)
    p(i).SD=[];
    for j=1:numel(p(i).spc)
        temp=char(p(i).spc{j});
        for k=1:length(temp)
            if isnan(str2double(temp(k)))
                if strcmp(temp(k),'#')==1
                    p(i).SD(end)=p(i).SD(end)+.4;
                else
                    p(i).SD(end)=p(i).SD(end)-.4;
                end
            else
                p(i).SD(end+1)=str2double(temp(k));
            end
        end
    end
end
allSD=[];
for i=1:numel(p)
    allSD=[allSD,p(i).SD];
end
uSD=unique(allSD);
allSDtrans=zeros(numel(uSD),numel(uSD)); % initialization of ...
% confusion matrix of all SD transitions
for i=1:numel(p)
    for j=2:numel(p(i).SD)
        for k=1:numel(uSD)
            if p(i).SD(j)==uSD(k) %consequent
                for l=1:numel(uSD)
                    if p(i).SD(j-1)==uSD(l) %antecedent
                        allSDtrans(l,k)=allSDtrans(l,k)+1;
                    end
                end
            end
        end
    end
end

% All scale-degree transitions by percent of pre-transition SD frequency
allSDtpct=allSDtrans;
for i=1:length(allSDtpct)
    allSDtpct(i,:)=allSDtpct(i,:)/sum(allSDtpct(i,:));
end

%How many pieces does a generic schema appear in?
schinv=zeros(27,21); %schema inventory
for i=1:27
    for j=1:y
        nump=numel(p(i).gen);
        for k=1:nump
            if strcmp(ugen(j),p(i).gen(k))==1
                schinv(i,j)=1;
                break
            end
        end
    end
end

%% Alterations to data

temp=schinv(:,3)+schinv(:,5); %combine cudworth (5) and/or complete cadence (3)
temp(temp==2)=1;
schinv(:,3)=temp;
schinv(:,5)=temp;

schinv(:,16)=0; %% set non-specific instances to zero

% sort by frequency
sumschinv=sum(schinv);
[B,I] = sort(sumschinv,'descend');
sumschinv=sumschinv';
sortugen=ugen(I); % list of generic schema from highest to lowest frequency

%% Create thresholded tables
for i=1:27 % i is the threshold
    pt=dat;
    for j=1:y % y is the number of unique generic schema
        if sumschinv(j)<i
            [n,~]=size(pt);
            for k=n:-1:1 % start at end of table
                if strcmp(ugen(j),pt{k,3})==1
                    pt(k,:)=[];
                end
            end
        end
    end
    writetable(pt,['GR_ACE_Thresh' num2str(i)]);
end

%% Prepare data for visualization
clear c
for i=1:27
    clear a
    b=numel(p(i).SD);
    midend=0;
    for j=1:numel(p(i).spc)
        temp(j)=p(i).dur(j);
    end
    piecedur=sum(temp);
    for j=1:numel(p(i).spc)
        a(1,j)=numel(erase(char(p(i).spc{j}),'#'))/b;
        a(1,j)=(p(i).dur(j))/piecedur;
        a(2,j)=p(i).stg(j);
    end
    c{i}=a;
end

for i=1:27
    d(i)=sum(c{i}(1,:));
    if d(i)~=1
        c{i}(1,:)=c{i}(1,:)/d(i);
    end
end

cmap=colormap(prism(6));
cmap(1,:)=[1,1,1];
%figure
hold on
for i=1:27
    for j=1:numel(c{i}(1,:))
        col=cmap(c{i}(2,j)+1,:);
        w(i,j)=c{i}(1,j)*100;
        if j>1
            n(i,j)=sum(c{i}(1,1:j-1))*100;
        else
            n(i,j)=0;
        end
        rectangle('Position',[n(i,j),i-.4,w(i,j),.8],'FaceColor',col,...
            'EdgeColor',[1 1 1]);
    end
end
axis ij
xlim([0 100])
ylim([.5 27.5])
xticks([0 100])
xticklabels({'Beginning','End'})
ticks=1:27;
yticks(ticks)
yticklabels(1:27);

%% Calculate percentage white space
whitespace=zeros(1,27);
for i=1:27
    whitespace(i)=sum(c{i}(1,c{i}(2,:)==0));
end
mws=mean(whitespace);
