Update to Brent and multi-param to put in bounds. Update to Brent and findSim to handle presettling.

Author: upibhalla

Brent_minimization.py

@@ -85,20 +85,30 @@ def enumerateFindSimFiles( location ):
         quit()
 
 class EvalFunc:
-    def __init__( self, objField, expts, weights, pool, modelFile ):
+    def __init__( self, objField, expts, weights, pool, modelFile, presettle = [] ):
         self.objField = objField
         self.expts = expts
         self.weights = weights
         self.pool = pool # pool of available CPUs
         self.modelFile = modelFile
+        self.presettle = presettle
 
     def doEval( self, x ):
         ret = []
         spl = self.objField.split( '.' )
         assert( len(spl) == 2 )
         obj, field = spl
+
+        settleDict = {}
+        if len( self.presettle ) == 3:
+            presettleTime = float( self.presettle[2] )
+            if presettleTime > 0:
+            #print("{}".format( self.presettle ) )
+                settleDict = findSim.innerMain( self.presettle[0], modelFile = self.presettle[1], hidePlot=True, silent=True, scaleParam=[obj,field,str(x)], settleTime = presettleTime )
+        #print( "Doing presettle, len = {}".format( len(settleDict) ) )
+
         for k in self.expts:
-            ret.append( self.pool.apply_async( findSim.innerMain, (k,), dict(modelFile = self.modelFile, hidePlot=True, silent=True, scaleParam=[obj,field,str(x)]), callback = reportReturn ) )
+            ret.append( self.pool.apply_async( findSim.innerMain, (k,), dict(modelFile = self.modelFile, hidePlot=True, silent=True, scaleParam=[obj,field,str(x)], settleDict=settleDict ), callback = reportReturn ) )
         score = [ i.get() for i in ret ]
         sumScore = sum([ s*w for s,w in zip(score, self.weights) if s>=0.0])
         sumWts = sum( [ w for s,w in zip(score, self.weights) if s>=0.0 ] )
@@ -111,8 +121,10 @@ def main():
     parser.add_argument( 'location', type = str, help='Required: Directory in which the scripts (in tsv format) are all located. OR: File in which each line is the filename of a scripts.tsv file, followed by weight to assign for that file.')
     parser.add_argument( '-n', '--numProcesses', type = int, help='Optional: Number of processes to spawn', default = 2 )
     parser.add_argument( '-m', '--model', type = str, help='Optional: Composite model definition file. First searched in directory "location", then in current directory.', default = "FindSim_compositeModel_1.g" )
-    parser.add_argument( '-p', '--parameter_sweep', nargs='*', default=[],  help='Does a parameter sweep in range 0.5-2x of each object.field pair.' )
-    parser.add_argument( '-f', '--file', type = str, help='Optional: File name for output of parameter sweep', default = "" )
+    parser.add_argument( '-p', '--parameter_optimize', nargs='*', default=[],  help='Does a parameter optimization for each specified object.field pair.' )
+    parser.add_argument( '-ps', '--presettle', nargs=3, default=[],  help='Arguments: tsv_file, model_file, settle_time. Obtains values of all concentrations after a specified settle-time, so that all calculations for the optimization runs can be initialized to this presettled value. The tsv_file is to specify which subset of the model_file to use. This option is most useful in costly multiscale models.' )
+    parser.add_argument( '-f', '--file', type = str, help='Optional: File name for output of parameter optimization', default = ""
+    )
     args = parser.parse_args()
     location = args.location
     if location[-1] != '/':
@@ -130,13 +142,14 @@ def main():
     pool = Pool( processes = args.numProcesses )
 
     results = {}
-    for i in args.parameter_sweep:
+    for i in args.parameter_optimize:
         print( "{}".format( i ) )
         spl = i.split( '.' )
         assert( len(spl) == 2 )
         obj, field = spl
-        ev = EvalFunc( i, fnames, weights, pool, modelFile )
-        results[i] = optimize.minimize_scalar( ev.doEval )
+        ev = EvalFunc( i, fnames, weights, pool, modelFile, args.presettle )
+        # Bounded method uses Bounded Brent method.
+        results[i] = optimize.minimize_scalar( ev.doEval, method = 'bounded', bounds = (0.0, 100.0) )
         print( "\n Finished optimizing for " + i)
     print( "\n---------------- Completed ----------------- " )
     dumpData = False

findSim.py

@@ -33,7 +33,7 @@
 **********************************************************************/
 
 '''
-
+from __future__ import print_function
 import heapq
 import pylab
 import numpy as np
@@ -545,7 +545,10 @@ def _scaleOneParam( self, params ):
 
         obj = self.findObj( '/model', params[0] )
         scale = float( params[2] )
-        assert( scale >= 0.0 and scale <= 100.0 )
+        if not ( scale >= 0.0 and scale <= 100.0 ):
+            print( "Error: Scale {} out of range".format( scale ) )
+            assert( False )
+        #assert( scale >= 0.0 and scale <= 100.0 )
         if params[1] == 'Kd':
             if not obj.isA[ "ReacBase" ]:
                 raise SimError( "scaleParam: can only assign Kd to a Reac, was: '{}'".format( obj.className ) )
@@ -1435,11 +1438,19 @@ def main():
     parser.add_argument( '-hs', '--hide_subplots', action="store_true", help='Hide subplot output of simulation. By default the graphs include dotted lines to indicate individual quantities (e.g., states of a molecule) that are being summed to give a total response. This flag turns off just those dotted lines, while leaving the main plot intact.' )
     parser.add_argument( '-o', '--optimize_elec', action="store_true", help='Optimize electrical computation. By default the electrical computation runs for the entire duration of the simulation. With this flag the system turns off the electrical engine except during the times when electrical stimuli are being given. This can be *much* faster.' )
     parser.add_argument( '-s', '--scale_param', nargs=3, default=[],  help='Scale specified object.field by ratio.' )
+    parser.add_argument( '-settle_time', '--settle_time', type=float, default=0,  help='Run model for specified settle time and return dict of {path,conc}.' )
     args = parser.parse_args()
-    innerMain( args.script, modelFile = args.model, dumpFname = args.dump_subset, paramFname = args.param_file, hidePlot = args.hide_plot, hideSubplots = args.hide_subplots, optimizeElec = args.optimize_elec, scaleParam = args.scale_param )
+    innerMain( args.script, modelFile = args.model, dumpFname = args.dump_subset, paramFname = args.param_file, hidePlot = args.hide_plot, hideSubplots = args.hide_subplots, optimizeElec = args.optimize_elec, scaleParam = args.scale_param, settleTime = args.settle_time )
+
 
+def innerMain( script, modelFile = "model/synSynth7.g", dumpFname = "", paramFname = "", hidePlot = True, hideSubplots = False, optimizeElec=True, silent = False, scaleParam=[], settleTime = 0, settleDict = {} ):
+    ''' If *settleTime* > 0, then we need to return a dict of concs of
+    all variable pools in the chem model obtained after loading in model, 
+    applying all modifications, and running for specified settle time.\n
+    If the *settleDict* is not empty, then the system goes through and 
+    matches up pools to assign initial concentrations.
+    '''
 
-def innerMain( script, modelFile = "model/synSynth7.g", dumpFname = "", paramFname = "", hidePlot = True, hideSubplots = False, optimizeElec=True, silent = False, scaleParam=[] ):
     global pause
     solver = "gsl"  # Pick any of gsl, gssa, ee..
     modelWarning = ""
@@ -1511,6 +1522,33 @@ def innerMain( script, modelFile = "model/synSynth7.g", dumpFname = "", paramFna
             for i in range( 10, 20 ):
                 moose.setClock( i, 0.1 )
 
+        ##############################################################
+        # Here we handle presettling. First to generate, then to apply
+        # the dict of settled values.
+        if settleTime > 0:
+            t0 = time.time()
+            moose.reinit()
+            #print settleTime
+            moose.start( settleTime )
+            w = moose.wildcardFind( modelId.path + "/##[ISA=PoolBase]" )
+            ret = {}
+            for i in w:
+                if not i.isBuffered:
+                    ret[i.path] = i.n
+                    #print( "{}.nInit =   {:.3f}".format( i.path, i.n ))
+            #print "-------------------- settle done -------------------"
+            moose.delete( modelId )
+            if moose.exists( '/library' ):
+                moose.delete( '/library' )
+            #print( "Done settling in {:.2f} seconds".format( time.time()-t0))
+            print( "s", end = '' )
+            sys.stdout.flush()
+            return ret
+
+        for key, value in settleDict.items():
+            moose.element( key ).nInit = value
+        ##############################################################
+
         t0 = time.time()
         score = runit( expt, model,stims, readouts, modelId )
         elapsedTime = time.time() - t0
@@ -1522,13 +1560,17 @@ def innerMain( script, modelFile = "model/synSynth7.g", dumpFname = "", paramFna
 
             pylab.show()
         moose.delete( modelId )
+        if moose.exists( '/library' ):
+            moose.delete( '/library' )
         return score
 
     except SimError as msg:
         if not silent:
             print( "Error: findSim failed for script {}: {}".format(script, msg ))
         if modelId:
             moose.delete( modelId )
+            if moose.exists( '/library' ):
+                moose.delete( '/library' )
         return -1.0
 # Run the 'main' if this script is executed standalone.
 if __name__ == '__main__':

multi_param_minimization.py

@@ -54,7 +54,7 @@
 
 resultCount = 0
 numIterations = 0
-TOLERANCE = 0.001 # Don't want to go too tight on tolarance.
+TOLERANCE = 0.0001 # Don't want to go too tight on tolarance.
 
 def reportReturn( result ):
     global resultCount
@@ -146,13 +146,18 @@ def main():
     pool = Pool( processes = args.numProcesses )
 
     params = []
+    bounds = []
     for i in args.parameters:
         print( "{}".format( i ) )
         spl = i.split( '.' )
         assert( len(spl) == 2 )
         params.append( i )
+        if spl[1] == 'Kd' or spl[1] == 'tau' or spl[1] == 'Km':
+            bounds.append( (0.03,30) )
+        else:
+            bounds.append( (0.0, 30 ) ) # Concs, Kfs and Kbs can be zero.
     ev = EvalFunc( params, fnames, weights, pool, modelFile )
-    results = optimize.minimize( ev.doEval, np.ones( len(params) ), tol = TOLERANCE, callback = optCallback )
+    results = optimize.minimize( ev.doEval, np.ones( len(params) ), method='L-BFGS-B', tol = TOLERANCE, callback = optCallback, bounds = bounds )
     print( "\n----------- Completed in {:.3f} sec ---------- ".format(time.time() - t0 ) )
     print( "\n----- Score= {:.4f} ------ ".format(results.fun ) )
     dumpData = False