lipidwrapper.py 137 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/python

''' Copyright (c) 2014, Jacob D. Durrant
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met: 

1. Redistributions of source code must retain the above copyright notice, this
   list of conditions and the following disclaimer. 
2. Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution. 

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

The views and conclusions contained in the software and documentation are those
of the authors and should not be interpreted as representing official policies, 
either expressed or implied, of the FreeBSD Project.'''

version = "1.14"

import sys
import math
import numpy
import scipy
from scipy.spatial import Delaunay
from scipy.spatial.distance import cdist
from scipy.spatial.distance import pdist
import getopt
import textwrap
import time
import random
import gc
import multiprocessing
import os
import shutil
import cPickle as pickle
import string
import platform
import gzip

################## MULTIPROCESSING ##################

class multi_threading():
    """A class for running calculations on multiple processors"""
    
    results = []
    
    def __init__(self, inputs, num_processors, task_class, params, progress_bar_prefix=''):
        """Launches a calculation on multiple processors
            
        Arguments:
        inputs -- A list, containing all the input required for the calculation
        num_processors -- An integer, the requested number of processors to use
        task_class -- An class, the class governing what calculations will be run on a given thread
        progress_bar_prefix -- An optional string, the prefix to append to the progress bar
        
        Returns:
        Nothing, though the objects self.results list is populated with the calculation results
        
        """
        
        # set up the progress bar
        self.results = []
        indices_to_star = []
        if len(inputs) < 50:
            indices_to_star = range(len(inputs))
        else:
            while len(indices_to_star) < 50:
                indx_to_add = random.choice(range(len(inputs)))
                if not indx_to_add in indices_to_star: indices_to_star.append(indx_to_add)

        if progress_bar_prefix != '':
            toadd = 78 - len(progress_bar_prefix) - 50
            progress_bar_prefix = progress_bar_prefix + (" "*toadd)
            sys.stdout.write(progress_bar_prefix)
        
        if num_processors == 1: # so just running on 1 processor, perhaps under windows
            single_thread = task_class()
            single_thread.total_num_tasks = len(inputs)
            single_thread.indices_to_star = indices_to_star

            single_thread.results = []
            for item in inputs: single_thread.value_func(item, None)
            
            self.results = single_thread.results
            
        else: # so it actually is running on multiple processors

            cpu_count = 1
            cpu_count = multiprocessing.cpu_count()
    
            # first, if num_processors <= 0, determine the number of processors to use programatically
            if num_processors <= 0: num_processors = cpu_count
    
            # reduce the number of processors if too many have been specified
            if len(inputs) < num_processors: num_processors = len(inputs)
            
            if len(inputs) == 0: # if there are no inputs, there's nothing to do.
                self.results = []
                return
            
            # now, divide the inputs into the appropriate number of processors
            inputs_divided = {}
            for t in range(num_processors): inputs_divided[t] = []
    
            for t in range(0, len(inputs), num_processors):
                for t2 in range(num_processors):
                    index = t + t2
                    if index < len(inputs): inputs_divided[t2].append(inputs[index])

            # now, run each division on its own processor
            running = multiprocessing.Value('i', num_processors)
            mutex = multiprocessing.Lock()
    
            arrays = []
            threads = []
            for i in range(num_processors):
                athread = task_class()
                athread.total_num_tasks = len(inputs)
                athread.indices_to_star = indices_to_star
                
                threads.append(athread)
                arrays.append(multiprocessing.Array('i',[0, 1]))
    
            results_queue = multiprocessing.Queue() # to keep track of the results
            
            processes = []
            for i in range(num_processors):
                p = multiprocessing.Process(target=threads[i].runit, args=(running, mutex, results_queue, inputs_divided[i]))
                p.start()
                processes.append(p)
    
            while running.value > 0: is_running = 0 # wait for everything to finish

            # compile all results into one list
            for thread in threads:
                chunk = results_queue.get()
                self.results.extend(chunk)
                
        if progress_bar_prefix != '': print # because the progress bar is now done

class general_task:
    """A parent class of others that governs what calculations are run on each thread"""
    
    results = []
    
    def print_star_if_appropriate(self, current_index):
        """If appropriate, prints an asterix as part of the progress bar
            
        Arguments:
        current_index -- An integer, the index of the current calculation
        
        """

        # if the current index is one of the ones that you should star, write out the asterix
        if current_index in self.indices_to_star: sys.stdout.write("*")
    
    def runit(self, running, mutex, results_queue, items):
        """Launches the calculations on this thread
            
        Arguments:
        running -- A multiprocessing.Value object
        mutex -- A multiprocessing.Lock object
        results_queue -- A multiprocessing.Queue() object for storing the calculation output
        items -- A list, the input data required for the calculation
        
        """

        for item in items: self.value_func(item, results_queue)
        
        mutex.acquire()
        running.value -= 1
        mutex.release()
        results_queue.put(self.results)

################## FUNCTIONS AND CLASSES TO EXTEND NUMPY ##################

class Quaternion:
    """A class supporting quaternion arithmetic"""
    
    def __init__(self, s, x, y, z):
        self.v = numpy.empty(4)
        self.v[0] = s
        self.v[1] = x
        self.v[2] = y
        self.v[3] = z
    
    def __str__(self):
        """String containing quaternion information in the form of s x y z
            
        Returns:
        A string, containing all information about this quaternion
            
        """
        
        return "" + str(self.v[0]) + "\t" + str(self.v[1]) + "\t" + str(self.v[2]) + "\t" + str(self.v[3])
    
    def copy_of(self):
        """Returns a copy of self"""
        
        return Quaternion(self.v[0], self.v[1], self.v[2], self.v[3])
    
    def load_from_mat(self, m):
        """Converts a rotation matrix that is pure orthogonal (det(matrix)=1) into a Quaternion. Adapted from http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/index.htm 
                
        Arguments:
        m -- A 2D numpy.array representing a pure orthogonal matrix
                
        """

        #Make sure m is a 3x3 array
        if m.shape[0] != 3 or m.shape[1] != 3:
            print "Could not load quaternion from matrix...size is not (3x3)"
            return
        
        #Check that matrix is orthogonal. m_T = m_inv
        if not numpy.array_equal(numpy.transpose(m),numpy.linalg.inv(m)):
            print "Load Quaternion error. Matrix is not orthogonal"
            return
        
        #Need to make sure that the matrix is special orthogonal
        if math.fabs(1-numpy.linalg.det(m)) > 0.000001: #Done for rounding errors
            print "Load Quaternion error.  Determinant is not 1"
            return
        
        #First calculate the sum of the diagonal elements
        t = m.trace()
        
        if t > 0:
            S = math.sqrt(t + 1.0) * 2
            self.v[0] = .25 * S
            self.v[1] = (m[2,1] - m[1,2]) / S
            self.v[2] = (m[0,2] - m[2,0]) / S
            self.v[3] = (m[1,0] - m[0,1]) / S
        elif m[0,0] > m[1,1] and m[0,0] > m[2,2]:
            S = math.sqrt(1.0 + m[0,0] - m[1,1] - m[2,2]) * 2
            self.v[0] = (m[2,1] - m[1,2]) / S
            self.v[1] = .25 * S
            self.v[2] = (m[0,1] + m[1,0]) / S
            self.v[3] = (m[0,2] + m[2,0]) / S
        elif m[1,1] > m[2,2]:
            S = math.sqrt(1.0 + m[1,1] - m[0,0] - m[2,2]) * 2
            self.v[0] = (m[0,2] - m[2,0]) / S
            self.v[1] = (m[0,1] + m[1,0]) / S
            self.v[2] = .25 * S
            self.v[3] = (m[2,1] + m[1,2]) / S
        else:
            S = math.sqrt(1.0) * 2
            self.v[0] = (m[1,0] - m[0,1]) / S
            self.v[1] = (m[0,2] + m[2,0]) / S
            self.v[2] = (m[2,1] + m[1,2]) / S
            self.v[3] = .25 * S
    
    def rep_as_44_matrix(self):
        """Creates a 4x4 matrix representation of the Quaternion.
            
        Returns:
        A 4x4 numpy array
        
        """
        
        n = self.normalize()
        qw = n.v[0]
        qx = n.v[1]
        qy = n.v[2]
        qz = n.v[3]
        
        return numpy.array([[qw,qx,qy,qz],[-qx,qw,-qz,qy],[-qy,qz,qw,-qx],[-qz,-qy,qx,qw]])

    def to_matrix(self):
        """Converts to a normalized 3x3 matrix
            
        Returns:
        A 3x3 numpy.array, corresponding to the quaternion
        
        """
        
        #First normalize
        n = self.normalize()
        qw = n.v[0]
        qx = n.v[1]
        qy = n.v[2]
        qz = n.v[3]
        return numpy.array(
                           [[1.0 - 2.0 * qy * qy - 2.0 * qz * qz, 2.0 * qx * qy - 2.0 * qz * qw, 2.0 * qx * qz + 2.0 * qy * qw],
                            [2.0 * qx * qy + 2.0 * qz * qw, 1.0 - 2.0 * qx * qx - 2.0 * qz * qz, 2.0 * qy * qz - 2.0 * qx * qw],
                            [2.0 * qx * qz - 2.0 * qy * qw,2.0 * qy * qz + 2.0 * qx * qw, 1.0 - 2.0 * qy * qy - 2.0 * qx * qx]])
                
    def add(self, q2):
        """Adds two quaternions
            
        Arguments:
        q2 -- A quaternion, to be added to self
        
        Returns:
        A Quaternion, with the values corresponding to self + q2
        
        """
        
        return Quaternion(self.v[0] + q2.v[0], self.v[1] + q2.v[1], self.v[2] + q2.v[2], self.v[3] + q2.v[3])
    
    def invert(self):
        """Takes the inverse of the quaternion for "division"
        
        Returns:
        A Quaternion, with the values corresponding to self^-1
        
        """
        
        return Quaternion(self.v[0], -1 * self.v[1], -1 * self.v[2], -1 * self.v[3])

    def minus(self, q2):
        """Multiplies two quaternions
        
        Arguments:
        q2 -- A quaternion, to be subtracted from self
        
        Returns:
        A Quaternion, with the values corresponding to self - q2
        
        """
        
        return Quaternion(self.v[0] - q2.v[0], self.v[1] - q2.v[1], self.v[2] - q2.v[2], self.v[3] - q2.v[3])
    
    def multiply(self, q2):
        """Multiplies two quaternions
            
        Arguments:
        q2 -- A quaternion, to be multiplied with self
        
        Returns:
        A Quaternion, with the values corresponding to self * q2
        
        """
        
        return Quaternion(self.v[0] * q2.v[0] - self.v[1] * q2.v[1] - self.v[2] * q2.v[2] - self.v[3] * q2.v[3],
                          self.v[1] * q2.v[0] + self.v[0] * q2.v[1] + self.v[2] * q2.v[3] - self.v[3] * q2.v[2],
                          self.v[0] * q2.v[2] - self.v[1] * q2.v[3] + self.v[2] * q2.v[0] + self.v[3] * q2.v[1],
                          self.v[0] * q2.v[3] + self.v[1] * q2.v[2] - self.v[2] * q2.v[1] + self.v[3] * q2.v[0])
    
    def normalize(self):
        """Normalizes the quaternion
            
        Returns:
        A normalized Quaternion
        
        """
        
        #First normalize
        n = math.sqrt(math.pow(self.v[0],2) + math.pow(self.v[1],2) + math.pow(self.v[2],2) + math.pow(self.v[3],2))
        
        return Quaternion(self.v[0] / n, self.v[1] / n, self.v[2] / n, self.v[3] / n)
                
    def scale(self, scalar):
        """Scales a quaternion
                
        Arguments:
        scalar -- the value to scale the quaternion by
        
        Returns:
        A Quaternion, with the values corresponding to self * scalar
            
        """
    
        return Quaternion(self.v[0] * scalar, self.v[1] * scalar, self.v[2] * scalar, self.v[3] * scalar)

def get_numpy_slice(numpy_array, indices):
    """A replacement for numpy's numpy_array[indices] functionality, where numpy_array and indices are both matrices.
    Unlike numpy's default behavior, this function returns an empty array if b is empty, rather than throwing an error.
        
    Arguments:
    numpy_array -- A numpy array
    indices -- A numpy array, the indices of the elements of numpy_array to return.
    
    Returns:
    A numpy array, numpy_array[indices] if indices is not empty, an empty numpy array otherwise.
    
    """

    try: return numpy_array[indices]
    except:
        if len(indices) == 0: return numpy.array([]) # why in the world isn't this numpy's default behavior?
        else: print "Error!"

################## MODIFYING AND MANIPULATING MOLECULAR MODELS ##################

class Molecule:
    """Loads, saves, and manupulates molecuar models."""
    
    def __init__ (self):
        """Initializes the variables of the Molecule class."""
        
        self.in_triangle_margin = True
        self.in_triangle_submargin = False
        self.headgroup_index = None
    
    def get_headgroup_index(self, lipid_headgroup_marker):
        """Get's the indices of the current molecule's headgroup
            
        Arguments:
        lipid_headgroup_marker -- A tuple of the form (chain, resname, resid, atomname) specifying the headgroup
        
        Returns:
        An integer, the index of this molecule's headgroup
        
        """

        if self.headgroup_index == None: self.headgroup_index = self.get_indices_of_mask_match(lipid_headgroup_marker)[0] # so calculate it only if it's never been calculated before
        return self.headgroup_index
    
    def load_pdb(self, filename):
        """Loads a PDB file into the current Molecule object from a file
        
        Arguments:
        filename -- a string, the name of the file to load

        """

        # Now load the file into a list
        file = open(filename,"r")
        lines = file.readlines()
        file.close()
        
        # load the molecule from the list
        self.load_pdb_from_lines(lines)
        
    def load_pdb_from_lines(self, lines):
        """Loads a PDB file into the current Molecule object from a list of PDB lines
        
        Arguments:
        lines -- a list, containing the PDB lines to load into the current object

        """

        self.__init__()

        gc.disable() # because appending objects slows down code if garbage collection turned on
        
        # set up the numpy arrays to store the data
        self.atom_inf_string_vals = numpy.empty((len(lines), 4), dtype='|S9') # chain, resname, atomname, id_keys 
        self.atom_inf_resids = numpy.empty(len(lines), dtype='|S4')
        self.all_atoms_numpy = numpy.empty((len(lines), 3))
        
        # read in the data from the lines
        count = 0
        for t in range(0,len(lines)):
            line=lines[t]
            if len(line) >= 7:
                if line[0:4]=="ATOM" or line[0:6]=="HETATM": # Load atom data (coordinates, etc.)
                    count = count + 1
                    
                    self.all_atoms_numpy[t][0] = float(line[30:38])
                    self.all_atoms_numpy[t][1] = float(line[38:46])
                    self.all_atoms_numpy[t][2] = float(line[46:54])
                    
                    resname = line[16:21].strip()
                    atomname = line[11:16].strip()

                    try: resid = line[22:26].strip()
                    except: resid = "0"

                    self.atom_inf_string_vals[t][0] = line[21:22].strip() # chain
                    self.atom_inf_string_vals[t][1] = resname # resname
                    self.atom_inf_string_vals[t][2] = atomname # atomname
                    self.atom_inf_string_vals[t][3] = resname + "_" + atomname # id_keys 
                
                    self.atom_inf_resids[t] = resid

        gc.enable()

        # now resize the array, cutting out bottom parts that were never populated
        self.atom_inf_string_vals = self.atom_inf_string_vals[:count]
        self.atom_inf_resids = self.atom_inf_resids[:count]
        self.all_atoms_numpy = self.all_atoms_numpy[:count]

    def save_pdb(self, filename):
        """Saves data to a PDB file.
        
        Arguments:
        filename -- A string, the filename to be written.
        
        """
        
        toprint = ""

        file = open(filename,"w")
        for index in range(len(self.all_atoms_numpy)): file.write(self.create_pdb_line(index) + "\n")
        file.close()
        
    def set_undo_point(self): # you can restore all atom positions to some undo point. This sets that point.
        """Sets ("saves") the undo point of all atoms. Any subsequent manipulations of atomic coordinates can be "undone" by reseting to this configuration."""
        
        self.all_atoms_numpy_undo = numpy.copy(self.all_atoms_numpy)

    def undo(self):
        """Resets the coordinates of all atoms to those saved using the set_undo_point function."""
        
        self.all_atoms_numpy = numpy.copy(self.all_atoms_numpy_undo)
        
    def rotate_mol_quat(self, rot_quat):
        """Support function that rotates a molecule according to a rotation quaternion
            
        Arguments:
        mol -- A molecule to be rotated in matrix format
        rot_quat -- Quaternion to rotate molecule
        
        """

        rot_mat = rot_quat.to_matrix()
        self.all_atoms_numpy = numpy.dot(self.all_atoms_numpy,rot_mat)
            
    def baseN(self, num, b, numerals="0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ"):
        """Return the value of a number in another base
            
        Arguments:
        num -- An integer, the number in base 10.
        b -- An integer, the number of the new base
        numerals -- An optional string, containing the numerals to use in the new base
        
        Returns:
        A string, the representation of the original integer, now in the specified base
        
        """

        return ((num == 0) and numerals[0]) or (self.baseN(num // b, b, numerals).lstrip(numerals[0]) + numerals[num % b])
    
    def create_pdb_line(self, index, output_index=None, output_resid=None):
        """Create a string formatted according to the PDB standard from the atomic information contained in this atom class.
        
        Arguments:
        index -- An integer, the index of the atom in the Molecule object.
        output_index -- An optional integer, the index to use in the PDB-line output. If not specified, index is used.
        output_resid -- An optional integer, the resid to use in the PDB-line output. If not specified, the existing resid is used.
        
        Returns:
        A string, formatted according to the PDB standard.
        
        """
        
        # use user-specified index if provided
        if output_index is None: output_index = str(index)
        else: output_index = str(output_index)
        
        # PDB format is fixed column, so if the index is too big just turn it into stars
        if len(output_index) >= 7: output_index = "******"
        
        # use the user-specified resid if provided
        if output_resid is None: output_resid = self.atom_inf_resids[index]
        else: output_resid = str(output_resid)
        
        # PDB format is fixed column, so if the resid is too big, switch over to a string identifier that is unique to each residue
        if len(output_resid) >= 5: # you need to start using letters in the resid after 9999
            # 2383280 is "a001" in base 62.
            # so map 10000 to 2383280 and convert to base 62.
            output_resid = self.baseN(int(output_resid) + 2373280, 62)
            # max using this method is 35999 residues
        
        # create the PDB line
        output = "ATOM "
        output = output + str(output_index).rjust(6) + self.atom_inf_string_vals[index][2].rjust(5) + self.atom_inf_string_vals[index][1].rjust(5) + self.atom_inf_string_vals[index][0].rjust(1) + output_resid.rjust(4)
        output = output + ("%.3f" % self.all_atoms_numpy[index][0]).rjust(12)
        output = output + ("%.3f" % self.all_atoms_numpy[index][1]).rjust(8)
        output = output + ("%.3f" % self.all_atoms_numpy[index][2]).rjust(8)
        
        return output

    def copy_of(self):
        """Create a copy of the current molecule
        
        Returns:
        A Molecule object, a copy of the current one
        
        """
        
        new = Molecule()

        new.atom_inf_string_vals=self.atom_inf_string_vals.copy()
        new.atom_inf_resids=self.atom_inf_resids.copy()
        new.all_atoms_numpy = self.all_atoms_numpy.copy()
        new.headgroup_index = self.headgroup_index
        
        return new

    def portion_of(self, list_of_indices):
        """Get a portion of the current molecule
        
        Arguments:
        list_of_indices -- A numpy array of integers, corresponding to the atoms in the current Molecule that should be in the return portion
        
        Returns:
        A Molecule object, containing only the atoms specified in list_of_indices
        
        """

        new = Molecule()
        
        new.atom_inf_string_vals=self.atom_inf_string_vals[list_of_indices]
        new.atom_inf_resids=self.atom_inf_resids[list_of_indices]
        new.all_atoms_numpy=self.all_atoms_numpy[list_of_indices]
        
        return new

    def get_indices_of_mask_match(self, masks):
        """Get the indices of atoms that match queries (masks)
        
        Arguments:
        masks -- A list of tuples. Each tuple is a mask/query of the format (chain, resname, resid, atomname). '', None, or -9999 all mean wildcard. For example: [('A', 'CHL1', 57, 'O3'), ('B', '', 783, 'P')]
        
        Returns:
        A numpy array of integers, containing the indices of the atoms that match the mask
        
        """

        indices = numpy.array([], dtype=int)
        
        for mask in masks:
            chain = mask[0]
            resname = mask[1]
            resid = mask[2]
            atomname = mask[3]
            
            # get all the indices of the ones that have the same resname
            if chain == "" or chain is None or chain == -9999: indices_of_ones_with_this_chain = numpy.array(range(len(self.atom_inf_string_vals))) # so it can be anything
            else: indices_of_ones_with_this_chain = numpy.nonzero(self.atom_inf_string_vals[:,0] == chain)[0]

            if resname == "" or resname is None or resname == -9999: indices_of_ones_with_this_resname = numpy.array(range(len(self.atom_inf_string_vals))) # so it can be anything
            else: indices_of_ones_with_this_resname = numpy.nonzero(self.atom_inf_string_vals[:,1] == resname)[0]
            
            if resid == "" or resid is None or resid == -9999 or resid == "-9999": indices_of_ones_with_this_resid = numpy.array(range(len(self.atom_inf_resids))) # so it can be anything
            else: indices_of_ones_with_this_resid = numpy.nonzero(self.atom_inf_resids == resid)[0]

            if atomname == "" or atomname is None or atomname == -9999: indices_of_ones_with_this_atomname = numpy.array(range(len(self.atom_inf_string_vals))) # so it can be anything
            else: indices_of_ones_with_this_atomname = numpy.nonzero(self.atom_inf_string_vals[:,2] == atomname)[0]
            
            # the intersection is the one that has both
            indices_in_all = numpy.intersect1d(indices_of_ones_with_this_chain, indices_of_ones_with_this_resname, assume_unique=True)
            indices_in_all = numpy.intersect1d(indices_in_all, indices_of_ones_with_this_atomname, assume_unique=True)
            indices_in_all = numpy.intersect1d(indices_in_all, indices_of_ones_with_this_resid, assume_unique=True)
            indices = numpy.union1d(indices, indices_in_all)
            
        return indices

class Triangle():
    """A class describing a triangle in three dimensions"""
    
    def __init__(self, pts):
        """Create a Triangle object.
        
        Arguments:
        pts -- A 3x3 numpy array containing the points of the triangle
        
        """

        self.points = pts

    def __getitem__(self, index):
        """Get one of the triangel points
        
        Arguments:
        index -- An integer, the index of the point
        
        Returns:
        A (1x3) numpy array, the coordinates of the requested point
        
        """

        return self.points[index]

    def center(self):
        """Get the triangle center
        
        Returns:
        A (1x3) numpy array, the coordinates of the triangle center
        
        """

        try: return self.center_pt
        except:
            self.center_pt = numpy.average(self.points, 0)
            return self.center_pt
        
    def radii(self):
        """Get the set of distances from the trangle center to each of its points (i.e., the "radii" of the triangle)
        
        Returns:
        A 1x3 numpy array, containing the distances/radii
        
        """

        try: return self.radii_lengths
        except:
            self.radii_lengths = cdist(numpy.array([self.center()]), self.points, 'euclidean')
            return self.radii_lengths

    def max_radius(self):
        """Returns a float, the maximum distance from the triangle center to any of the contituent points"""

        try: return self.radius_length
        except:
            self.radius_length = numpy.amax(self.radii())
            return self.radius_length
        
    def project_points_onto_triangle(self, pts):
        """Projects a series of points onto the triangle plane
        
        Arguments:
        pts -- An nx3 numpy array, containing the points to be projected
        
        Returns:
        An nx3 numpy array, the coordinates of the requested point now projected onto the triangle plane
        
        """

        # define the triangle plane
        AB = self.points[1]-self.points[0]
        AC = self.points[2]-self.points[0]
        normal_to_plane = numpy.cross(AB, AC)
        normal_to_plane = normal_to_plane/numpy.linalg.norm(normal_to_plane) # so it's a unit vector now
        
        # project the headgroups onto the triangle plane
        return pts - ((numpy.transpose([numpy.dot(pts - self.points[1], normal_to_plane)])) * normal_to_plane)

    def get_indices_of_points_within_triangle_boundaries(self, pts):
        """For a set of points that are coplanar with the current triangle, identify the indices of the points that are within the triangle boundaries
        
        Arguments:
        pts -- An nx3 numpy array, containing the points to be evaluated
        
        Returns:
        A numpy array, containing the indices of the points that fall within the triangle boundaries
        
        """
        
        # if the triangle is really just a line or a point, return an empty list
        if numpy.array_equal(self.points[0], self.points[1]) or numpy.array_equal(self.points[0], self.points[2]) or numpy.array_equal(self.points[1], self.points[2]): return numpy.array([])
        
        # some error has occurred previously, return an empty list
        if numpy.isnan(self.points).any(): return numpy.array([])
        
        # get bounding box
        tri_min = numpy.min(self.points,0) - numpy.array([1e-6, 1e-6, 1e-6]) # subtract a little to avoid rounding errors
        tri_max = numpy.max(self.points,0) + numpy.array([1e-6, 1e-6, 1e-6])
        
        # identify points that couldn't possibly be in the triangle because they're outside the box
        pts_not_in_triangle = (pts < tri_min).any(1)
        pts_not_in_triangle = numpy.logical_or(pts_not_in_triangle, (pts > tri_max).any(1))
        pts_potentially_in_triangle = numpy.logical_not(pts_not_in_triangle)
        
        # get the indices of the ones that could possibly be inside the triangle
        indices_of_pts_potentially_in_triangle = numpy.nonzero(pts_potentially_in_triangle)[0]
        
        # verify which ones really are in the triangle
        indices_to_keep = []
        for t in indices_of_pts_potentially_in_triangle:
            
            # calculate three vectors from the triangle verticies to the projection
            t_v1 = self.points[0] - pts[t]
            t_v2 = self.points[1] - pts[t]
            t_v3 = self.points[2] - pts[t]
            
            # get the appropriate angles
            angle1 = angle_between(t_v1, t_v2)
            angle2 = angle_between(t_v1, t_v3)
            angle3 = angle_between(t_v2, t_v3)
            
            # sometimes, if a triangle is small and the comparison point is very far away,
            # two of the vectors can end up being the same, especially after normalization.
            # Inevitably, in this case the point is not in the triangle.
            # we should account for that.
            if angle1 == "NORMALIZED VECTORS EQUAL!" or angle2 == "NORMALIZED VECTORS EQUAL!" or angle3 == "NORMALIZED VECTORS EQUAL!": continue
                
            if math.fabs(angle1 + angle2 + angle3 - 2 * math.pi) < 0.01: # it's inside the triangle
                indices_to_keep.append(t)

        return numpy.array(indices_to_keep)
    
    def new_triangle_expanded_by_margin(self, margin):
        """Return a triangle that is larger or smaller than the current one. The triangle is not simply scaled. A new triangle is carefully constructed such that each of its edges and the edges of the original triangle are a user-specified distance apart.
        
        Arguments:
        margin -- A float, the distance between the edges of the new triangle and the corresponding edges of the original triangle
        
        Returns:
        A Triangle object, the new triangle
        
        """
        
        # first, if this triangle is already a single point and the margin is negative, you can't collapse it further
        if margin < 0.0:
            if numpy.array_equal(self.points[0], self.points[1]) and numpy.array_equal(self.points[0], self.points[2]):
                return Triangle(self.points)
        
        # get the centers of each side
        side_center_1 = numpy.average(self.points[[0,1]],0)
        side_center_2 = numpy.average(self.points[[1,2]],0)
        side_center_3 = numpy.average(self.points[[0,2]],0)
        side_centers = numpy.array([side_center_1, side_center_2, side_center_3])
        
        # get the vectors from the triangle center to each side center
        center_to_center_vec = side_centers - self.center()
        old_lengths = numpy.apply_along_axis(numpy.linalg.norm, 1, center_to_center_vec)
        new_lengths = old_lengths + margin
        
        # sanity check. None of the new lengths should be negative.
        # if any of them are, just set the whole triangle to a point
        # at the old triangle's center
        if new_lengths[0] < 0.0 or new_lengths[1] < 0.0 or new_lengths[2] < 0.0: return Triangle(numpy.array([self.center(), self.center(), self.center()]))
        
        # extend (or contract) these vectors farther out
        center_to_center_vec_normalized = (center_to_center_vec.T/old_lengths).T
        center_to_center_vec_new_length = self.center() + (center_to_center_vec_normalized.T * new_lengths).T
        
        # get an additional point on the extended side of the triangle
        second_points_on_new_side = center_to_center_vec_new_length + self.points - side_centers
        
        def seg_intersect(a1,a2, b1,b2):
            """Identify (or approximate) the point where two lines in 3D space intersect
            
            Arguments:
            a1 -- A 3x1 numpy array, a point on the first line
            a2 -- A 3x1 numpy array, a second point on the first line
            b1 -- A 3x1 numpy array, a point on the second line
            b2 -- A 3x1 numpy array, a second on the second line
            
            Returns:
            A 1x3 numpy array, the point (or an approximation thereof when an exact solution is not possible) where the defined lines intersect
            
            """
        
            # first, define the lines from the provided points
            pt1 = a1
            vec1 = a2-a1
            
            pt2 = b1
            vec2 = b2-b1
            
            # now get the points on the lines that are closest to each other
            coeffs = numpy.vstack((vec2, -vec1)).T
            best_sol_all = numpy.linalg.lstsq(coeffs, pt1-pt2)
            best_sol = best_sol_all[0]
            
            if best_sol_all[1][0] == 0.0: # an exact solution because the lines intersect
                return vec1 * best_sol[1] + pt1
            else: # return the average pt of the two points that are closest to each other
                close_pt1 = vec1 * best_sol[1] + pt1
                close_pt2 = vec2 * best_sol[0] + pt2
                
                return (close_pt1 + close_pt2) * 0.5 # return the average pt

        # get the corners of the new triangle
        pt1 = seg_intersect( center_to_center_vec_new_length[0],second_points_on_new_side[0], center_to_center_vec_new_length[1],second_points_on_new_side[1])
        pt2 = seg_intersect( center_to_center_vec_new_length[0],second_points_on_new_side[0], center_to_center_vec_new_length[2],second_points_on_new_side[2])
        pt3 = seg_intersect( center_to_center_vec_new_length[2],second_points_on_new_side[2], center_to_center_vec_new_length[1],second_points_on_new_side[1])
        
        new_pts = numpy.vstack((pt1, pt2, pt3))
        
        return Triangle(new_pts)
    
    def near_other_triangle(self, other_triangle, params):
        """Determines if another triangle is near this one
        
        Arguments:
        other_triangle -- A Triangle object, the other triangle to be evaluated
        params -- A dictionary, the user-specified command-line parameters
        
        Returns:
        A boolean, True of the two triangles are near each other, False otherwise
        
        """
        
        # check if the two triangles share a corner
        dists_between_triangle_pts = cdist(self.points, other_triangle.points)
        if True in (dists_between_triangle_pts == 0.0): return True # so they are adjacent
        
        # check if the distance to their centers is too close as well
        dist_between_center_points = numpy.linalg.norm(self.center() - other_triangle.center())
        if dist_between_center_points < params['triangle_center_proximity_cutoff_distance']: return True

############### GET COMMANDLINE PARAMETERS ###############

def get_commandline_parameters(argv):
    """Get the user-defined command-line parameters
    
    Returns:
    A dictionary, the user-specified command-line parameters
    
    """

    # first check if the user has requested the help file
    if '--help' in [t.lower() for t in argv]:
        print
        print "The initial lipid model"
        print "======================="
        print
        print textwrap.fill("--lipid_pdb_filename: This parameter specifies a PDB file containing an all-atom model of a planar lipid bilayer. LipidWrapper will wrap this lipid around the user-generated mesh. Example: --lipid_pdb_filename lipid.pdb", 70, subsequent_indent = "      ")
        print textwrap.fill("--lipid_headgroup_marker: A unique atom representing the headgroup of each lipid residue must be specified. The --lipid_headgroup_marker accepts a comma-separated lists of atom specifications (RESNAME_ATOMNAME). If either RESNAME or ATOMNAME is omitted, any value will be accepted. By default, LipidWrapper identifies lipid headgroups by looking for any atom named \"P\" (_P) or any atom named \"O3\" belonging to a cholesterol molecule (CHL1_O3). Example: --lipid_headgroup_marker \"_P,CHL1_O3\"", 70, subsequent_indent = "      ")
        print
        print "Methods for creating a surface mesh"
        print "==================================="
        print
        print textwrap.fill("--surface_equation: Generate a surface mesh from a python-formatted equation defining z, given x and y. The --min_x, --max_x, --min_y, and --max_y parameters are used to specify the region over which the function should be evaluated. The --step_x and --step_y parameters define the x-y distance between adjacent points. Python functions from the math, numpy, and scipy modules can be used. Example: --surface_equation \"z = 250*numpy.sin(x*x/60000 +y*y/60000)\"", 70, subsequent_indent = "      ")
        print textwrap.fill("--surface_filename: If this parameter specifies a file with the PDB extension, a surface mesh is generated from the coordinates of the PDB atoms. Example: --surface_filename mymesh.pdb", 70, subsequent_indent = "      ")
        print textwrap.fill("--surface_filename: If this parameter specifies a file with the DAE extension, the mesh points and triangulations will be taken from the file. Example: --surface_filename mymodel.dae", 70, subsequent_indent = "      ")
        print textwrap.fill("--surface_filename: If this parameter specifies a file that does not have the PDB extension, the file is assumed to be a gray-scale image, where black represents regions that are topologically low, and white represents regions that are topologically high. The --min_x, --max_x, --min_y, and --max_y parameters are used to specify the region where the mesh should be generated. The --step_x and --step_y parameters define the x-y distance between adjacent points. The --max_height parameter determines the height of the bilayer model at those locations where the image is white; black regions are assigned a height of 0. This feature is only available if the python PIL module has been installed on your system. Example: --surface_filename mymesh.png", 70, subsequent_indent = "      ")
        print
        print "Methods for resolving lipid clashes"
        print "==================================="
        print
        print textwrap.fill("--delete_clashing_lipids: It's common for lipids to sterically clash at the interface of two adjacent surface-mesh tessellated triangles. If this parameter is set to TRUE, any clashing lipids are deleted. Example: --delete_clashing_lipids TRUE", 70, subsequent_indent = "      ")
        print textwrap.fill("--clash_cutoff: If you do choose to delete clashing lipids, this parameter determines how close two atoms must be (in Angstroms) to constitute a steric clash. Example: --clash_cutoff 2.0", 70, subsequent_indent = "      ")
        print textwrap.fill("--clashing_potential_margin: Lipid clashes occur at the edges of adjacent tessellated triangles. If these triangles are very large, it's faster to only check for clashes and holes near the triangle edges. This variable specifies how far from the edges, in Angstroms, that LipidWrapper should look for clashes and holes. Example: --clashing_potential_margin 25.0", 70, subsequent_indent = "      ")
        print textwrap.fill("--fill_holes: Deleting lipids often leaves holes in the membrane. If this parameter is set to TRUE, LipidWrapper tries to fill the hole. Example: --fill_holes TRUE", 70, subsequent_indent = "      ")
        print textwrap.fill("--very_distant_lipids_cutoff: LipidWrapper determines if two lipids clash by comparing the distance between every atom in the first lipid with every atom in the second lipid. This can be computationally expensive. However, sometimes two lipids are so distant from each other, that it's obvious there are no clashes, making the pair-wise comparison unnecessary. Before performing this expensive pair-wise comparison, LipidWrapper calculates the distance between one atom of each lipid. If this distance is greater than this user-specified cutoff, the program will simply assume there are no clashes. WARNING: Remember to consider the width of your lipid bilayer when choosing this value. Adjacent lipids on opposite sides of the bilayer can seem distant when considering the distance between their headgroups, for example. Example: --very_distant_lipids_cutoff 50.0", 70, subsequent_indent = "      ")
        print textwrap.fill("--triangle_center_proximity_cutoff_distance: Lipid steric clashes/holes typically occur between lipids that belong to adjacent tessellated triangles. However, if tessellated triangles are small enough, clashes are possible between lipids that belong to non-adjacent triangles as well. Consequently, in addition to checking for adjacency, LipidWrapper also checks the distance between the triangle centers, using this user-specified value as a cutoff. Example: --triangle_center_proximity_cutoff_distance 50.0", 70, subsequent_indent = "      ")
        print textwrap.fill("--fill_hole_exhaustiveness: Essentially, how long LipidWrapper should try to fill the holes. Example: --fill_hole_exhaustiveness 10", 70, subsequent_indent = "      ")
        print textwrap.fill("--memory_optimization_factor: When the tessellated triangles are very large and consequently contain many individual lipids, the extensive pairwise distance comparisons required can result in memory errors. This parameter tells lipid Wrapper to divide the list of atoms being compared into smaller chunks. The pairwise distance comparison is performed piecewise on each chunk-chunk pair and so uses less memory, albeit at the expensive of speed. Only increase the value of this parameter if you run into memory errors. Example: --memory_optimization_factor 1", 70, subsequent_indent = "      ")
        print
        print "Additional options"
        print "=================="
        print
        print textwrap.fill("--number_of_processors: Using multiple processors can significantly increase the speed of the LipidWrapper algorithm. Example: --number_of_processors 8", 70, subsequent_indent = "      ")
        print textwrap.fill("--show_grid_points: Aside from producing PDB coordinates for lipid atoms, additional coordinates will be appended to the bottom of the output containing \"atoms\" named \"X\" that specify the location of the surface mesh points. Example: --show_grid_points TRUE", 70, subsequent_indent = "      ")
        print textwrap.fill("--create_triangle_tcl_file: A separate file named \"triangles.tcl\" will be generated containing a tcl script that can be run in VMD to visualize the mesh surface. Example: --create_triangle_tcl_file TRUE", 70, subsequent_indent = "      ")
        print textwrap.fill("--output_directory: If an output directory is specified, all LipidWrapper output files, as well as additional files representing the intermediate steps required to build the final bilayer, will be saved in that directory. Example: --output_directory ./my_output/", 70, subsequent_indent = "      ")
        print textwrap.fill("--use_disk_instead_of_memory: For very large systems, storing the growing model in memory can be problematic. If this parameter is set to TRUE, the growing model will be stored on the hard disk instead. However, expect longer execution times if this parameter is set to TRUE. Example: --use_disk_instead_of_memory TRUE", 70, subsequent_indent = "      ")
        print textwrap.fill("--compress_output: Depending on the user options selected, LipidWrapper output can require a lot of disk space. If this parameter is set to TRUE, the output will be automatically compressed using the gzip algorithm (Lempel-Ziv coding LZ77). The files can be uncompressed with the UNIX gunzip utility, or similar Windows-based packages. Example: --compress_output TRUE", 70, subsequent_indent = "      ")
        print
        print "Example"
        print "======="
        print
        print textwrap.fill('python lipidwrapper.py --surface_equation "z = 250*numpy.sin(x*x/60000 +y*y/60000) * (-numpy.sqrt(x*x+y*y)/(560 * numpy.sqrt(2)) + 1)" --min_x 500 --max_x 1000 --min_y 500 --max_y 1000 --step_x 25 --step_y 25 --lipid_pdb_filename lipid.pdb --lipid_headgroup_marker "_P,CHL1_O3" --delete_clashing_lipids TRUE --clash_cutoff 1.0 --fill_holes TRUE --fill_hole_exhaustiveness 10 > lipid_model.pdb', 70, subsequent_indent = "      ")
        print
        sys.exit(0)
    
    # defaults
    params = {}
    params['surface_filename'] = '' # could be a PDB or image file, depending on surface_source value
    params['surface_equation'] = 'z = 100*numpy.sin(x*x/60000 +y*y/60000) * (-numpy.sqrt(x*x+y*y)/(560 * numpy.sqrt(2)) + 1)' # used if surface_source is set to "EQUATION"
    params['min_x'] = 500 # used if surface_source is PNG or EQUATION
    params['max_x'] = 750 # used if surface_source is PNG or EQUATION
    params['min_y'] = 500 # used if surface_source is PNG or EQUATION
    params['max_y'] = 750 # used if surface_source is PNG or EQUATION
    params['step_x'] = 30 # used if surface_source is PNG or EQUATION
    params['step_y'] = 30 # used if surface_source is PNG or EQUATION
    params['max_height'] = 0 # used if surface_source is PNG
    params['lipid_pdb_filename'] = '' # the filename containing the small, planar lipid model
    params['lipid_headgroup_marker'] = '_P,CHL1_O3' # by default, any phosphate atom is considered a marker for the lipid headgroup, and also any O3 atom belonging to a cholesterol
    params['show_grid_points'] = 'FALSE'
    params['create_triangle_tcl_file'] = 'FALSE'
    params['delete_clashing_lipids'] = 'FALSE'
    params['use_disk_instead_of_memory'] = 'FALSE'
    params['clash_cutoff'] = 2.0
    params['fill_holes'] = 'TRUE'
    params['output_directory'] = ''
    params['fill_hole_exhaustiveness'] = 10
    params['number_of_processors'] = 1
    params['clashing_potential_margin'] = 25.0
    params['triangle_center_proximity_cutoff_distance'] = 50.0
    params['memory_optimization_factor'] = 1
    params['very_distant_lipids_cutoff'] = 50.0
    params['compress_output'] = "FALSE"    

    # get commandline parameters
    options, remainder = getopt.getopt(argv[1:], '', [ 'surface_filename=', 'surface_equation=', 'min_x=', 'max_x=', 'min_y=', 'max_y=', 'step_x=', 'step_y=', 'max_height=', 'lipid_pdb_filename=', 'lipid_headgroup_marker=', 'show_grid_points=', 'create_triangle_tcl_file=', 'delete_clashing_lipids=', 'clash_cutoff=', 'fill_holes=', 'fill_hole_exhaustiveness=', 'output_directory=', 'number_of_processors=', 'use_disk_instead_of_memory=', 'clashing_potential_margin=', 'triangle_center_proximity_cutoff_distance=', 'memory_optimization_factor=', 'very_distant_lipids_cutoff=', 'compress_output='])
    
    # set parameters to variables
    params_string = ['compress_output', 'surface_filename', 'surface_equation', 'lipid_pdb_filename', 'lipid_headgroup_marker', 'show_grid_points', 'create_triangle_tcl_file', 'delete_clashing_lipids', 'fill_holes', 'output_directory', 'use_disk_instead_of_memory']
    params_floats = ['very_distant_lipids_cutoff', 'memory_optimization_factor', 'triangle_center_proximity_cutoff_distance', 'clashing_potential_margin', 'min_x', 'max_x', 'min_y', 'max_y', 'step_x', 'step_y', 'max_height', 'clash_cutoff', 'fill_hole_exhaustiveness', 'number_of_processors']
    
    for opt, arg in options:
        opt = opt.replace('-','')
        if opt in params_floats: arg = float(arg)
        params[opt] = arg
    
    # some parameters should be integers
    params['fill_hole_exhaustiveness'] = int(params['fill_hole_exhaustiveness'])
    params['number_of_processors'] = int(params['number_of_processors'])
    params['memory_optimization_factor'] = int(params['memory_optimization_factor'])
    
    # directories should end in / or \ (depending on os)
    if params['output_directory'] != '' and params['output_directory'][-1:] != "/": params['output_directory'] = params['output_directory'] + "/"
    
    # check if running windows. If so, you can only use one processor
    if platform.system().lower() == "windows" and params['number_of_processors'] > 1:
        print "REMARK WARNING: Use of multiple processors is only supported on Linux and OS X."
        params['number_of_processors'] = 1

    # Print out header
    toprint = []
For faster browsing, not all history is shown. View entire blame