#utilities
def readgenfile(filename):
    count = 0 
    with open(filename,'r') as f:
        header = ''
        gotheader = False
        genomes = []
        for line in f:
            #print(line[0])
            if gotheader:
                genomes.append({'header':header.rstrip(), 'sequence':line.rstrip()})
                gotheader = False
            elif line[0]=='>':
                header = line[1:]
                gotheader = True        
            count += 1
    print ("lines read: ", count)
    print ("sequences read: ", len(genomes))
    return genomes

def diffheaders(fname1, fname2):
    genomes1 = readgenfile(fname1)
    genomes2 = readgenfile(fname2)
    headers1 = [genome["header"] for genome in genomes1 ]
    headers2 = [genome["header"] for genome in genomes2 ]
    for header in headers1:
        if not (header in headers2):
            print( header, "from", fname1, "does not exist in", fname2)
    for header in headers2:
        if not (header in headers1):
            print( header, "from", fname2, "does not exist in", fname1)
    print("diff finished!")
    
def savegenfile(filename, genomes):
    count = 0
    with open(filename,'w') as f:
        for genome in genomes:
            f.write('>'+genome["header"])
            f.write('\n')
            f.write(genome["sequence"])
            f.write('\n')
            count += 1
    print('sequences saved:', count)

def filtergenomes(genomes, filter):
    return [genome for genome in genomes if filter in genome["header"] ]

def filtergenomesbyfile(genomes, filterfilename):
    '''
    returns filtered list of genomes with headers containing lines from filterfilename
    '''
    filteredgenomes = []
    with open(filterfilename, "r") as f:
        for line in f:
            line = line.rstrip()
            filtered = filtergenomes(genomes,line)
            for genome in filtered:
                filteredgenomes.append(genome)
            if len(filtered)==0:    
                print("Not found! ", line)            
    return filteredgenomes

def filtergensize(genomes, direction, size):
    if direction=="greater": 
        return [genome for genome in genomes if len(genome["sequence"])>size ]
    elif direction=="smaller":
        return [genome for genome in genomes if len(genome["sequence"])<size ]
    else:
        return []
    
def getgensampleparams(genomes):
    '''
    returns mean and variance of sequence samples length
    '''
    glen = len(genomes)
    if glen>0:
        total = 0
        for genome in genomes: 
            total += len(genome["sequence"])
        mean = total/glen
        total = 0
        for genome in genomes: 
            diff = len(genome["sequence"])-mean
            #print(diff*diff)
            total += diff*diff
        variance = total/glen
    else:
        mean = variance = 0.0
    return mean, variance