source: branches/2929_PrioritizedGrammarEnumeration/HeuristicLab.Algorithms.DataAnalysis.PGE/3.3/src/go-pge/pge/pge_search.go @ 16620

Last change on this file since 16620 was 16620, checked in by hmaislin, 9 months ago

#2929: Reorganized folder structure for make script, removed explicit marshalling, erased go-side logging

File size: 15.2 KB
Line 
1package pge
2
3import (
4  "fmt"
5  "io/ioutil"
6  "math"
7  "strconv"
8  "strings"
9
10  levmar "go-levmar"
11  config "go-pge/config"
12  probs "go-pge/problems"
13  expr "go-symexpr"
14)
15
16type PgeConfig struct {
17  // search params
18  maxGen        int
19  pgeRptEpoch   int
20  pgeRptCount   int
21  pgeArchiveCap int
22
23  simprules expr.SimpRules
24  treecfg   *probs.TreeParams
25
26  // PGE specific options
27  peelCnt     int
28  sortType    probs.SortType
29  zeroEpsilon float64
30
31  initMethod string
32  growMethod string
33
34  evalrCount int
35}
36
37func pgeConfigParser(field, value string, config interface{}) (err error) {
38
39  PC := config.(*PgeConfig)
40
41  switch strings.ToUpper(field) {
42  case "MAXGEN":
43    PC.maxGen, err = strconv.Atoi(value)
44  case "PGERPTEPOCH":
45    PC.pgeRptEpoch, err = strconv.Atoi(value)
46  case "PGERPTCOUNT":
47    PC.pgeRptCount, err = strconv.Atoi(value)
48  case "PGEARCHIVECAP":
49    PC.pgeArchiveCap, err = strconv.Atoi(value)
50
51  case "PEELCOUNT":
52    PC.peelCnt, err = strconv.Atoi(value)
53
54  case "EVALRCOUNT":
55    PC.evalrCount, err = strconv.Atoi(value)
56
57  case "SORTTYPE":
58    switch strings.ToLower(value) {
59    case "paretotrainerror":
60      PC.sortType = probs.PESORT_PARETO_TRN_ERR
61    case "paretotesterror":
62      PC.sortType = probs.PESORT_PARETO_TST_ERR
63
64    default:
65    }
66
67  case "ZEROEPSILON":
68    PC.zeroEpsilon, err = strconv.ParseFloat(value, 64)
69
70  default:
71    // check augillary parsable structures [only TreeParams for now]
72    if PC.treecfg == nil {
73      PC.treecfg = new(probs.TreeParams)
74    }
75    found, ferr := probs.ParseTreeParams(field, value, PC.treecfg)
76    if ferr != nil {
77      return ferr
78    }
79    if !found {
80    }
81
82  }
83  return
84}
85
86type PgeSearch struct {
87  id   int
88  cnfg PgeConfig
89  prob *probs.ExprProblem
90  iter int
91  stop bool
92
93  // comm up
94  commup *probs.ExprProblemComm
95
96  // comm down
97
98  // best exprs
99  Best *probs.ReportQueue
100
101  // training data in C format
102  c_input  []levmar.C_double
103  c_ygiven []levmar.C_double
104
105  // logs
106
107  // equations visited
108  Trie  *IpreNode
109  Queue *probs.ReportQueue
110
111  // eval channels
112  eval_in  chan expr.Expr
113  eval_out chan *probs.ExprReport
114
115  // genStuff
116  GenRoots   []expr.Expr
117  GenLeafs   []expr.Expr
118  GenNodes   []expr.Expr
119  GenNonTrig []expr.Expr
120
121  // FFXish stuff
122  ffxBases []expr.Expr
123
124  // statistics
125  neqns    int
126  ipreCnt  int
127  maxSize  int
128  maxScore int
129  minError float64
130}
131
132func (PS *PgeSearch) GetMaxIter() int {
133  return PS.cnfg.maxGen
134}
135func (PS *PgeSearch) GetPeelCount() int {
136  return PS.cnfg.peelCnt
137}
138func (PS *PgeSearch) SetMaxIter(iter int) {
139  PS.cnfg.maxGen = iter
140}
141func (PS *PgeSearch) SetPeelCount(cnt int) {
142  PS.cnfg.peelCnt = cnt
143}
144func (PS *PgeSearch) SetInitMethod(init string) {
145  PS.cnfg.initMethod = init
146}
147func (PS *PgeSearch) SetGrowMethod(grow string) {
148  PS.cnfg.growMethod = grow
149}
150func (PS *PgeSearch) SetEvalrCount(cnt int) {
151  PS.cnfg.evalrCount = cnt
152}
153
154func (PS *PgeSearch) ParseConfig(filename string) {
155  fmt.Printf("Parsing PGE Config: %s\n", filename)
156  data, err := ioutil.ReadFile(filename)
157  if err != nil {
158  }
159  err = config.ParseConfig(data, pgeConfigParser, &PS.cnfg)
160  if err != nil {
161  }
162}
163
164func (PS *PgeSearch) SetSort(typen int) {
165  if typen == 1 {
166    PS.Best.SetSort(probs.GPSORT_PARETO_TST_ERR)
167    PS.Best.SetSort(probs.GPSORT_PARETO_TST_ERR)
168  } else {
169    PS.Best.SetSort(probs.PESORT_PARETO_TST_ERR)
170    PS.Best.SetSort(probs.PESORT_PARETO_TST_ERR)
171  }
172}
173
174func (PS *PgeSearch) Init(done chan int, prob *probs.ExprProblem, logdir string, input interface{}) {
175
176  fmt.Printf("Init'n PGE\n")
177  // setup data
178
179  // open logs
180  //PS.initLogs(logdir)
181
182  PS.stop = false
183
184  // copy in common config options
185  PS.prob = prob
186  if PS.cnfg.treecfg == nil {
187    PS.cnfg.treecfg = PS.prob.TreeCfg.Clone()
188  }
189  srules := expr.DefaultRules()
190  srules.ConvertConsts = true
191  PS.cnfg.simprules = srules
192
193  fmt.Println("Roots:   ", PS.cnfg.treecfg.RootsS)
194  fmt.Println("Nodes:   ", PS.cnfg.treecfg.NodesS)
195  fmt.Println("Leafs:   ", PS.cnfg.treecfg.LeafsS)
196  fmt.Println("NonTrig: ", PS.cnfg.treecfg.NonTrigS)
197
198  PS.GenRoots = make([]expr.Expr, len(PS.cnfg.treecfg.Roots))
199  for i := 0; i < len(PS.GenRoots); i++ {
200    PS.GenRoots[i] = PS.cnfg.treecfg.Roots[i].Clone()
201  }
202  PS.GenNodes = make([]expr.Expr, len(PS.cnfg.treecfg.Nodes))
203  for i := 0; i < len(PS.GenNodes); i++ {
204    PS.GenNodes[i] = PS.cnfg.treecfg.Nodes[i].Clone()
205  }
206  PS.GenNonTrig = make([]expr.Expr, len(PS.cnfg.treecfg.NonTrig))
207  for i := 0; i < len(PS.GenNonTrig); i++ {
208    PS.GenNonTrig[i] = PS.cnfg.treecfg.NonTrig[i].Clone()
209  }
210
211  PS.GenLeafs = make([]expr.Expr, 0)
212  for _, t := range PS.cnfg.treecfg.LeafsT {
213    switch t {
214    case expr.TIME:
215      PS.GenLeafs = append(PS.GenLeafs, expr.NewTime())
216
217    case expr.VAR:
218      fmt.Println("Use Vars: ", PS.cnfg.treecfg.UsableVars)
219      for _, i := range PS.cnfg.treecfg.UsableVars {
220        PS.GenLeafs = append(PS.GenLeafs, expr.NewVar(i))
221      }
222
223    case expr.SYSTEM:
224      for i := 0; i < PS.prob.Train[0].NumSys(); i++ {
225        PS.GenLeafs = append(PS.GenLeafs, expr.NewSystem(i))
226      }
227
228    }
229  }
230
231  /*** FIX ME
232  PS.GenLeafs = make([]expr.Expr, len(PS.cnfg.treecfg.Leafs))
233  for i := 0; i < len(PS.GenLeafs); i++ {
234    PS.GenLeafs[i] = PS.cnfg.treecfg.Leafs[i].Clone()
235  }
236  ***/
237
238  //fmt.Println("Roots:   ", PS.GenRoots)
239  //fmt.Println("Nodes:   ", PS.GenNodes)
240  //fmt.Println("Leafs:   ", PS.GenLeafs)
241  //fmt.Println("NonTrig: ", PS.GenNonTrig)
242
243  // setup communication struct
244  PS.commup = input.(*probs.ExprProblemComm)
245
246  // initialize bbq
247  PS.Trie = new(IpreNode)
248  PS.Trie.val = -1
249  PS.Trie.next = make(map[int]*IpreNode)
250
251  PS.Best = probs.NewReportQueue()
252  PS.Best.SetSort(probs.GPSORT_PARETO_TST_ERR)
253
254  PS.Queue = PS.GenInitExpr()
255  PS.Queue.SetSort(probs.PESORT_PARETO_TST_ERR)
256
257  PS.neqns = PS.Queue.Len()
258
259  PS.minError = math.Inf(1)
260
261  PS.eval_in = make(chan expr.Expr, 4048)
262  PS.eval_out = make(chan *probs.ExprReport, 4048)
263
264  fmt.Printf("Init5a\n")
265  for i := 0; i < PS.cnfg.evalrCount; i++ {
266    go PS.Evaluate()
267  }
268}
269
270func (PS *PgeSearch) Stop() {
271  PS.stop = true
272}
273
274func (PS *PgeSearch) Evaluate() {
275
276  for !PS.stop {
277    e := <-PS.eval_in
278    if e == nil {
279      continue
280    }
281    //re :=
282    //fmt.Printf("reg: %v\n", re)
283    PS.eval_out <- RegressExpr(e, PS.prob)
284  }
285
286}
287
288func (PS *PgeSearch) Run() {
289  fmt.Printf("Running PGE\n")
290
291  PS.Loop()
292
293  fmt.Println("PGE exitting")
294
295  PS.Clean()
296  PS.commup.Cmds <- -1
297}
298
299func (PS *PgeSearch) Loop() {
300
301  PS.checkMessages()
302  for !PS.stop {
303
304    fmt.Println("in: PS.step() ", PS.iter)
305    PS.Step()
306
307    // if PS.iter%PS.cnfg.pgeRptEpoch == 0 {
308    PS.reportExpr()
309    // }
310
311    // report current iteration
312    PS.commup.Gen <- [2]int{PS.id, PS.iter}
313    PS.iter++
314
315    PS.Clean()
316
317    PS.checkMessages()
318
319  }
320
321  // done expanding, pull the rest of the regressed solutions from the queue
322  p := 0
323  for PS.Queue.Len() > 0 {
324    e := PS.Queue.Pop().(*probs.ExprReport)
325
326    bPush := true
327    if len(e.Coeff()) == 1 && math.Abs(e.Coeff()[0]) < PS.cnfg.zeroEpsilon {
328      // fmt.Println("No Best Push")
329      bPush = false
330    }
331
332    if bPush {
333      // fmt.Printf("pop/push(%d,%d): %v\n", p, PS.Best.Len(), e.Expr())
334      PS.Best.Push(e)
335      p++
336    }
337
338    if e.TestScore() > PS.maxScore {
339      PS.maxScore = e.TestScore()
340    }
341    if e.TestError() < PS.minError {
342      PS.minError = e.TestError()
343      fmt.Printf("EXITING New Min Error:  %v\n", e)
344    }
345    if e.Size() > PS.maxSize {
346      PS.maxSize = e.Size()
347    }
348  }
349
350  fmt.Println("PGE sending last report")
351  PS.reportExpr()
352
353}
354
355type PeelResult struct {
356  Es *probs.ExprReport
357
358  Nobestpush    bool
359  BestNewMinErr bool
360
361  Bestlen1 int
362  Bestlen2 int
363
364  Coeff     []float64
365  TestScore int
366
367  Expre   expr.Expr
368  ExpreRe *probs.ExprReport
369}
370
371func (PS *PgeSearch) Step() int {
372
373  loop := 0
374  eval_cnt := 0 // for channeled eval
375
376  es := PS.peel()
377
378  ex := PS.expandPeeled(es)
379
380  cnt_ins := 0
381
382  inserts := make(probs.ExprReportArray, 0)
383  for cnt := range ex {
384    E := ex[cnt]
385
386    if E == nil {
387      continue
388    }
389
390    for _, e := range E {
391      if e == nil {
392        continue
393      }
394      if !PS.cnfg.treecfg.CheckExpr(e) {
395        continue
396      }
397
398      // check ipre_trie
399      serial := make([]int, 0, 64)
400      serial = e.Serial(serial)
401      ins := PS.Trie.InsertSerial(serial)
402      if !ins {
403        continue
404      }
405
406      // for serial eval
407      //re := RegressExpr(e, PS.prob) //needed for TestScore calc via RegressExpr
408      //inserts = append(inserts, re)
409
410      // start channeled eval
411      PS.eval_in <- e
412      eval_cnt++
413    }
414  }
415
416  for i := 0; i < eval_cnt; i++ {
417    re := <-PS.eval_out
418    // end channeled eval
419
420    // check for NaN/Inf in re.error  and  if so, skip
421    if math.IsNaN(re.TestError()) || math.IsInf(re.TestError(), 0) {
422      // fmt.Printf("Bad Error\n%v\n", re)
423      continue
424    }
425
426    if re.TestError() < PS.minError {
427      PS.minError = re.TestError()
428    }
429
430    // check for coeff == 0
431    doIns := true
432    for _, c := range re.Coeff() {
433      // i > 0 for free coeff
434      if math.Abs(c) < PS.cnfg.zeroEpsilon {
435        doIns = false
436        break
437      }
438    }
439
440    //fmt.Printf("StepQueue.Push(): %v\n", re)
441    //fmt.Printf("StepQueue.Push(): %v\n", re.Expr())
442
443    if doIns {
444      re.SetProcID(PS.id)
445      re.SetIterID(PS.iter)
446      re.SetUnitID(loop)
447      re.SetUniqID(PS.neqns)
448      loop++
449      PS.neqns++
450      // fmt.Printf("Queue.Push(): %v\n%v\n\n", re.Expr(), serial)
451      // fmt.Printf("Queue.Push(): %v\n", re)
452      // fmt.Printf("Queue.Push(): %v\n", re.Expr())
453      cnt_ins++
454      PS.Queue.Push(re)
455      //fmt.Printf("Testscore: %v\n", re.TestScore())
456      //PS.commup.Rpts <- &re     //sort missing! > 3 ergs
457    }
458  }
459
460  PS.Queue.Sort() //3 besten werden beim naechsten peel ausgegeben
461
462  for p := 0; p < PS.cnfg.peelCnt && PS.Queue.Len() > 0; p++ {
463    val := PS.Queue.Pop().(*probs.ExprReport)
464    inserts = append(inserts, val)
465  }
466
467  for _, e := range inserts {
468    PS.Queue.Push(e)
469  }
470
471  PS.Queue.Sort()
472
473  PS.commup.Res <- &inserts
474
475  PS.reportExpr()
476  PS.iter++
477
478  return len(inserts)
479}
480
481func (PS *PgeSearch) peel() []*probs.ExprReport {
482  es := make([]*probs.ExprReport, PS.cnfg.peelCnt)
483
484  rpt := make(probs.ExprReportArray, 0)
485
486  for p := 0; p < PS.cnfg.peelCnt && PS.Queue.Len() > 0; p++ {
487
488    e := PS.Queue.Pop().(*probs.ExprReport)
489
490    bPush := true
491    if len(e.Coeff()) == 1 && math.Abs(e.Coeff()[0]) < PS.cnfg.zeroEpsilon {
492      fmt.Println("No Best Push")
493      p--
494      continue
495    }
496
497    if bPush {
498      fmt.Printf("BEST PUSH/push(%d,%d): %v\n", p, PS.Best.Len(), e.Expr())
499      PS.Best.Push(e)
500      rpt = append(rpt, e)
501    }
502
503    es[p] = e
504
505    if e.TestScore() > PS.maxScore {
506      fmt.Printf("Testscore: %v\n", e.TestScore())
507      PS.maxScore = e.TestScore()
508    }
509    if e.TestError() < PS.minError {
510      PS.minError = e.TestError()
511      fmt.Printf("Best New Min Error:  %v\n", e)
512    }
513    if e.Size() > PS.maxSize {
514      PS.maxSize = e.Size()
515    }
516
517  }
518
519  //fmt.Printf("sand %d best pushes in peel\n", len(rpt))
520  //PS.commup.Res <- &rpt
521
522  _ = rpt
523
524  return es
525}
526
527func (PS *PgeSearch) expandPeeled(es []*probs.ExprReport) [][]expr.Expr {
528  eqns := make([][]expr.Expr, PS.cnfg.peelCnt)
529  for p := 0; p < PS.cnfg.peelCnt; p++ {
530    if es[p] == nil {
531      continue
532    }
533    // fmt.Printf("expand(%d): %v\n", p, es[p].Expr())
534    if es[p].Expr().ExprType() != expr.ADD {
535      add := expr.NewAdd()
536      add.Insert(es[p].Expr())
537      add.CalcExprStats()
538      es[p].SetExpr(add)
539    }
540    eqns[p] = PS.Expand(es[p].Expr())
541    // fmt.Printf("Results:\n")
542    // for i, e := range eqns[p] {
543    //  fmt.Printf("%d,%d:  %v\n", p, i, e)
544    // }
545    // fmt.Println()
546  }
547  fmt.Println("\n")
548  return eqns
549}
550
551//Hier ist das pgeRpt problem, liefert immer die besten pgerpt cnt ergebnisse zurueck
552// //*
553func (PS *PgeSearch) reportExpr() {
554
555  cnt := PS.cnfg.pgeRptCount
556  PS.Best.Sort()
557
558  // repot best equations
559  rpt := make(probs.ExprReportArray, cnt)
560  if PS.Best.Len() < cnt {
561    cnt = PS.Best.Len()
562  }
563  copy(rpt, PS.Best.GetQueue()[:cnt])
564
565  errSum, errCnt := 0.0, 0
566  for _, r := range rpt {
567    if r != nil && r.Expr() != nil {
568      errSum += r.TestError()
569      errCnt++
570    }
571  }
572
573  //PS.commup.Rpts <- &rpt
574
575}
576
577func (PS *PgeSearch) FirstPeel() {
578  PS.peel()
579}
580
581func (PS *PgeSearch) Clean() {
582
583}
584
585func (PS *PgeSearch) initLogs(logdir string) {
586
587}
588
589func (PS *PgeSearch) checkMessages() {
590
591  // check messages from superior
592  select {
593  case cmd, ok := <-PS.commup.Cmds:
594    if ok {
595      if cmd == -1 {
596        fmt.Println("PGE: stop sig recv'd")
597        PS.stop = true
598        return
599      }
600    }
601  default:
602    return
603  }
604}
605
606var c_input, c_ygiven []levmar.C_double
607
608func RegressExpr(E expr.Expr, P *probs.ExprProblem) (R *probs.ExprReport) {
609
610  guess := make([]float64, 0)
611  guess, eqn := E.ConvertToConstants(guess)
612
613  var coeff []float64
614  if len(guess) > 0 {
615
616    // fmt.Printf("x_dims:  %d  %d\n", x_dim, x_dim2)
617
618    // Callback version
619    coeff = levmar.LevmarExpr(eqn, P.SearchVar, P.SearchType, guess, P.Train, P.Test)
620
621    // Stack version
622    // x_dim := P.Train[0].NumDim()
623    // if c_input == nil {
624    //  ps := P.Train[0].NumPoints()
625    //  PS := len(P.Train) * ps
626    //  x_tot := PS * x_dim
627
628    //  c_input = make([]levmar.C_double, x_tot)
629    //  c_ygiven = make([]levmar.C_double, PS)
630
631    //  for i1, T := range P.Train {
632    //    for i2, p := range T.Points() {
633    //      i := i1*ps + i2
634    //      c_ygiven[i] = levmar.MakeCDouble(p.Depnd(P.SearchVar))
635    //      for i3, x_p := range p.Indeps() {
636    //        j := i1*ps*x_dim + i2*x_dim + i3
637    //        c_input[j] = levmar.MakeCDouble(x_p)
638    //      }
639    //    }
640    //  }
641    // }
642    // coeff = levmar.StackLevmarExpr(eqn, x_dim, guess, c_ygiven, c_input)
643
644    // serial := make([]int, 0)
645    // serial = eqn.StackSerial(serial)
646    // fmt.Printf("StackSerial: %v\n", serial)
647    // fmt.Printf("%v\n%v\n%v\n\n", eqn, coeff, steff)
648  }
649
650  R = new(probs.ExprReport)
651  R.SetExpr(eqn) /*.ConvertToConstantFs(coeff)*/
652  R.SetCoeff(coeff)
653  R.Expr().CalcExprStats()
654
655  // hitsL1, hitsL2, evalCnt, nanCnt, infCnt, l1_err, l2_err := scoreExpr(E, P, coeff)
656  _, _, _, trnNanCnt, _, trn_l1_err, _ := scoreExpr(E, P, P.Train, coeff)
657  _, _, tstEvalCnt, tstNanCnt, _, tst_l1_err, tst_l2_err := scoreExpr(E, P, P.Test, coeff)
658
659  R.SetTrainScore(trnNanCnt)
660  R.SetTrainError(trn_l1_err)
661
662  R.SetPredScore(tstNanCnt)
663  R.SetTestScore(tstEvalCnt)
664  R.SetTestError(tst_l1_err)
665  R.SetPredError(tst_l2_err)
666
667  return R
668}
669
670func scoreExpr(e expr.Expr, P *probs.ExprProblem, dataSets []*probs.PointSet, coeff []float64) (hitsL1, hitsL2, evalCnt, nanCnt, infCnt int, l1_err, l2_err float64) {
671  var l1_sum, l2_sum float64
672  for _, PS := range dataSets {
673    for _, p := range PS.Points() {
674      y := p.Depnd(P.SearchVar)
675      var out float64
676      if P.SearchType == probs.ExprBenchmark {
677        out = e.Eval(0, p.Indeps(), coeff, PS.SysVals())
678      } else if P.SearchType == probs.ExprDiffeq {
679        out = e.Eval(p.Indep(0), p.Indeps()[1:], coeff, PS.SysVals())
680      }
681
682      if math.IsNaN(out) {
683        nanCnt++
684        continue
685      } else if math.IsInf(out, 0) {
686        infCnt++
687        continue
688      } else {
689        evalCnt++
690      }
691
692      diff := out - y
693      l1_val := math.Abs(diff)
694      l2_val := diff * diff
695      l1_sum += l1_val
696      l2_sum += l2_val
697
698      if l1_val < P.HitRatio {
699        hitsL1++
700      }
701      if l2_val < P.HitRatio {
702        hitsL2++
703      }
704    }
705  }
706
707  if evalCnt == 0 {
708    l1_err = math.NaN()
709    l2_err = math.NaN()
710  } else {
711    fEvalCnt := float64(evalCnt + 1)
712    l1_err = l1_sum / fEvalCnt
713    l2_err = math.Sqrt(l2_sum / fEvalCnt)
714  }
715
716  return
717}
718
719func (PS *PgeSearch) CreateDS(maxGen int, pgeRptEpoch int, pgeRptCount int, pgeArchiveCap int, peelCnt int, evalrCount int, zeroEpsilon float64, initMethod string, growMethod string, sortType int) {
720  PS.id = 0
721
722  var PC PgeConfig
723
724  PC.pgeRptEpoch = pgeRptEpoch
725  PC.pgeRptCount = pgeRptCount
726  PC.pgeArchiveCap = pgeArchiveCap
727  PC.peelCnt = peelCnt
728  PC.evalrCount = evalrCount
729  PC.zeroEpsilon = zeroEpsilon
730  PC.initMethod = initMethod
731  PC.growMethod = growMethod
732
733  if sortType == 1 {
734    PC.sortType = probs.PESORT_PARETO_TRN_ERR
735  } else {
736    PC.sortType = probs.PESORT_PARETO_TST_ERR
737  }
738
739  PC.maxGen = maxGen
740  PS.cnfg = PC
741  PS.stop = false
742  PS.iter = 0
743}
744
745func (PS *PgeSearch) SetProb(ep *probs.ExprProblem) {
746  PS.prob = ep
747}
748func (PS *PgeSearch) GetIter() int {
749  return PS.iter
750}
Note: See TracBrowser for help on using the repository browser.