function [score, m, n, c, g,totLen, N,minIdx,maxIdx,alignSeq1Ge] = JPSalignmentscore(iFile,totLen)

%epsilon = 0.0001; %To handle divide by zero
%iFile = './mafft-win/hsa-miR-286-align.fasta';
fastaDtl = fastaread(iFile);
seq1 = fastaDtl(1,1).Sequence;
seq2 = fastaDtl(2,1).Sequence;
[m, n, c, g, N,minIdx,maxIdx] = getSeqCompareDetails(seq1,seq2);
alignSeq1Ge = seq1(1,minIdx:maxIdx);

score = computeScore(m, n, c, g, totLen);
end

function [numMatchPosition, noMatchPosition, maxSeqCount, gaps, totLength,minIdx,maxIdx] = getSeqCompareDetails(seq1,seq2)
    seq1_int = nt2int(seq1);
    seq2_int = nt2int(seq2);
    
    idx = find(seq2_int < 5);
    minIdx = min(idx);
    maxIdx = max(idx);
    
    seqXOR = computeXOR(seq1_int(minIdx:maxIdx),seq2_int(minIdx:maxIdx));
    %seq1_int(minIdx:maxIdx)
    %seq2_int(minIdx:maxIdx)
    %seqXOR
    numMatchPosition = length(find(seqXOR == 0));    
    totLength = size(seqXOR,2);
    noMatchPosition = totLength - numMatchPosition;
    %matchPercentage = numMatchPosition/totLength;
    seqCount = 0;
    maxSeqCount = 0;
    for i = 1:(maxIdx - minIdx + 1)
        if isequal(seqXOR(i),0)
            seqCount = seqCount + 1;
        else            
            seqCount = 0;
        end
        if (maxSeqCount < seqCount), maxSeqCount = seqCount; end
    end
    gaps = numberOfGaps(seqXOR);
end

function [gaps] = numberOfGaps(seq)
    seqLen = size(seq,2);
    
    if seq(1)>0
        gaps = 1;
    else
        gaps = 0;
    end
    
    for i = 2:seqLen
        if (seq(i)>0) && (seq(i) ~= seq(i-1))
            gaps = gaps + 1;
        end
    end
end

function [s] = computeScore(m, n, c, g, N)
    temp1 = 1/(g+1);
    temp2 = (m * c)/N;
    temp3 = log(temp2 + 1);
    temp4 = N/(n + 1);
    temp5 = temp1 * temp4;
    s = temp5 * temp3;
end

function [returnXOR] = computeXOR(seq1, seq2)
    returnXOR = [];
    for i = 1:size(seq1,2)
        if isequal(seq1(i),seq2(i)) 
            returnXOR(i) = 0;
        else
            returnXOR(i) = 1;
        end
    end
end