@@ -25,16 +25,21 @@ import io.github.mandar2812.dynaml.utils
2525import scala .util .Random
2626
2727/**
28+ * Skeleton implementation of the Coupled Simulated Annealing algorithm
2829 * @author mandar datum 01/12/15.
29- *
30- * Implementation of the Coupled Simulated Annealing algorithm
31- * for global optimization.
3230 */
33- class AbstractCSA [M <: GloballyOptimizable ](model : M )
34- extends AbstractGridSearch [M ](model : M ){
31+ abstract class AbstractCSA [M <: GloballyOptimizable , M1 ](model : M )
32+ extends AbstractGridSearch [M , M1 ](model : M ) {
3533
3634 protected var MAX_ITERATIONS : Int = 10
3735
36+ protected var variant = AbstractCSA .MuSA
37+
38+ def setVariant (v : String ) = {
39+ variant = v
40+ this
41+ }
42+
3843 def setMaxIterations (m : Int ) = {
3944 MAX_ITERATIONS = m
4045 this
@@ -55,104 +60,185 @@ class AbstractCSA[M <: GloballyOptimizable](model: M)
5560 this
5661 }
5762
58- protected val acceptance = (energy : Double ,
59- coupling : Double ,
60- temperature : Double ) => {
61- val prob = math.exp(- 1.0 * energy/ temperature)
62- prob/ (prob+ coupling)
63- }
63+ var iTemp = 1.0
6464
65- protected val mutate = (config : Map [String , Double ], temperature : Double ) => {
66- config.map((param) => {
67- val dist = new CauchyDistribution (0.0 , temperature)
68- val mutated = param._2 + dist.sample()
69- (param._1, math.abs(mutated))
70- })
71- }
65+ var alpha = 0.05
66+
67+ private def computeDesiredVariance = AbstractCSA .varianceDesired(variant) _
68+
69+ private def computeCouplingFactor = AbstractCSA .couplingFactor(variant) _
70+
71+ private def computeAcceptanceProb = AbstractCSA .acceptanceProbability(variant) _
72+
73+ protected val mutate : (Map [String , Double ], Double ) => Map [String , Double ] =
74+ (config : Map [String , Double ], temperature : Double ) => {
75+ logger.info(" Mutating configuration: \n " + GlobalOptimizer .prettyPrint(config)+ " \n " )
76+
77+ config.map((param) => {
78+ val dist = new CauchyDistribution (0.0 , temperature)
79+ val mutated = param._2 + dist.sample()
80+ (param._1, math.abs(mutated))
81+ })
82+ }
7283
7384 def acceptanceTemperature (initialTemp : Double )(k : Int ): Double =
7485 initialTemp/ math.log(k.toDouble+ 1.0 )
7586
7687 def mutationTemperature (initialTemp : Double )(k : Int ): Double =
7788 initialTemp/ k.toDouble
7889
79- override def optimize (initialConfig : Map [String , Double ],
80- options : Map [String , String ] = Map ()) = {
90+ protected def performCSA (
91+ initialConfig : Map [String , Double ],
92+ options : Map [String , String ] = Map ()): List [(Double , Map [String , Double ])] = {
93+
94+ logger.info(" -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-" )
95+ logger.info(" Coupled Simulated Annealing (CSA): " + AbstractCSA .algorithm(variant))
96+ logger.info(" -+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-" )
8197
82- // create grid
83- val iTemp = 1.0
8498 var accTemp = iTemp
8599 var mutTemp = iTemp
86- // one list for each key in initialConfig
87- val hyper_params = initialConfig.keys.toList
88- val scaleFunc = if (logarithmicScale) (i : Int ) => math.exp((i+ 1 ).toDouble* step) else
89- (i : Int ) => (i+ 1 ).toDouble* step
90100
91- val gridvecs = initialConfig.map((keyValue) => {
92- (keyValue._1, List .tabulate[Double ](gridsize)(scaleFunc))
93- })
101+ // Calculate desired variance
102+ val sigmaD = computeDesiredVariance(math.pow(gridsize, initialConfig.size).toInt)
94103
95- val grid = utils.combine(gridvecs.map(_._2)).map(x => DenseVector (x.toArray))
96-
97- val energyLandscape = grid.map((config) => {
98- val configMap = List .tabulate(config.length){i => (hyper_params(i), config(i))}.toMap
99- logger.info(" Evaluating Configuration: " + configMap)
100-
101- val configEnergy = system.energy(configMap, options)
102-
103- logger.info(" Energy = " + configEnergy+ " \n " )
104-
105- (configEnergy, configMap)
106- }).toMap
107-
108- var currentEnergyLandscape = energyLandscape
109- var newEnergyLandscape = energyLandscape
110-
111- (1 to MAX_ITERATIONS ).foreach((iteration) => {
112- logger.info(" **************************" )
113- logger.info(" CSA Iteration: " + iteration)
114- // mutate each element of the grid with
115- // the generating distribution
116- // and accept using the acceptance distribution
117- mutTemp = mutationTemperature(iTemp)(iteration)
118- accTemp = acceptanceTemperature(iTemp)(iteration)
119- val couplingFactor = currentEnergyLandscape.map(c => math.exp(- 1.0 * c._1/ accTemp)).sum
120- // Now mutate each solution and accept/reject
121- // according to the acceptance probability
122-
123- newEnergyLandscape = currentEnergyLandscape.map((config) => {
124- // mutate this config
125- val new_config = mutate(config._2, mutTemp)
126- val new_energy = system.energy(new_config, options)
127- val ans = if (new_energy < config._1) {
128- (new_energy, new_config)
129- } else {
130- val acc = acceptance(new_energy, couplingFactor, accTemp)
131- if (Random .nextDouble <= acc) (new_energy, new_config) else config
132- }
133- ans
134- })
104+ val initialEnergyLandscape = getEnergyLandscape(initialConfig, options, meanFieldPrior)
135105
136- currentEnergyLandscape = newEnergyLandscape
106+ val gamma_init = computeCouplingFactor(initialEnergyLandscape.map(_._1), accTemp)
137107
108+ var acceptanceProbs : List [Double ] = initialEnergyLandscape.map(c => {
109+ computeAcceptanceProb(c._1, c._1, gamma_init, accTemp)
138110 })
139111
112+ val hyp = initialConfig.keys
113+
114+ val usePriorFlag : Boolean = hyp.forall(meanFieldPrior.contains)
115+
116+ def CSATRec (eLandscape : List [(Double , Map [String , Double ])], it : Int ): List [(Double , Map [String , Double ])] =
117+ it match {
118+ case 0 => eLandscape
119+ case num =>
120+ logger.info(" **************************" )
121+ logger.info(" CSA Iteration: " + (MAX_ITERATIONS - it+ 1 ))
122+ // mutate each element of the grid with
123+ // the generating distribution
124+ // and accept using the acceptance distribution
125+ mutTemp = mutationTemperature(iTemp)(it)
126+ accTemp = variant match {
127+ case AbstractCSA .MwVC =>
128+ val (_,variance) = utils.getStats(acceptanceProbs.map(DenseVector (_)))
129+
130+ if (variance(0 ) < sigmaD)
131+ accTemp * (1 - alpha)
132+ else accTemp * (1 + alpha)
133+ case _ =>
134+ acceptanceTemperature(iTemp)(it)
135+ }
136+
137+ val maxEnergy = eLandscape.map(_._1).max
138+
139+ val couplingFactor = computeCouplingFactor(eLandscape.map(t => t._1 - maxEnergy), accTemp)
140+
141+ // Now mutate each solution and accept/reject
142+ // according to the acceptance probability
143+ val (newEnergyLandscape,probabilities) = eLandscape.map((config) => {
144+ // mutate this config
145+ val new_config = mutate(config._2, mutTemp)
146+
147+ val priorEnergy =
148+ if (usePriorFlag)
149+ new_config.foldLeft(0.0 )(
150+ (p_acc, keyValue) => p_acc - meanFieldPrior(keyValue._1).underlyingDist.logPdf(keyValue._2)
151+ )
152+ else 0.0
153+
154+ val new_energy = system.energy(new_config, options) + priorEnergy
155+
156+ logger.info(" New Configuration: \n " + GlobalOptimizer .prettyPrint(new_config))
157+ logger.info(" Energy = " + new_energy)
158+
159+ // Calculate the acceptance probability
160+ val acceptanceProbability =
161+ computeAcceptanceProb(
162+ new_energy - maxEnergy, config._1,
163+ couplingFactor, accTemp)
164+
165+ val ans = if (new_energy < config._1) {
166+
167+ logger.info(" Status: Accepted\n " )
168+ ((new_energy, new_config), acceptanceProbability)
169+
170+ } else {
171+
172+ if (Random .nextDouble <= acceptanceProbability) {
173+ logger.info(" Status: Accepted\n " )
174+ ((new_energy, new_config), acceptanceProbability)
175+ } else {
176+ logger.info(" Status: Rejected\n " )
177+ (config, acceptanceProbability)
178+ }
179+ }
180+
181+ ans
182+ }).unzip
183+
184+ acceptanceProbs = probabilities
185+ CSATRec (newEnergyLandscape, it- 1 )
186+ }
187+
188+ CSATRec (initialEnergyLandscape, MAX_ITERATIONS )
189+ }
140190
141- val optimum = currentEnergyLandscape.keys.min
191+ }
142192
143- logger.info(" Optimum value of energy is: " + optimum+
144- " \n Configuration: " + currentEnergyLandscape(optimum))
193+ object AbstractCSA {
145194
146- system.energy(currentEnergyLandscape(optimum), options)
147- (system, currentEnergyLandscape(optimum))
195+ val MuSA = " CSA-MuSA"
196+ val BA = " CSA-BA"
197+ val M = " CSA-M"
198+ val MwVC = " CSA-MwVC"
199+ val SA = " SA"
200+
201+ def algorithm (variant : String ): String = variant match {
202+ case MuSA => " Multi-state Simulated Annealing"
203+ case BA => " Blind Acceptance"
204+ case M => " Modified CSA"
205+ case MwVC => " Modified CSA with Variance Control"
148206 }
149- }
150207
151- object AbstractCSA {
152- def apply [M <: GloballyOptimizable ](model : M ,
153- initialConfig : Map [String , Double ],
154- options : Map [String , String ] = Map ()): M = {
155- new AbstractCSA [M ](model).optimize(initialConfig, options)._1
208+ def couplingFactor (variant : String )(
209+ landscape : Seq [Double ],
210+ Tacc : Double ): Double = {
211+
212+ if (variant == MuSA || variant == BA )
213+ landscape.map(energy => math.exp(- 1.0 * energy/ Tacc )).sum
214+ else if (variant == M || variant == MwVC )
215+ landscape.map(energy => math.exp(energy/ Tacc )).sum
216+ else 1.0
217+
218+ }
219+
220+ def acceptanceProbability (variant : String )(
221+ energy : Double , oldEnergy : Double ,
222+ gamma : Double , temperature : Double ) = {
223+
224+ if (variant == MuSA )
225+ math.exp(- 1.0 * energy/ temperature)/ (math.exp(- 1.0 * energy/ temperature)+ gamma)
226+ else if (variant == BA )
227+ 1.0 - (math.exp(- 1.0 * oldEnergy/ temperature)/ gamma)
228+ else if (variant == M || variant == MwVC )
229+ math.exp(oldEnergy/ temperature)/ gamma
230+ else gamma/ (1.0 + math.exp((energy - oldEnergy)/ temperature))
231+
232+ }
233+
234+ def varianceDesired (variant : String )(m : Int ): Double = {
235+ if (variant == MuSA || variant == BA )
236+ 0.99
237+ else
238+ 0.99 * (m- 1 )/ math.pow(m, 2.0 )
239+
156240 }
241+
157242}
158243
244+
0 commit comments