/pliant/math/curven.pli
 
 1  module "/pliant/language/unsafe.pli" 
 2  module "/pliant/math/functions.pli" 
 3   
 4  constant max_dim 32 
 5   
 6   
 7  type Curven 
 8    field Int p n 
 9    field (Array Array:Float) grid 
 10    field Array:Int factors 
 11    field Array:Float points 
 12   
 13   
 14  method c resize p n grid 
 15    arg_rw Curven c ; arg Int n ; arg (Array Array:Float) grid 
 16    check n=grid:size 
 17    check n<=max_dim 
 18    := n ; := p 
 19    grid := grid 
 20    c:factors size := n 
 21    var Int := p 
 22    for (var Int d) n-1 
 23      c:factors := shunt d=p  (c:factors d-1)*(c:grid d-1):size 
 24      *= c:grid:size 
 25    c:points size := t 
 26   
 27   
 28  method c define params point 
 29    arg_rw Curven c ; arg Array:Float params point 
 30    check params:size=c:n 
 31    check point:size=c:p 
 32    var Int := 0 
 33    for (var Int d) c:n-1 
 34      var Int index := 0 
 35      while index+1<c:grid:d:size and params:d>=(c:grid:index+1) 
 36        index += 1 
 37      if params:d<>c:grid:d:index 
 38        error error_id_unexpected "Point is not on the grid" 
 39      += index*c:factors:d 
 40    for (var Int p) c:p-1 
 41      c:points i+:= point p 
 42   
 43   
 44  method c in i d -> f 
 45    arg Curven c ; arg Int d ; arg Float f 
 46    var Int index := (i\c:factors:d)%c:grid:d:size 
 47    check index>and index<c:grid:d:size-1 
 48    var Float slope := ((c:points i+c:factors:d)-(c:points i-c:factors:d))/((c:grid:index+1)-(c:grid:index-1)) 
 49    := ((c:grid:index-1)-(c:grid:index))*slope 
 50   
 51  method c out i d -> f 
 52    arg Curven c ; arg Int d ; arg Float f 
 53    var Int index := (i\c:factors:d)%c:grid:d:size 
 54    check index>and index<c:grid:d:size-1 
 55    var Float slope := ((c:points i+c:factors:d)-(c:points i-c:factors:d))/((c:grid:index+1)-(c:grid:index-1)) 
 56    := ((c:grid:index+1)-(c:grid:index))*slope 
 57   
 58  method c pos i -> f 
 59    arg Curven c ; arg Int i ; arg Float f 
 60    := c:points i 
 61   
 62  function weigh x -> w 
 63    arg Float w 
 64    := 1-(cos x*pi/2) 
 65   
 66  method c apply params -> point 
 67    arg Curven c ; arg Array:Float params point 
 68    check params:size=c:n 
 69    var (Array Int max_dim) indices 
 70    var (Array Float max_dim) remains wremains w1remains 
 71    var Int base := 0 
 72    for (var Int d) c:n-1 
 73      var Int index := 0 
 74      while index+2<c:grid:d:size and params:d>=(c:grid:index+1) 
 75        index += 1 
 76      base += index*c:factors:d 
 77      indices := index 
 78      remains := (params:d-c:grid:d:index)/((c:grid:index+1)-(c:grid:index)) 
 79      wremains := weigh remains:d 
 80      w1remains := weigh 1-remains:d 
 81    point size := p 
 82    for (var Int p) point:size-1 
 83      point := 0 
 84    var Float := 0 
 85    for (var Int k) 0 2^c:n-1 
 86      var Float := 1 
 87      var Int := base 
 88      part try 
 89        for (var Int d) c:n-1 
 90          if (.and. 2^d)<>0 
 91            if (indices:d=c:grid:d:size-2) 
 92              leave try 
 93            eif indices:d=0 
 94              *= remains d 
 95            else 
 96              *= wremains d 
 97            += c:factors d 
 98          else 
 99            if indices:d=0 
 100              leave try 
 101            eif (indices:d=c:grid:d:size-2) 
 102              *= 1-remains:d 
 103            else 
 104              *= w1remains d 
 105        for (var Int p) c:p-1 
 106          point += w*(pos i+p) 
 107          for (var Int d) c:n-1 
 108            point += w*(shunt (.and. 2^d)<>0 (1-remains:d)*(in i+d) remains:d*(out i+d)) 
 109        += w 
 110    for (var Int k) 0 2^c:n-1 
 111      var Float := 1-t 
 112      var Int := base 
 113      for (var Int d) c:n-1 
 114        if (.and. 2^d)<>0 
 115          *= remains d 
 116          += c:factors d 
 117        else 
 118          *= 1-remains:d 
 119      for (var Int p) c:p-1 
 120        point += w*(pos i+p) 
 121   
 122   
 123  export Curven '. resize' '. define' '. apply' 
 124   
 125   
 126  if false 
 127    function test 
 128      var (Array Array:Float) grid 
 129      grid size := 1 
 130      grid:0 size := 5 
 131      for (var Int i) 0 grid:0:size-1 
 132        grid:0 i := i/(grid:0:size-1)*pi/2 
 133      var Curven c 
 134      c resize 2 1 grid 
 135      var Array:Float params 
 136      params size := 1 
 137      var Array:Float point 
 138      point size := 2 
 139      for (var Int i) 0 grid:0:size-1 
 140        params 0 := grid:0:i 
 141        point 0 := sin grid:0:i 
 142        point 1 := cos grid:0:i 
 143        c define params point 
 144      var Float m := 0 
 145      for (var Int i) 0 100 
 146        params 0 := i/100*pi/2 
 147        point := c apply params 
 148        var Float d0 := abs point:0-(sin params:0) 
 149        var Float d1 := abs point:1-(cos params:0) 
 150        console i " -> " d0 " " d1 " (" point:0 " " point:1 ")" eol 
 151        m := max m (max d0 d1) 
 152      console "maximum delta is " m eol 
 153    test