forked from jcforbes/gidget
-
Notifications
You must be signed in to change notification settings - Fork 0
/
exper.py
3160 lines (2821 loc) · 179 KB
/
exper.py
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
import copy
import os, argparse
import sys,glob
import shutil
import subprocess
import time
import math
import numpy as np
import pdb
import random
from bolshoireader import bolshoireader
# This is a script to allow you to run the GIDGET code with
# a wide variety of systematically varied parameters. The functions here are described
# in some detail, but if you skip to the end of the script, there are plenty of
# examples of what you can run from the command line, e.g.
# $ python exper.py ex1
# will execute the example used in the README file.
t0=time.time()
allModels={} # a global dictionary of all experiments created.
allProcs=[] # a global list of all processes we've started
allStartTimes=[] # a global list of the start times of all of these processes
globalBolshoiReader = bolshoireader('rf_registry4.txt',3.0e11,3.0e14, os.environ['GIDGETDIR']+'/bolshoi/')
bolshoiSize = len(globalBolshoiReader.keys)
def HowManyStillRunning(procs):
''' Given a list of processes created with subprocess.Popen, (each of which
in this case will be a single run of the gidget code), count up how
many of those runs are still going.'''
nStillRunning=0
for proc in procs:
if(proc.poll()==None):
nStillRunning+=1
return nStillRunning
def successCode(filename):
''' Given a file assumed to contain the standard error output of
a run of the gidget code, return 0 for success (file is empty)
1 for a time step below floor result, 2 for a problem with Qst,
3 for a rejected accr. history, or 4 for a miscellaneous error code. '''
if (not os.path.exists(filename)):
return -1 # if the file *stde.txt doesn't exist, this particular run has not been started.
if os.stat(filename).st_size == 0:
return 0 # the *_stde.txt file is empty, i.e. there were no errors (so far) for this run
with open(filename,'r') as f:
contents = f.read()
if 'floor' in contents:
return 1 # the time step has gotten very small. Typically
# this is caused by 1 cell in the gas column density
# distribution becoming very small.
if 'Qst' in contents:
return 2 # a spike in the stellar column density distribution
# which in turn causes an interpolation error
if 'merger' in contents:
return 3 # accretion history rejected because of a merger
return 4 # none of the most common types of failures
class experiment:
''' This class is a container to generate and store the necessary information to run
a series of GIDGET models with parameters varied in a systematic way.'''
def __init__(self,name):
''' All we need here is a name. This will be the directory containing and the prefix for
all files created in all the GIDGET runs produced as part of this experiment.'''
# fiducial model
self.p=[name,200,1.5,.01,4.0, \
1,1,.01,1,10, \
220.0,20.0,7000.0,2.5,.5, \
1.0,2.0,5000.0,int(1e9),1.0e-3, \
1.0,-2.0/3.0,0,.5,2.0, \
2.0,0,0.0,1.5,1, \
2.0,1.0,1.0e12,5.0,3.0, \
0,2.0,-1.0,2.5,0.0, \
0.54,0.1,200,.30959,0.38, \
-0.25,1.0,1.0,0.3,1.0, \
0,0.0,.1,.03,2.0, \
.002,1.0,1.0,0.3,0.0,.054,0.0,0.0, -1.0/3.0, 10.0, 2.0, 1.0e12, \
0.1, 2.0, 2.0, 1.4, 0.5, 1.0, -1.0, 1.0, 0.5]
self.p_orig=self.p[:] # store a copy of p, possibly necessary later on.
self.pl=[self.p[:]] # define a 1-element list containing a copy of p.
# store some keys and the position to which they correspond in the p array
# When adding new variables, put them /before/ bolshoiWeight, leaving bolshoiWeight as the final
# element of this list!
self.names=['name','nx','eta','epsff','tauHeat', \
'analyticQ','cosmologyOn','xmin','NActive','NPassive', \
'vphiR','R','gasTemp','Qlim','fg0', \
'phi0','zstart','tmax','stepmax','TOL', \
'muNorm','muColScaling','b','innerPowerLaw','softening', \
'diskScaleLength','whichAccretionHistory','alphaMRI','thickness','migratePassive', \
'fixedQ','kappaMetals','Mh0','minSigSt','NChanges', \
'dbg','accScaleLength','zquench','zrelax','xiREC', \
'RfREC','deltaOmega','Noutputs','accNorm','accAlphaZ', \
'accAlphaMh','accCeiling','fscatter','invMassRatio','fcool', \
'whichAccretionProfile','alphaAccretionProfile','widthAccretionProfile','fH2Min','tDepH2SC', \
'ZIGM','fg0mult','ZIGMfac','chiZslope','deltaBeta','yREC','concentrationRandomFactor','muFgScaling', 'muMhScaling', 'ksuppress', 'kpower', 'MQuench', \
'epsquench', 'muQuench', 'stScaleReduction', 'gaScaleReduction', 'ZMix', 'energyInjectionFactor', 'CloudHeatingRate', 'AccretionHeatingRate', 'bolshoiWeight']
assert len(self.p)==len(self.names)
self.keys={}
ctr=0
for n in self.names:
self.keys[n]=ctr
ctr=ctr+1
self.expName=name
self.covariables=np.zeros(len(self.p),int)
# store the location of various expected subdirectories in the gidget distribution.
#self.base=os.getcwd() # Assume we are in the base directory - alter this to /path/to/gidget/directory if necessary
self.base = os.environ['GIDGETDIR']
self.src=self.base+'/src'
self.analysis=self.base+'/analysis'
self.bin=self.base+'/bin'
self.out=self.base+'/output'
self.xgrid=self.base+'/xgrid'
allModels[name] = self
def changeName(self,newName):
''' If this object is copied to form the starting point for a separate experiment,
the self.expName, the names of each GIDGET run, and the entry in allModels will
need to be updated. This method takes care of all of these tasks.'''
oldName = self.expName
self.p[0] = newName
self.p_orig[0] = newName
for p in self.pl:
p[0] = newName + p[0][len(oldName):] # new prefix + old suffixes
self.expName = newName
allModels[newName] = self
def vary(self,name,mmin,mmax,n,log,cov=0):
'''Set up to vary a particular parameter, name, over the
range [mmin,mmax] with n steps spaced linearly if log=0 or
logarithmically if log=1. The cov flag means that all variables
with the same value of cov, e.g. accScaleLength and whichAccretionHistory,
will be varied together rather than independently. '''
if(n>1 and log==0):
self.p[self.keys[name]]=[mmin + (mmax-mmin)*type(self.p_orig[self.keys[name]])(i)/type(self.p_orig[self.keys[name]])(n-1) for i in range(n)]
elif(n>1 and log==1 and mmin!=0 and mmin!=0.0):
rat = (float(mmax)/float(mmin))
typ = type(self.p_orig[self.keys[name]])
self.p[self.keys[name]]=[typ(float(mmin) * (rat)**(float(i)/float(n-1))) for i in range(n)]
elif(n==1):
self.p[self.keys[name]]=mmin # a mechanism for changing the parameter from its default value.
self.covariables[self.keys[name]] = cov
return self.p[self.keys[name]]
def multiply(self,name,multiple):
keynum = self.keys[name]
if(type(self.p[keynum]) == type([])): # we're dealing with a list
self.p[keynum] = [a*multiple for a in self.p[keynum]]
else:
self.p[keynum] = [self.p[keynum] * multiple]
return self.p[keynum]
def irregularVary(self,name,values,cov=0):
''' Similar to vary, except, rather than generating the
array of parameter values to use for parameter "name",
this list is specified explicitly in "values"'''
keyNum = self.keys[name] # the element of the list self.p which we're editing today.
typ = type(self.p_orig[keyNum]) # what type of variable is this in a single run of gidget?
if(type(values) == type([1]) or type(values)==type(np.array([1.0,2.0]))): # if values is a list..
# replace the given index of p with the list given by values.
self.p[keyNum]=[typ(value) for value in values]
else:
# even though we weren't given a list, put our new value in a single-parameter list.
try:
self.p[keyNum]=[typ(values)]
except:
pdb.set_trace()
# Next, we want to pass along the information we've been given on whether this parameter
# should be covaried:
self.covariables[self.keys[name]] = cov
# and finally, give back what we've done.
return self.p[keyNum]
def ConsistencyCheck(self):
''' Here we check whether all parameters in a covarying set have
the same length. Otherwise it makes little sense to vary them together.
'''
consistent = True
distinctCovSetIndices = set(self.covariables)
for ind in distinctCovSetIndices:
covIndicesForThisSet = np.where(self.covariables == ind)
lengths=[]
for var in covIndicesForThisSet:
lengths.append(len(self.p[var]))
if(len(set(lengths)) != 1):
consistent=False
print( "Problem found in the ", ind, " set of covarying variables." )
print( "This corresponds to the variables ", covIndicesForThisSet )
print( "A set of variables you have asked to covary do not have the same lengths!" )
def generatePl(self):
'''vary() will change a certain element of self.p into a list
of values for that variable. This method assumes that some
(perhaps many) variables have been varied in this manner, and
expands p into a list of well-defined parameter lists, each
one of which corresponds to a run of the program. This list can
then be sent to the xgrid or run on your machine. See other
methods for details on how to do that.'''
self.pl=[self.p_orig[:]] # in case generatePl() has already been run, reset pl.
# for each separate parameter (e.g. number of cells, dissipation rate, etc.)
for j in range(len(self.p)):
param = self.p[j]
# if instead of a single value of the parameter, we have a whole list, create a set of runs
# such that for every previous run we had, we now have len(param) more, each with a different
# value from the list param.
if(type(param)==list):
covaryOtherVarsWithJ=False
varyThisVariable=True # i.e. the variable corresponding to j.
cov = self.covariables[j] # the flag for covarying variables.
if (cov!=0): # if this parameter (j) is part of a set of covarying parameters...
covIndices = np.where(self.covariables == cov)[0]
# if our current variable, j, is the first in a set of variables we're covarying,
# let's go ahead and vary other the other variables in the same set.
if( covIndices[0] == j ):
covaryOtherVarsWithJ = True
# whereas if this is not the case, then this variable has already been varied
else:
varyThisVariable = False
if(varyThisVariable):
# make len(param) copies of the current pl
pl2=[]
for i in range(len(param)):
f=copy.deepcopy(self.pl)
pl2.append(f) # pl2 ~ [pl,pl,pl]~[[p,p,p],[p,p,p],[p,p,p]]
# in each copy, set the jth parameter in p to param[i]
for a_p in pl2[i]: # each element is a list of parameters for
a_p[j]=param[i]
if(covaryOtherVarsWithJ):
for covIndex in covIndices:
if(covIndex != j): # already taken care of with a_p[j]=...
#print covIndex, i, len(a_p), len(self.p), len(self.p[covIndex])
a_p[covIndex] = self.p[covIndex][i]
# in each copy, append to the name a...z corresponding to
# which copy is currently being edited
if(i<=25):
base='a'
app=''
# If there are more than 26 variations, add a numeral for each
# time we have looped through the alphabet when this number is
# greater than zero.
else:
base='a'
app=str((i-(i%26))/26)
# avoid adding more letters to distinguish individual models
# if such a distinction is unnecessary.
if(len(param) > 1):
# but if such a distinction is necessary, by all means:
a_p[self.keys['name']]+=chr(ord(base)+(i%26))+app
# end loop over p's in pl2[i]
# end loop over range of this varied parameter
# collapse pl to be a 1d array of p's.
self.pl=[element for sub in pl2 for element in sub]
# end of param-being-a-list contingency
else : #just a regular non-varying parameter
# If the user has used vary() but with N=1, the new argument will
# not be a list! In this case, we set this parameter in each element
# of pl to be the value the user picked, instead of the default value.
for i in range(len(self.pl)):
if(self.pl[i][j]!=param):
self.pl[i][j]=param
else:
pass # nothing to do- just use default value
# end of non-varying parameter contingency
# end of for loop over each parameter.
def write(self,name):
'''Write a text file which can be handled by GridStuffer to
run the experiment. This consists of a list of lists, each
one being a set of command line arguments to give to the
executable. Note that before running this, you need to run
generatePl() to generate this list of lists, and before that,
you should run vary() to set up which parameter(s) to change
in each run. Once you have written this file, to use it on
your local xgrid, you use GridStuffer to create a new metajob
with the "Input File" as the text file, and the "Output Folder"
as gidget/output'''
self.generatePl()
with open(name,'w') as f:
for a_p in self.pl:
f.write('xgrid -job submit -in '+self.bin+' -out '+self.analysis+'/'+self.expName+' ./gidget')
for param in a_p:
f.write(' '+repr(param))
if(self.pl[-1]!=a_p):
f.write('\n') # no return after the last line
if(os.path.exists(self.analysis+self.expName)):
print("Warning: this experiment already exists! It would be a good idea to rename your experiment or delete the pre-existing one manually.")
os.rename(name,self.xgrid+'/'+name) # move the text file to gidget/xgrid
os.mkdir(self.xgrid+'/output/'+self.expName) # make gidget/xgrid/output/name to store stde and stdo files.
print ("Prepared file ",name," for submission to xGrid.")
def ExamineExperiment(self):
self.generatePl()
varied=[] # which elements of p are lists
Nvaried=[] # how many variations for each such list.
theLists=[] # a copy of each list.
for j in range(len(self.p)):
param = self.p[j]
if(type(param)==list ):
varied.append(j)
Nvaried.append(len(param))
theLists.append(param[:])
successTable = np.ndarray(shape=tuple(Nvaried[:]),dtype=np.int64)
successTableKeys = np.zeros(tuple(Nvaried[:]))
for a_p in self.pl: # for each run..
successTableKey =[]
for i in range(len(varied)): # for each parameter that was varied
j = varied[i]
successTableKey.append(theLists[i].index(a_p[j]))
successTable[tuple(successTableKey)] = successCode(self.analysis+'/'+self.expName+'/'+a_p[0]+"_stde.txt")
successTableKeys[tuple(successTableKey)] = 1#tuple(successTableKey)
print()
print( "Success Table:")
print (successTable.T)
print()
return successTable.T
def makeList(self):
""" This function is under construction. The idea is that when we want to run a large
number of experiments in series, we don't want to wait for the first experiment to
finish if, say, it is still using up 2 processors but there are 2 other processors free.
"""
self.generatePl()
ctr=0
binary = self.bin+'/gidget'
expDir = self.analysis+'/'+self.expName
if(os.path.exists(expDir) and startAt == 0):
print( "Cancelling this run! ",self.expName)
return ([],[],[])
elif(os.path.exists(expDir) and startAt != 0):
pass
else:
os.mkdir(expDir)
stdo=[]
stde=[]
cmds=[]
expDirs=[]
for a_p in self.pl[startAt:]:
ctr += 1
tmpap = a_p[:]
stdo.append(expDir+'/'+a_p[self.keys['name']]+'_stdo.txt')
stde.append(expDir+'/'+a_p[self.keys['name']]+'_stde_aux.txt')
cmds.append([binary]+tmpap[:1]+[repr(el) for el in tmpap[1:]])
expDirs.append(expDir)
return (cmds,stdo,stde,expDirs)
def localRun(self,nproc,startAt,maxTime=7200.0,overwrite=False):
''' Run the specified experiment on this machine,
using no more than nproc processors, and starting
at the "startAt"th run in the experiment '''
self.generatePl()
procs=[]
startTimes=[]
ctr=0
global allProcs
binary=self.bin+'/gidget'
expDir=self.analysis+'/'+self.expName #directory for output files
#if(os.path.exists(expDir) and startAt == 0 and not overwrite):
#print "dbg: ",startAt,overwrite,glob.glob(expDir+'/*comment.txt'),os.path.exists(expDir)
if( len(glob.glob(expDir+'/*comment.txt'))!=0 and startAt == 0 and not overwrite ):
print( "************")
print( "This directory already contains output. CANCELLING this run!")
print()
return
elif(os.path.exists(expDir) and (startAt != 0 or overwrite or len(glob.glob(expDir+'/*comment.txt'))==0)):
pass # let the user overwrite whatever runs they so desire.
else:
os.mkdir(expDir)
for a_p in self.pl[startAt:]:
ctr+=1
tmpap=a_p[:]
globalBolshoiReader.copynearest(expDir, a_p[self.keys['name']]+'_inputRandomFactors.txt', a_p[self.keys['bolshoiWeight']], np.log10(a_p[self.keys['Mh0']]))
try:
shutil.copy('Lacey84_table_K.txt', expDir+'/'+'Lacey84_table_K.txt')
shutil.copy('Lacey84_table_L.txt', expDir+'/'+'Lacey84_table_L.txt')
except:
#pdb.set_trace()
print( "WARNING: failed to copy Lacey84 files. May need to add a set_trace in exper.py::experiment::localRun")
with open(expDir+'/'+a_p[self.keys['name']]+'_stdo.txt','w') as stdo:
with open(expDir+'/'+a_p[self.keys['name']]+'_stde_aux.txt','w') as stde:
print ("Sending run #",ctr,"/",len(self.pl[startAt:])," , ",tmpap[0]," to a local core.")
#print "Parameters: "
#print [binary]+tmpap[:1]+[repr(el) for el in tmpap[1:]]
os.chdir(expDir)
# Here we exclude the final element of tmpap since it corresponds to the variable bolshoiWeight, which
# is only used in the python code, not by the C code.
procs.append(subprocess.Popen([binary]+tmpap[:1]+[repr(el) for el in tmpap[1:-1]],stdout=stdo,stderr=stde))
cmd = ''
for el in tmpap:
cmd+=str(el)+' '
print ("Using cmd: ",cmd)
allProcs.append(procs[-1])
startTimes.append(time.time())
allStartTimes.append(time.time())
nPrinted=True
# we've started a process off and running, but there are
# probably more processes waiting to go. We do not want
# to exceed the number of processors nproc the user
# is willing to let run at once, so let's check what
# our processes are up to:
while True:
# Check whether each process has exceeded the max time:
for k,proc in enumerate(procs):
if proc.poll()==None:
# The process is still running.
if(time.time() - startTimes[k] > maxTime):
print ("WARNING: process has reached maximum allowed time of ",maxTime," seconds! Sending the kill signal.")
print (" If you're sure the process should be running longer, increase the maxTime argument of localRun.")
print (" Usually a run this long means something has gone horribly wrong.")
with open(expDir+'/'+a_p[self.keys['name']]+'_stde_aux.txt','a') as stde:
stde.write('Reached max allowed time '+str(maxTime)+'\n')
proc.kill()
nStillRunning=HowManyStillRunning(allProcs)
nStillRunning2=HowManyStillRunning(procs)
#print "dbg2",nStillRunning,nStillRunning2,expDir
# if our processors are all booked, wait a minute and try again
if(nStillRunning >= nproc):
time.sleep(60) # wait a little bit
if(nPrinted):
print ("Waiting for a free processor...")
nPrinted=False # only print the above message if we haven't seen it since
# a new job has been submitted
# whereas if we have free processors, loop back around to submit the next run.
# If there are no more runs, we're done!
else:
break
def LocalRun(runBundle,nproc):
cmds,stdo,stde,expDirs = runBundle
for i in range(len(cmds)):
with open(stdo[i],'w') as stdo_file:
with open(stde[i],'w') as stde_file:
os.chdir(expDirs[i])
procs.append(subprocess.Popen(cmds[i],stdout=stdo_file,stderr=stde_file))
nPrinted=True
while True:
nStillRunning = HowManyStillRunning(procs)
# if our processors are all booked, wait a minute and try again
if(nStillRunning >= nproc):
time.sleep(10) # wait a little bit
if(nPrinted):
print ("Waiting for a free processor...")
nPrinted=False # only print the above message if we haven't seen it since
# a new job has been submitted
# whereas if we have free processors, loop back around to submit the next run.
# If there are no more runs, we're done!
else:
break
# now all of our processes have been sent off
nPrev = 0
while True:
nStillRunning=HowManyStillRunning(procs)
# has anything changed since the last time we checked?
if(nStillRunning == nPrev and nStillRunning != 0):
# do nothing except wait a little bit
time.sleep(5)
else:
nPrev=nStillRunning
if(nPrev == 0):
break # we're done!
print ("Still waiting for ",nPrev, " processes to finish; I'll check every few seconds for changes.")
print ("Local run complete!")
def GetScaleLengths(N,median=0.045,scatter=0.5,Mh0=1.0e12,sd=100,lower=0.0,upper=1.0e3,multiple=1.0):
''' Return a vector of scale lengths in kpc such that the haloes will have
spin parameters of median, with a given scatter in dex. These
are also scaled by halo mass, since the Virial radius goes as M_h^(1/3).
sd seeds the random number generator so that this function is deterministic.'''
random.seed(sd)
scLengths=[]
while len(scLengths) < N:
spinParameter = median * (10.0 ** random.gauss(0.0,scatter)) # lognormal w/ scatter in dex
if(type(Mh0) == type([1,2,3])): # if Mh0 is a list
if(len(Mh0) == N):
Mh=Mh0[len(scLengths)]
else:
print ("You've given GetScaleLengths the wrong number of halo masses")
Mh=1.0e12
else:
Mh=Mh0
length = multiple * (311.34 * (Mh/1.0e12)**(1.0/3.0) * spinParameter/(2.0**.5))
#length = 2.5*(Mh/1.0e12)**(1./3.) * (spinParameter / .045)
if(length > lower and length < upper):
scLengths.append(length)
return scLengths
def PrintSuccessTables(successTables):
''' Take the results of all experiments named when this script was called
and sum them up such that it's obvious which runs succeeded and which failed.
A schematic example output for a situation where 5 different experiments were
run, each of which varied the same 2 variables, 2 values for 1, 3 for the other:
[ [ 00000, 00000 ], [ 00010, 00010], [00011, 00011] ]
indicates that, as far as whether the runs crash or succeed,
the 2-value variable doesn't make a difference, but for increasing
values of the 3-value variable, fewer runs succeed. The numbers which
appear here are explained in the successCode function above. '''
sumOfSuccessTables = np.zeros(shape=successTables[0].shape,dtype=np.int64)
for i in range(len(successTables)):
table = successTables[i]
sumOfSuccessTables += (table * (10**i))
strSumOfSuccessTables = np.array(sumOfSuccessTables, dtype=str)
it = np.nditer(sumOfSuccessTables, flags=['multi_index'])
while not it.finished:
strSumOfSuccessTables[it.multi_index] = str(it[0]).zfill(len(successTables)-1)
it.iternext()
print()
print ("Here is the sum of all success tables, in homogenized form")
print()
print (strSumOfSuccessTables)
def letter(i, caps=False):
if caps:
return chr(ord("A")+i)
else:
return chr(ord("a")+i)
def NewSetOfExperiments(copyFrom, name, N=1):
if(type(copyFrom)==type([1])):
# if(N!=len(copyFrom)):
# print "WARNING: you asked for ",N,"experiments copied from ",\
# name," containing ",len(copyFrom),"experiments. Ignoring ",N,"."
theList = [copy.deepcopy(copyFrom[i]) for i in range(len(copyFrom))]
else:
theList=[copy.deepcopy(copyFrom) for i in range(N)]
[theList[i].changeName(name+letter(i)) for i in range(len(theList))]
return theList
def experSuiteFromBroadMCMC(list_of_emceeparams, list_of_news, name, list_of_masses, ngal=10, accHistories=None):
#offsets = [ 3.0, 2.0, 2.0, 2.0, 2.0, # raccVir, rstarRed, rgasRed, fg0mult, muColScaling
# 2.0, 10.0, 1.0, 3.0, 0.28, # muFgScaling, muNorm, muMhScaling, ZIGMfac, zmix
# 2.0, 2.0, 2.0, 10.0, 0.28, # eta, Qf, alphaMRI, epsquench, accCeiling
# 0.3, 3.0, 0.23, 2.0, 0.5, # conRF, kZ, xiREC, epsff, initialSlope
# 3.0, 3.0, 0.1] # mquench, enInjFac, chiZslope
#logVars = [1, 1, 1, 1,
# 1, 0, 0, 1,
# 0, 1, 0, 1,
# 1, 1, 1, 0,
# 0, 1, 0, 1,
# 0, 1, 1, 0]
logVars = [1,1,1,1,
0,0,1,0,
1,0,1,1,
1,1,0,0,
1,0,1,0,
1,1,0,0]
# example\:
#q50_11 = [ -0.75225373, 0.45321502, 0.46782298, 0.76024269,
# -0.08407033, 0.29444897, -0.7869128, -0.94169641,
# -0.07408208, 0.18528492, 0.25329861, 0.20123192,
# -1.33669819, -2.75274658, 0.96226019, -0.09557928,
# -1.26589114, 0.54808871, -1.88592761, -0.03303261,
# 12.16116276, -0.51294004, 0.24137744, 1.00588415]
#from broad_svm import lnprior
experList = []
for k in range(len(list_of_emceeparams)):
emceeparams = list_of_emceeparams[k]
news = list_of_news[k]
for i in range(len(emceeparams)):
if logVars[i]==1:
news[i] = 10.0**news[i]
emceeparams[i] = 10.0**emceeparams[i]
experList.append(experFromBroadMCMC(emceeparams, name=name+'fid_k'+str(k), ngal=10, accHistories=accHistories, mass=list_of_masses[k]))
for i in range(len(news)):
thisX = emceeparams[:]
thisX[i] = news[i]
#if logVars[i]==1:
# thisX[i] = thisX[i]*offsets[i]
#else:
# thisX[i] = thisX[i]+offsets[i]
#if not np.isfinite(lnprior(thisX)):
#if logVars[i]==1:
# thisX[i] = thisX[i]/offsets[i]/offsets[i]
#else:
# thisX[i] = thisX[i] - 2.0*offsets[i]
#print "WARNING: Going the opposite direction for variable ",i
#if not np.isfinite(lnprior(thisX)):
# pdb.set_trace() # tried adjusting in both directions and still ended up outside the allowed range?!
experList.append(experFromBroadMCMC(thisX, name=name+'_k'+str(k)+'_'+chr(65+i), ngal=10, accHistories=accHistories, mass=list_of_masses[k]))
return experList
def experFromBroadMCMC(emceeparams, name=None, ngal=30, accHistories=None, mass=12.0):
raccRvir, rstarRed, rgasRed, fg0mult, muColScaling, muFgScaling, muNorm, muMhScaling, ZIGMfac, zmix, eta, Qf, alphaMRI, epsquench, accCeiling, conRF, kZ, xiREC, epsff, initialSlope, mquench, enInjFac, chiZslope = emceeparams
if hasattr( mass, '__len__' ):
assert np.all(mass>1.0e5)
Mhz0 = list(mass)
else:
assert mass<20
Mhz0 = list(np.power(10.0, np.linspace(mass, mass+0.001,ngal)))
# create experiment
assert not name is None, "Need a name!"
#if name is None:
#name = basename+str(runnumber).zfill(3)+'_'+str(rank).zfill(5)+'_'+str(proccounter).zfill(5)
thisexper = experiment(copy.copy(name))
# we need to put the random factors for the accretion history into a file for gidget to read.
#directory = chaindir[0:-len(chaindirrel)]+'/'+name+'/'
#if not os.path.exists(directory):
# os.makedirs(directory)
#proccounter+=1
#### Paramters to vary:
# Mhz0, raccRvir, rstarRed, rgasRed, fg0mult, muColScaling, muFgScaling, muNorm, ZIGMfac, zmix, eta, Qf, alphaMRI, epsquench, accCeiling, conRF, kZ, xiREC
thisexper.irregularVary('kappaMetals',kZ)
thisexper.irregularVary('xiREC',xiREC)
Mhz0 = np.array(Mhz0)
Mhz4 = Mhz0/(.75*33.3 * np.power(Mhz0/1.5e13,0.359)) # very roughly... We know for Mh0 = 3e11, Mhz4 =4e10, and same for Mh0=2e13->Mhz4=6e11. Using those 4 numbers and assuming a powerlaw, we get this relation.
def Moster(Mh, mparams):
M10, M11, N10, N11, beta10, beta11, gamma10, gamma11 = mparams
zti = 4.0
logM1z = M10 + M11*zti/(zti+1.0)
Nz = N10 + N11*zti/(zti+1.0)
betaz = beta10 + beta11*zti/(zti+1.0)
gammaz = gamma10 + gamma11*zti/(zti+1.0)
M1 = np.power(10.0, logM1z)
eff = 2.0*Nz / (np.power(Mh/M1,-betaz) + np.power(Mh/M1,gammaz))
return eff
central = np.array([11.590, 1.195, 0.0351, -0.0247, 1.376+initialSlope, -0.826, 0.608, 0.329])
eff = Moster(Mhz4,central)
mst = eff*Mhz4 # mstar according to the moster relation. ## at z=4 !!
f0 = 1.0/(1.0 + np.power(mst/10.0**9.15,0.4)) # from Hayward & Hopkins (2015) eq. B2
tau4 = 12.27/(12.27+1.60) # fractional lookback time at z=4
fgz4 = f0*np.power(1.0 - tau4*(1.0-np.power(f0,1.5)), -2.0/3.0)
reff4 = 5.28*np.power(mst/1.0e10, 0.25)*np.power(1.0+4.0,-0.6) # kpc (eq B3) at z=4
reff411 = 5.28*np.power(1.0e11/1.0e10, 0.25)*np.power(1.0+4.0,-0.6) # kpc (eq B3) at z=4
weight1 = np.exp(-Mhz0/1.0e12) #the ad-hoc reductions in fg and fcool make things bad for high-mass galaxies. weight1 is for tiny galaxies.
weight2 = 1-weight1
rat4 = 1.0/(1.0/fgz4 - 1.0)
rat4 *= fg0mult # double the gas content
# even though we give fg0 to the gidget here, it is overridden later by giving gidget fg0mult and dbg switch 10, which internally computes in gidget's own initialization what the fg0 should be. Essentially this is because at this moment in the python code we don't actually know what the halo mass is at the initial redshift, and hence we don't know what the stellar mass is, so we don't know where we are in the "HH15" relation of fg vs. stellar mass. Within gidget we redo this ^ calculation with the actual initial halo,stellar masses.
fgUsed =1.0/(1.0/rat4 + 1.0)*(1.0*weight1+1.0*weight2)
if not np.all(fgUsed<=1.0):
pdb.set_trace()
thisexper.irregularVary( 'fg0', list(fgUsed), 5)
thisexper.irregularVary( 'Mh0', list(Mhz0), 5)
fcools = 1.0/(1.0-fgUsed) * mst/(0.17*Mhz4) * (1.0*weight1+1.0*weight2) # The factor of 0.3 is a fudge factor to keep the galaxies near Moster for longer, and maybe shift them to be a bit more in line with what we would expect for e.g. the SF MS.
if not np.all(fcools<=1.0):
print ("WARNING: ", name, "fcools exceed unity: ", fcools)
thisexper.irregularVary('fcool', list(fcools), 5)
thisexper.irregularVary('NChanges', 1001)
width = 0.001
asls = raccRvir*np.power(10.0, np.random.normal(0,1,len(Mhz0))*width) # * np.power(Mhz0/1.0e12,alpharmh)
#if asls<0.005:
# asls=0.005
asls = np.clip(asls, 0.005, np.inf)
thisexper.irregularVary('accScaleLength', list(asls), 5)
thisexper.irregularVary( 'R', list(np.power(reff4/reff411, 1.0)*50* asls/0.042) , 5)
thisexper.irregularVary('xmin',0.001)
if accHistories is None:
bolweights = list( np.random.random(len(Mhz0)) )
else:
assert len(Mhz0) == len(accHistories)
bolweights = list(accHistories[:])
thisexper.irregularVary('bolshoiWeight', bolweights ,5)
thisexper.irregularVary('Noutputs',200) ## why so many??
thisexper.irregularVary('zstart',3.98)
thisexper.irregularVary('zrelax',4.0)
thisexper.irregularVary('tauHeat', 1.0)
thisexper.irregularVary('muNorm', muNorm)
thisexper.irregularVary('muFgScaling', muFgScaling)
thisexper.irregularVary('muColScaling', muColScaling)
thisexper.irregularVary('muMhScaling', muMhScaling)
thisexper.irregularVary('fscatter', 1.0)
thisexper.irregularVary('accCeiling',accCeiling)
thisexper.irregularVary('NPassive',10)
thisexper.irregularVary('eta',eta)
thisexper.irregularVary('epsff',epsff)
thisexper.irregularVary('yREC',0.03)
thisexper.irregularVary( 'nx', 256 ) # FFT o'clock.
ZLee = np.power(10.0, 5.65 + chiZslope*np.log10(mst/1.0e10) + 0.3*np.log10(1.0e10) - 8.7 ) # Z in Zsun
thisexper.irregularVary('ZIGM', list(ZIGMfac*0.05*ZLee*0.02), 5)
thisexper.irregularVary('concentrationRandomFactor', conRF) ## units of dex! # 0.7
thisexper.irregularVary('ZIGMfac', ZIGMfac)
thisexper.irregularVary('chiZslope', chiZslope)
thisexper.irregularVary('deltaBeta', initialSlope)
thisexper.irregularVary('fg0mult', fg0mult)
thisexper.irregularVary('whichAccretionHistory', -112)
thisexper.irregularVary('ksuppress', 10.0 )
thisexper.irregularVary('kpower', 2.0)
thisexper.irregularVary('dbg', 2**0 + 2**1 + 2**2 + 2**4 + 2**5 + 2**6 + 2**7 + 2**10 + 2**16 )
thisexper.irregularVary('fscatter', 1.0)
thisexper.irregularVary('MQuench', mquench)
thisexper.irregularVary('epsquench', epsquench)
thisexper.irregularVary('muQuench', 0.0)
thisexper.irregularVary('gaScaleReduction', rgasRed )
thisexper.irregularVary('stScaleReduction', rstarRed )
thisexper.irregularVary('ZMix', zmix)
thisexper.irregularVary('alphaMRI', alphaMRI)
thisexper.irregularVary('Qlim', Qf+0.05)
thisexper.irregularVary('fixedQ', Qf)
thisexper.irregularVary('TOL', 6.0e-4)
thisexper.irregularVary('energyInjectionFactor', enInjFac)
thisexper.irregularVary('minSigSt', 10.0)
thisexper.irregularVary('AccretionHeatingRate', 0.0)
return thisexper
if __name__ == "__main__":
# The structure here is straightforward:
# 1) First, parse the arguments
# 2) Then define a bunch of experiments
# 3) Run the experiments the user told us to run in the command line arguments
# Read in the arguments. To see the results of this bit of code, just try to run
# this script without any arguments, i.e. $ python exper.py
parser = argparse.ArgumentParser(description='Analyze data and/or run experiments associated with a list of experiment names.')
parser.add_argument('models', metavar='experiment', type=str, nargs='+',
help='a string contained in any experiment to be run, e.g. if rs04a rs04b rs05a are the models defined, then exper.py rs04 will run the first two, exper.py a will run the first and last, exper.py rs will run all three, exper.py rs04b rs05 will run the last two, etc.')
parser.add_argument('--nproc',type=int,help="maximum number of processors to use (default: 16)",default=16)
parser.add_argument('--start',metavar='startingModel',type=int,
help='The number of the model in the experiment (as ordered by GeneratePl) at which we will start sending experiments to the processors to run. (default: 0)',default=0)
parser.add_argument('--xgrid',type=bool,help="run on an xGrid (requires the user to submit the generated file to an xGrid (default: False)",default=False)
args = parser.parse_args()
# Store the names of the experiments the user has told us to run.
modelList=[]
for modelName in args.models:
modelList.append(modelName)
# Begin defining experiments!
re01=experiment("re01")
re01.irregularVary("R",40)
re01.irregularVary('accScaleLength', .05) # now a fraction of r200
re01.irregularVary('NPassive',3)
re01.irregularVary('dbg',2)
re01.irregularVary('xmin',.01)
re01.irregularVary('alphaMRI',.01)
re01.irregularVary('fcool',0.6)
re01.irregularVary('nx',200)
re01.irregularVary('kappaMetals',1.0)
re01.irregularVary('xiREC',0.0)
re01.irregularVary('muNorm',1)
re02=NewSetOfExperiments(re01,'re02')
re02[0].irregularVary('muNorm',[.1,1,10])
re03=NewSetOfExperiments(re01,'re03')
re03[0].irregularVary('accScaleLength', [.01,.03,.1])
re04=NewSetOfExperiments(re01,'re044')
re04[0].irregularVary('NPassive',1)
re04[0].irregularVary('zstart',[4.8],3)
re04[0].irregularVary('zrelax',[5.0],3)
re04[0].irregularVary('accScaleLength',[.01,.031,.1,.25],4)
re04[0].irregularVary('R',[15,20,35,40],4)
re04[0].irregularVary('alphaMRI',0)
re04[0].irregularVary('muNorm',0.0)
re05=NewSetOfExperiments(re01,'re05')
re05[0].irregularVary('zstart',[4.9],3)
re05[0].irregularVary('zrelax',[5.0],3)
re05[0].irregularVary('accScaleLength',[.01,.031,.1],4)
re05[0].irregularVary('R',[15,30,45],4)
re05[0].irregularVary('muNorm',[.1,1,10])
re05[0].irregularVary('eta',[.5,4.5])
re05[0].irregularVary('fixedQ',[1.3,2.0,3.0])
re05[0].irregularVary('fcool',[.1,.9])
re05[0].irregularVary('fg0',[.3,.6,.9])
re05[0].irregularVary('Mh0',[3.0e11,1.0e12,3.0e12])
# Playing around with mass dependence
re06=NewSetOfExperiments(re01,'re06666')
re06[0].irregularVary('NPassive',1)
re06[0].irregularVary('zstart',[4.99],3)
re06[0].irregularVary('zrelax',[5.0],3)
re06[0].irregularVary('accScaleLength',.05)
re06[0].vary('Mh0',1.0e10,1.0e13,9,1,4)
re06[0].vary('R',5,50,9,1,4)
re06[0].irregularVary('fg0',.9999)
re06[0].irregularVary('fcool',.05)
re06[0].irregularVary('alphaMRI',0)
re06[0].irregularVary('muNorm',0.0)
re06[0].irregularVary('ZIGM',.00002)
re06[0].irregularVary('whichAccretionProfile', 2)
re06[0].irregularVary('widthAccretionProfile', 0.6)
re06[0].irregularVary('dbg',2**4+2**1+2**0)
# Adding some stochasticity before messing around with previous expt.
re07=NewSetOfExperiments(re01,'re07')
re07[0].irregularVary('NPassive',1)
re07[0].irregularVary('zstart',[4.8],3)
re07[0].irregularVary('zrelax',[5.0],3)
re07[0].irregularVary('accScaleLength',.05)
re07[0].vary('Mh0',1.0e10,1.0e13,200,1,4)
re07[0].vary('R',5,50,200,1,4)
re07[0].vary('whichAccretionHistory',-300,-101,0,4)
re07[0].irregularVary('alphaMRI',0)
re07[0].irregularVary('muNorm',0.0)
# Use params from toy model.
re08=NewSetOfExperiments(re01,'re08')
re08[0].irregularVary('NPassive',1)
re08[0].irregularVary('zstart',[4.99],3)
re08[0].irregularVary('zrelax',[5.0],3)
re08[0].irregularVary('accScaleLength',.05)
re08[0].vary('Mh0',1.0e10,1.0e13,9,1,4)
re08[0].vary('R',5,50,9,1,4)
re08[0].irregularVary('fg0',.9999)
re08[0].irregularVary('fcool',.05)
re08[0].irregularVary('alphaMRI',0)
re08[0].vary('muNorm',22,.21,9,1,4)
re08[0].irregularVary('ZIGM',.00002)
re08[0].irregularVary('dbg',2**4+2**1+2**0)
re09=NewSetOfExperiments(re08[0],'re09')
re09[0].irregularVary('dbg',2**4+2**1+2**0+2**12)
# Modified the code to use mu = .5 (Mh/10^12)^(-2/3)
re10=NewSetOfExperiments(re08[0],'re10')
re10[0].irregularVary('accCeiling',0.6)
re10[0].vary('R',15,100,9,1,4)
re10[0].irregularVary('accScaleLength',0.1)
re10[0].vary('fcool',0.05/2,0.05,9,1,4)
#re10[0].irregularVary('accNorm',.1)
#re10[0].irregularVary('dbg',2**4+2**1+2**0+2**
# sample to study MS.
re11=NewSetOfExperiments(re01,'re11')
re11[0].vary('whichAccretionHistory',-300,-101,200,0,4)
re11[0].vary('R', 15, 100, 200, 1, 4)
re11[0].vary('fcool', 0.05/2, 0.05, 200, 1, 4)
re11[0].vary('Mh0',1.0e10, 1.0e13, 200, 1, 4)
re11[0].vary('muNorm', 22, .21, 200, 1, 4)
re11[0].irregularVary('fscatter', .45)
re11[0].irregularVary('NChanges', list(np.random.binomial(40,.25,size=200)) ,4)
re11[0].irregularVary('accCeiling',0.6)
re11[0].irregularVary('accScaleLength',0.1)
re11[0].irregularVary('NPassive',1)
re11[0].irregularVary('zstart',[4.99],3)
re11[0].irregularVary('zrelax',[5.0],3)
re11[0].irregularVary('fg0',.9999)
re11[0].irregularVary('alphaMRI',0)
re11[0].irregularVary('ZIGM',.00002)
re11[0].irregularVary('dbg',2**4+2**1+2**0)
# more reasonable ZIGM - are galaxies more in equilibrium?
re12=NewSetOfExperiments(re10[0], 're12')
re12[0].irregularVary('ZIGM',.002)
re12[0].irregularVary('accCeiling',0.4)
re12[0].irregularVary('accScaleLength',0.1)
# MW progenitors
#re13=NewSetOfExperiments(re01[0], 're13')
#
re14=NewSetOfExperiments(re01, 're14')
re14[0].irregularVary('xmin', [.00031, .001, .0031, .01, .031])
re14[0].irregularVary('zstart',[4.95],3)
re14[0].irregularVary('zrelax',[5.0],3)
re14[0].irregularVary('dbg',2**4+2**1+2**0)
re15 = NewSetOfExperiments(re01,'re15')
re15[0].irregularVary('zstart',[4.95],3)
re15[0].irregularVary('zrelax',[5.0],3)
re15[0].irregularVary('dbg',2**4+2**1+2**0)
re15[0].irregularVary('Mh0',[1.0e11,1.0e12,1.0e13,1.0e14])
re15[0].irregularVary('accScaleLength',[.01,.03,.05,.1,.2])
re15[0].irregularVary('concentrationRandomFactor',[-.3,-.1,0.0,.1,.3])
#from mcmc_individual import emceeParameterSpaceToGidgetExperiment
#favoriteParams = [ 2.83269124e-01, 9.30881877e-03, 1.82351732e-01, 1.50469698e+00,
# 3.24674422e-01, 1.85339104e+00, 7.43781644e-03, 1.94868075e-01,
# 4.05944370e+12, 1.38101233e-01, -2.36399004e-01, -5.50757004e-01,
# 5.94584202e-01, -6.06652441e-01, 8.88202578e+00, 3.92808177e-01,
# -2.84901649e+00]
#re16 = emceeParameterSpaceToGidgetExperiment(favoriteParams,'re16')[0]
#re16.irregularVary("Noutputs",200)
#allModels['re16']=re16
# Barro plot...
re17=NewSetOfExperiments(re01,'re17')
re17[0].irregularVary('zstart',4.95)
re17[0].irregularVary('zrelax',5.0)
re17[0].irregularVary('dbg',2**4+2**1+2**0)
re17[0].vary('whichAccretionHistory',-300,-101,200,0,4)
re17[0].irregularVary('R', 10)
re17[0].irregularVary('fcool', .8)
re17[0].vary('Mh0',1.0e12, 1.0e13, 200, 1, 4)
re17[0].irregularVary('muNorm', 0.5)
re17[0].irregularVary('fscatter', .45)
re17[0].irregularVary('NChanges', list(np.random.binomial(40,.25,size=200)) ,4)
re17[0].irregularVary('accCeiling',0.6)
re17[0].irregularVary('accScaleLength',0.1)
re17[0].irregularVary('NPassive',1)
re17[0].irregularVary('eta',0.5)
re18=NewSetOfExperiments(re17,'re18')
re18[0].irregularVary('R',30)
re18[0].irregularVary('Noutputs',300)
re18[0].vary('Mh0',3.0e12, 1.0e13, 200, 1, 4)
re18[0].irregularVary('yREC',.02)
re18[0].irregularVary('kappaMetals',.5)
re19a=NewSetOfExperiments(re18,'re19a')
re19a[0].irregularVary('accScaleLength', .01)
re19a[0].irregularVary('R',10)
re19b=NewSetOfExperiments(re18,'re19b')
re19b[0].irregularVary('accScaleLength', .04)
re19b[0].irregularVary('R',20)
re19c=NewSetOfExperiments(re18,'re19c')
re19c[0].irregularVary('accScaleLength', .09)
re19c[0].irregularVary('R',30)
re20=NewSetOfExperiments(re01,'re20')
re20[0].irregularVary('Noutputs',600)
re20[0].irregularVary('zstart',3.95)
re20[0].irregularVary('zrelax',4.0)
re20[0].irregularVary('dbg',2**4+2**1+2**0 )
asls = np.array([.012,.027,.042,.057,.072])
re20[0].irregularVary('accScaleLength',list(asls), 5)
re20[0].irregularVary('R', list(asls * 20/0.012), 5)
re20[0].irregularVary('fcool', .8)
re20[0].irregularVary('Mh0', [1.0e13])
re20[0].irregularVary('muNorm', 1.5)
re20[0].irregularVary('muFgScaling', 0.4)
re20[0].irregularVary('muColScaling', 0)
re20[0].irregularVary('fscatter', .45)
re20[0].irregularVary('accCeiling',0.6)
re20[0].irregularVary('NPassive',1)
re20[0].irregularVary('eta',0.5)
re21=NewSetOfExperiments(re20,'re21')
re21[0].irregularVary('dbg', 2**4+2**1+2**0 + 2**12 )
re22=NewSetOfExperiments(re20,'re22')
re22[0].irregularVary('accScaleLength', list(asls*4), 5)
re22[0].irregularVary('fixedQ',1.5)
re22[0].irregularVary('whichAccretionProfile', 2)
re22[0].irregularVary('widthAccretionProfile', 0.5)
re22[0].irregularVary('TOL',1.0e-4)
re23=NewSetOfExperiments(re20,'re23')
re23[0].irregularVary('Qlim',.05)
re23[0].irregularVary('fixedQ',.05)
re30=NewSetOfExperiments(re20,'re30')
re30[0].irregularVary('Noutputs',600)
re30[0].irregularVary('zstart',3.95)
re30[0].irregularVary('zrelax',4.0)
re30[0].irregularVary('dbg',2**4+2**1+2**0 + 2**14 )
asls = np.array([.012,.027,.042,.057,.072])
re30[0].irregularVary('accScaleLength',.042 )
mhl = np.power(10.0, np.linspace(10,13,200))
re30[0].irregularVary('R', list(.042 * 20/0.012 * np.power(mhl/1.0e12,1.0/3.0) ), 6)
re30[0].irregularVary('fcool', .8)
re30[0].irregularVary('Mh0', list(mhl), 6)
re30[0].irregularVary('muNorm', list(1.5*np.power(mhl/1.0e12, -2.0/3.0)), 6)
re30[0].irregularVary('muFgScaling', 0.4)
re30[0].irregularVary('muColScaling', 0)
re30[0].irregularVary('fscatter', .45)
re30[0].irregularVary('accCeiling',0.6)