import sys import os import time from datetime import datetime from collections import namedtuple, defaultdict import logging from Queue import Empty import math from random import randrange, uniform, choice, random from pygene.gene import FloatGeneMax, FloatGeneExchange from pygene.organism import MendelOrganism from pygene.population import Population import numpy as np from runphoton import BBox, Resolution from utils.materials import loadMaterials, materialsAsXml, Brdf from utils.xml import FileToTree, prettyPrint, StringToTree, subelement, TreeToFile from utils.path import path from utils import daemonthread from utils import accelerated import simulate as simulator log = logging.getLogger('lws') DEVICE_ADAPTION = [-90, -85, -80, -75,-70, -65, -60, -55, -50, -45, -40,-30, -20, 0] DEVICE_ADAPTION = np.arange(-95, -20, 3.0) #~ DEVICE_ADAPTION = [-90, -80, -70, -60, -50, -40, -30, 0] DA_GENE_NAMES = ['da_%s' % -i for i in DEVICE_ADAPTION] # set via server.py - holds lws instance app = None class PowerGene(FloatGeneExchange): """ Gene which represents the numbers used in our organism """ # genes get randomly generated within this range # initialized in startOptimizeRun() randMin = 0 randMax = 0 # probability of mutation mutProb = 0.2 # degree of mutation mutAmt = 0.2 def mutate(self): """ Mutate this gene's value by a random amount within the range, which is determined by multiplying self.mutAmt by the distance of the gene's current value from either endpoint of legal values perform mutation IN-PLACE, ie don't return mutated copy """ if random() < 0.5: # mutate downwards: fact = 1.0 - uniform(0, self.mutAmt) newval = self.value * fact self.value = max(self.randMin, newval) else: # mutate upwards: fact = 1.0 + uniform(0, self.mutAmt) newval = self.value * fact self.value = min(self.randMax, newval) def randomValue(self): """ Generates a plausible random value for this gene. """ values = [] v = self.randMin while v < self.randMax: values.append(v) v = v * 1.3 return choice(values) class NormalizerGene(FloatGeneExchange): """ Gene which represents the numbers used in our organism """ # genes get randomly generated within this range randMin = 89.0 randMax = 91.0 # probability of mutation mutProb = 0.2 # degree of mutation mutAmt = 0.2 class DeviceAdaptionGene(FloatGeneExchange): # genes get randomly generated within this range randMin = -8.0 randMax = 8.0 mutProb = 0.1 mutAmt = 0.1 class MaterialGene(FloatGeneExchange): """ Gene which represents the numbers used in our organism """ # genes get randomly generated within this range randMin = 0.0 randMax = 1.0 # probability of mutation mutProb = 0.2 # degree of mutation mutAmt = 0.1 def optimizeDeviceAdaption(station2apid_locid2measurement, apid_locid2estimated, max_generations): stations = station2apid_locid2measurement.keys() shape = station2apid_locid2measurement.values()[0].shape class DeviceAdaptionConverger(MendelOrganism): """ Implements the organism which tries to converge a function """ genome = {'%s_%s' % (station, n): DeviceAdaptionGene for n in DA_GENE_NAMES for station in stations} def genes2values(self): da_values = {} for station in station2apid_locid2measurement: da_values[station] = [(DEVICE_ADAPTION[i], self['%s_%s' % (station, n)]) for i, n in enumerate(DA_GENE_NAMES)] return da_values def fitness(self): if self.fitness_cache is not None: return self.fitness_cache da_values = self.genes2values() self.used_values = da_values avg_delta = accelerated.compareForDeviceAdaption(da_values, estimated_3d, station2apid_locid2measurement) #~ deltas = [] #~ for i in range(shape[0]): #~ for j in range(shape[1]): #~ res = a[i, j, 0] - self.apid_locid2measurement[i, j] #~ if np.isfinite(res): #~ deltas.append(abs(res)) #~ avg_delta = sum(deltas) / len(deltas) fitness = avg_delta self.fitness_cache = fitness return fitness #~ print time.time() - t def __repr__(self): s = [] for station in ['iconia']: davalues = [(DEVICE_ADAPTION[i], self['%s_%s' % (station, n)]) for i, n in enumerate(DA_GENE_NAMES)] s.append(station + ': ' + ', '.join('%s: %.1f' % e for e in davalues)) return " avg-delta=%.2f\n%s" % (self.fitness_cache, '\n'.join(s)) ## ====================== # reshape to 3 dim array estimated_3d = apid_locid2estimated.reshape(shape[0], shape[1], 1).astype(np.float32) childcull = 50 childcount = 100 pop = Population(species=DeviceAdaptionConverger, init=300) t = time.time() #~ print DeviceAdaptionConverger().fitness() gen_idx = 0 while gen_idx < max_generations: # execute a generation pop.gen(childcull, childcount) log.info('generation of genes: %s [population: %s] [fitness: %.2f]' % (gen_idx, len(pop.organisms), pop.best().fitness_cache)) #~ if gen_idx-1 % 25 == 0: #~ log.info('%r' % pop.best()) gen_idx += 1 best = pop.best() return best.used_values, best.fitness() class RadioPropConverger(MendelOrganism): """ Implements the organism which tries to converge a function """ genome = dict() ORGANISM_COUNT = 0 def __init__(self, *args, **kwargs): RadioPropConverger.ORGANISM_COUNT += 1 self.organismcount = RadioPropConverger.ORGANISM_COUNT # our id MendelOrganism.__init__(self, *args, **kwargs) def prepare_fitness(self): if not app.optimizer.running: return if not getattr(self, 'prepared', False): self.runs = set() materials = {} for matname in app.optimizer.adjustable_materials: materials[matname] = Brdf(reflect=self['%s_reflect' % matname], alpha=self['%s_alpha' % matname]) queued_ids = [] # only for logging for apid in app.optimizer.activeAPs: apcfg = app.config['aps'][apid] power = self['%s_power' % apcfg['powerid']] ap = simulator.AccessPoint(id=apid, x=apcfg['x'], y=apcfg['y'], z=apcfg['z'], power=power) density = app.optimizer.density bbox = app.optimizer.bbox resolution = app.optimizer.resolution scenename = app.optimizer.scene numphotons = app.optimizer.numphotons disabledObjects = app.optimizer.disabledObjects run = simulator.Run(scenename, app.activeScene['objfile'], materials, ap, bbox, resolution, density, numphotons, code='inspect_locids', locidinfo=app.known_locations, disabledObjects=disabledObjects) # if we use self.organismcount as priority, the organisms are evaluated pretty much sequentially # and that is good for parallelizablity priority = self.organismcount app.simulator.queue(run, priority=priority) self.runs.add(run) queued_ids.append(run.id) log.info('organism %s, queued %s runs: %s' % (self.organismcount, len(queued_ids), ','.join(queued_ids))) self.prepared = True def fitness_complete(self): remaining = len(self.runs - app.simulator.finished_jobs) return remaining == 0 def fitness(self): if self.fitness_cache is not None: return self.fitness_cache self.prepare_fitness() while True: if not app.optimizer.running: #~ log.info('optimizer down, skipping fitness calculation') return 10**9 # return something really unfit if self.fitness_complete(): break else: waiting, running, finished = app.simulator.queryRuns(simulator.Run.known) remaining = len(self.runs - finished) log.info('organism %s fitness calc: %s remaining [waiting: %s, running %s]' % (self.organismcount, remaining, len(waiting), len(running))) time.sleep(3.0) #~ da_values = {} #~ for station in app.optimizer.stations: #~ da_values[station] = [(DEVICE_ADAPTION[i], self['%s_%s' % (station, n)]) for i, n in enumerate(DA_GENE_NAMES)] #~ rssieq = self['rssi_eq'] max_locid = max(app.known_locations.keys()) apid_locid2estimated = np.zeros(shape=(len(app.optimizer.activeAPs), max_locid+1)) + np.inf for run in self.runs: apidx = app.optimizer.activeAPs.index(run.ap.id) for locid in app.known_locations.keys(): estimated = run.locid2rssi[str(locid)] #~ print estimated apid_locid2estimated[apidx, locid] = estimated OPTIMIZE_DAV = False if OPTIMIZE_DAV: da_values, avg_delta_opt = optimizeDeviceAdaption(app.optimizer.station2measurements, apid_locid2estimated, 20) else: da_values = {} for station in app.optimizer.station2measurements: da_values[station] = [] avg_delta_opt = None deltas = [] values = defaultdict(list) # locid, measurements, estimated, delta for each ap for run in self.runs: apidx = app.optimizer.activeAPs.index(run.ap.id) for station in da_values: apid_locid2measurement = app.optimizer.station2measurements[station] for locid in app.known_locations.keys(): measured = apid_locid2measurement[apidx, locid] # read simulated rss value for locid # that has been transfered from worker node via worker_code_inspect_locids.py estimated_raw = run.locid2rssi[str(locid)] if estimated_raw == 0: continue estimated = math.log(estimated_raw, 10) * 10 normalized = accelerated.normalize(estimated_raw, da_values[station], isInDB=False) delta = abs(measured - normalized) if np.isnan(delta) or not np.isfinite(delta) or delta > 100: continue if not normalized > -100: continue deltas.append(delta) values[run.ap.id].append((station, locid, measured, estimated, normalized, delta)) # all runs have the same material params # use the params of the last run materials = run.materials #~ res = run(x) avg_delta = np.sum(deltas) / float(len(deltas)) if avg_delta_opt is not None: assert round(avg_delta_opt, 4) == round(avg_delta, 4), '%s != %s' % (avg_delta_opt, avg_delta) powers = [(g[:-6], self[g]) for g in self.genome.keys() if g.endswith('_power')] app.optimizer.newSimulationResult(self.organismcount, materials, powers, da_values, avg_delta, values) fitness = avg_delta self.fitness_cache = fitness return fitness def __repr__(self): return "organism %s avg-delta=%.2f" % (self.organismcount, self.fitness_cache if self.fitness_cache is not None else -1.0) #~ Iteration = namedtuple('Iteration', 'iter avg_delta paramxml path') #~ def loadData(REP_DIR): #~ dirs = os.listdir(REP_DIR) #~ def parse(name): #~ try: #~ return datetime.strptime(name, '%Y-%m-%dT%H_%M_%S') #~ except Exception: #~ return None #~ optruns = [] #~ for name in dirs: #~ if parse(name) is None: #~ continue #~ subdirs = [os.path.join(REP_DIR, name, p) for p in os.listdir(os.path.join(REP_DIR, name))] #~ if len(subdirs) == 1: #~ session = os.path.split(subdirs[0])[-1] #~ bestdir = os.path.join(subdirs[0], 'best') #~ best_iterations = [] #~ avg_delta = None #~ iter = 0 #~ if os.path.exists(bestdir): #~ for n in os.listdir(bestdir): #~ iter = int(n.split('_')[0]) #~ avg_delta = float(n.split('_')[1]) #~ paramxml = open(os.path.join(bestdir, n, 'params.xml')).read() #~ best_iterations.append(Iteration(iter, avg_delta, paramxml, os.path.join(bestdir, n))) #~ lessThanOneDayOld = os.stat(subdirs[0]).st_mtime > time.time() - 3600 * 24 * 7 #~ interesting = (avg_delta is not None and avg_delta < 4.0) or (lessThanOneDayOld and iter > 1) #~ if interesting: #~ optruns.append((name, parse(name), session, best_iterations)) #~ # return sorted by ts ascending #~ return list(sorted(optruns, key=lambda x: x[1])) class Optimizer(object): def __init__(self): self.runner = None self.running = False self.startts = None self.timelimit = 0 self.genlimit = 0 self.sessiondir = None # directory for the current optsession with count: optsession_12345 self.optsession = None # current optsession self.station = None # list of stations in string form self.stations = None # list of stations that should be included in optimize run self.gen_idx = 0 # current generation index self.min_avg_delta = 10**9 # best avg delta of current optsession (== best fitness) # optsession init params self.power = self.density = self.numphotons = self.resolution = self.bbox = self.scene = None self.tmpdir = app.config['tmp_optruns'] if not self.tmpdir.exists(): self.tmpdir.mkdir() # loaded measurements for each active station self.station2measurements = {} self.activeAPs = [] daemonthread(target=self._checkTimelimit, name="optimizer_timelimit_check") def _checkTimelimit(self): while True: try: if self.running and self.genlimit > 0 and self.gen_idx > self.genlimit: log.warn('GENERATION LIMIT reached, stopping!') self.running = False elif (self.running and self.timelimit > 0 and time.time() - self.startts > self.timelimit): log.warn('TIMELIMIT reached, stopping!') self.running = False except Exception: log.exception('error during _checkTimelimit') try: time.sleep(1) except Exception: pass # ignore error during shutdown cleanup #~ if self.startts is not None: #~ log.info('TIMELIMIT: %s seconds remaining' % (self.timelimit - (time.time() - self.startts))) def loadPopulation(self, fname): tree = FileToTree(fname) #~ print 'load', fname print prettyPrint(tree) #~ tree.xpath() def newSimulationResult(self, orgid, materials, powers, da_values, avg_delta, deltavalues): minsfile = self.sessiondir.joinpath('mins.txt') minsfile.touch() allfile = self.sessiondir.joinpath('population_%s.txt' % self.gen_idx) allfile.touch() dt = datetime.now().strftime("%Y-%m-%dT%H:%M:%S") def storeToFile(f): s = '%s gen:%s org:%s avg-delta:%.2f ' % (dt, self.gen_idx, orgid, avg_delta) for matname, brdf in sorted(materials.items()): s += '%s:%.3f,%.3f ' % (matname, brdf.reflect, brdf.alpha) for powid, power in powers: s += '%s:%.3e ' % (powid, power) #~ s += 'n/rssieq:%.3f ' % rssi_eq for station, dav in sorted(da_values.items()): for da_rssi, v in dav: s += 'da/%s_%s:%.1f ' % (station, da_rssi, v) open(f, 'a').write(s + '\n') storeToFile(allfile) if avg_delta < self.min_avg_delta: log.info('new minimum reached: %.2f' % avg_delta) storeToFile(minsfile) self.min_avg_delta = avg_delta minsdir = self.sessiondir.joinpath('mins') if not minsdir.exists(): minsdir.mkdir() deltafile = minsdir.joinpath('deltas_%s_%s.txt' % (self.gen_idx, orgid)) with deltafile.open('w') as f: for apid, values in sorted(deltavalues.items()): f.write('ap: %s\n' % apid) f.write('locid measured simulated normalized delta\n') for station, locid, measured, estimated, normalized, delta in values: f.write('%s %s %.1f %.1f %.1f %.1f\n' % (station, locid, measured, estimated, normalized, delta)) def startOptimizeRun(self, optsession, station, power, density, numphotons, resolution, bbox, scene, startpop, childcount, childcull, mutprob, mutamt, resume=None, neworgs=0, timelimit=0, genlimit=0): if self.running: log.error('an optimizer run is still running') pass else: self.station = station self.stations = station.split(',') self.timelimit = int(timelimit) self.genlimit = int(genlimit) self.power = power self.density = density self.numphotons = numphotons self.resolution = resolution self.bbox = bbox self.scene = scene self.startpop = startpop self.childcount = childcount self.childcull = childcull # setup PowerGene PowerGene.randMin = power * 0.05 # 1/10 of initial power = lower bound PowerGene.randMax = power * 50.0 # 10 times initial power = upper bound PowerGene.mutProb = mutprob PowerGene.mutAmt = mutamt MaterialGene.mutProb = mutprob MaterialGene.mutAmt = mutamt sesscount = max([0] + [int(d.name.split('_')[1]) for d in self.tmpdir.dirs() if d.name.split('_')[0] == optsession]) self.optsession = '%s_%s' % (optsession, sesscount+1) self.sessiondir = self.tmpdir.joinpath(self.optsession) self.sessiondir.mkdir() ## HACKISH for eval self.adjustable_materials = ('MatConcrete', 'MatLightWalls') self.adjustable_materials = ('MatConcrete',) self.adjustable_materials = app.activeScene['materials'].scalars # [so.name for so in self.scene_objects if not ('light' in so.name.lower() or 'concrete' in so.name.lower())] self.disabledObjects = ['Iron_Doors_ID1', 'Fascade_OG2_M6.002', 'Cupboards_Plane.002', 'Fascade_EG_M6.001', 'Glass_Doors_Plane', 'Fascade_OG1_M6', 'Railing_Plane.Railing', 'Glass_Windows_OG1_GW01', 'Normal_Doors_ND1', 'Glass_Windows_OG2_GW01.002', 'Glass_Windows_EG1_GW01.001', 'Stairs_Plane.Stairs', 'Tables_Plane.004', 'Hardware_Plane.001'] self.disabledObjects = [] # see also simulator remapMaterials initfile = self.sessiondir.joinpath('init.txt') s = 'session: %s\n' % optsession s += 'timelimit: %s\n' % timelimit s += 'genlimit: %s\n' % genlimit s += 'station: %s\n' % station s += 'power: %.1e\n' % power s += 'density: %.3f\n' % density s += 'numphotons: %s\n' % numphotons s += 'resolution: %s,%s,%s\n' % resolution s += 'bbox: %s,%s,%s,%s,%s,%s\n' % bbox s += 'scene: %s\n' % scene s += 'startpop: %s\n' % startpop s += 'childcount: %s\n' % childcount s += 'childcull: %s\n' % childcull s += 'mutprob: %s\n' % mutprob s += 'mutamt: %s\n' % mutamt initfile.write_text(s) apids = [apid for apid in app.config['aps'] if app.config['aps'][apid]['optimize']] self.activeAPs[:] = apids self._setupGenes() if resume: d = self.tmpdir.joinpath(resume) if not d.exists(): raise RuntimeError('no session %s available' % resume) popfiles = d.files('population_*.xml') if len(popfiles) == 0: raise RuntimeError('no population for session %s available' % resume) popfiles.sort(key=lambda f: f.mtime) organisms = self.xml2population(FileToTree(popfiles[-1])) for i in range(neworgs): organisms.append(RadioPropConverger()) else: organisms = None self.min_avg_delta = 10**9 self.running = True self.startts = time.time() self.runner = daemonthread(target=self._optimizeRun, name="optimizeRun", args=(organisms, )) def stopOptimizeRun(self): if self.running: self.running = False app.simulator.finishAll() def xml2population(self, tree): self._setupGenes() organisms = [] for orgEl in tree.xpath('//organism'): genes = {} for geneEl in orgEl.xpath('genepair'): # get class from current module cls = globals()[geneEl.attrib['class']] g1 = float(geneEl.attrib['g1']) g2 = float(geneEl.attrib['g2']) genes[geneEl.attrib['name']] = (cls(g1), cls(g2)) organism = RadioPropConverger(**genes) organisms.append(organism) return organisms def population2xml(self, pop, **kwds): tree = StringToTree('') for k, v in kwds.items(): tree.getroot().attrib[k] = str(v) unique_genes = set() for org in pop.organisms: orgEl = subelement(tree.getroot(), 'organism', attribs={'fitness': org.fitness(), 'count': org.organismcount}) for name, cls in org.genome.items(): pairEl = subelement(orgEl, "genepair", attribs={'name': name, 'class': cls.__name__,}) pair = org.genes[name] pairEl.attrib['g1'] = str(pair[0].value) pairEl.attrib['g2'] = str(pair[1].value) unique_genes.add(cls) for cls in unique_genes: attribs = {'class': cls.__name__, 'mutamt': cls.mutAmt, 'mutprob': cls.mutProb, 'min': cls.randMin, 'max': cls.randMax} subelement(tree.getroot(), 'genespec', attribs=attribs) return tree def _setupGenes(self): aps = app.config['aps'] powids = {aps[apid]['powerid'] for apid in self.activeAPs} RadioPropConverger.genome.clear() RadioPropConverger.genome.update(dict( [('%s_reflect' % matname, MaterialGene) for matname in self.adjustable_materials] + [('%s_alpha' % matname, MaterialGene) for matname in self.adjustable_materials] + [('%s_power' % powid, PowerGene) for powid in powids] )) def _loadMeasurements(self): self.station2measurements.clear() max_locid = max(app.known_locations.keys()) for station in self.stations: # initialize with infinity a = self.station2measurements[station] = np.zeros(shape=(len(self.activeAPs), max_locid+1)) + np.inf locid2measurecount = defaultdict(int) remainingActive = [] for i, apid in enumerate(self.activeAPs): locid2measurement = app.getMeasurements(station, apid) if locid2measurement is not None: #~ locid2measurement_overlay = app.getMeasurements(station, apid, fromMeasurementStore=False) for locid, (x, y, z, measured, all_rssis) in locid2measurement.items(): a[i, locid] = float(measured) locid2measurecount[locid] += 1 # overwrite with values from tracked paths #~ if locid in locid2measurement_overlay: #~ x, y, z, measured, all_rssis = locid2measurement_overlay[locid] #~ a[i, locid] = float(measured) remainingActive.append(apid) ### reduce to set of location with the N = 50 sorted_locs = sorted(locid2measurecount.items(), key=lambda e: e[1]) for l, c in sorted_locs: print l, c topNlocids = set([locid for locid, count in sorted_locs][-N:]) for i, apid in enumerate(self.activeAPs): locid2measurement = app.getMeasurements(station, apid) if locid2measurement is not None: for locid, (x, y, z, measured, all_rssis) in locid2measurement.items(): if locid not in topNlocids: # clear measurement by setting to inf a[i, locid] = np.inf log.info('%s of %s remaining active aps due to available measurements' % (len(remainingActive), i)) self.activeAPs[:] = remainingActive def _optimizeRun(self, startorganisms): try: self._loadMeasurements() if startorganisms is not None: pop = Population(*startorganisms, init=self.startpop) else: pop = Population(species=RadioPropConverger, init=self.startpop) for org in pop.organisms: org.prepare_fitness() self.gen_idx = 0 while self.running: log.info('generation of genes: %s [population: %s]' % (self.gen_idx, len(pop.organisms))) # execute a generation pop.gen(self.childcull, self.childcount) f = self.tmpdir.joinpath(self.sessiondir, 'population_%s.xml' % self.gen_idx) tree = self.population2xml(pop, init=self.startpop, childCount=self.childcount, childCull=self.childcull) TreeToFile(tree, f) self.gen_idx += 1 except Exception: log.exception('error during optrun') finally: self.running = False def getDeviceAdaptionFromOptrun(self, optsession, device): f = self.tmpdir.joinpath(optsession, 'mins.txt') davalues = defaultdict(list) if not f.exists(): log.error('got no mins for optrun %s - empty davalues' % optsession) return davalues # take last minimum mindata = f.text().strip().split('\n')[-1].split() for s in mindata[1:]: name, value = s.split(':') if name.startswith('da/'): device_darssi = name[3:].split('_') if len(device_darssi) == 1: # backward compat davalues[device].append((float(device_darssi[0]), float(value))) else: dev, darssi = device_darssi davalues[dev].append((float(darssi), float(value))) return davalues