/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 p n ; arg (Array Array:Float) grid
16	check n=grid:size
17	check n<=max_dim
18	c n := n ; c p := p
19	c grid := grid
20	c:factors size := n
21	var Int t := p
22	for (var Int d) 0 n-1
23	c:factors d := shunt d=0 p  (c:factors d-1)*(c:grid d-1):size
24	t *= c:grid:d 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 i := 0
33	for (var Int d) 0 c:n-1
34	var Int index := 0
35	while index+1<c:grid:d:size and params:d>=(c:grid:d 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	i += index*c:factors:d
40	for (var Int p) 0 c:p-1
41	c:points i+p := point p
42
43
44	method c in i d -> f
45	arg Curven c ; arg Int i d ; arg Float f
46	var Int index := (i\c:factors:d)%c:grid:d:size
47	check index>0 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:d index+1)-(c:grid:d index-1))
49	f := ((c:grid:d index-1)-(c:grid:d index))*slope
50
51	method c out i d -> f
52	arg Curven c ; arg Int i d ; arg Float f
53	var Int index := (i\c:factors:d)%c:grid:d:size
54	check index>0 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:d index+1)-(c:grid:d index-1))
56	f := ((c:grid:d index+1)-(c:grid:d index))*slope
57
58	method c pos i -> f
59	arg Curven c ; arg Int i ; arg Float f
60	f := c:points i
61
62	function weigh x -> w
63	arg Float x w
64	w := 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) 0 c:n-1
73	var Int index := 0
74	while index+2<c:grid:d:size and params:d>=(c:grid:d index+1)
75	index += 1
76	base += index*c:factors:d
77	indices d := index
78	remains d := (params:d-c:grid:d:index)/((c:grid:d index+1)-(c:grid:d index))
79	wremains d := weigh remains:d
80	w1remains d := weigh 1-remains:d
81	point size := c p
82	for (var Int p) 0 point:size-1
83	point p := 0
84	var Float t := 0
85	for (var Int k) 0 2^c:n-1
86	var Float w := 1
87	var Int i := base
88	part try
89	for (var Int d) 0 c:n-1
90	if (k .and. 2^d)<>0
91	if (indices:d=c:grid:d:size-2)
92	leave try
93	eif indices:d=0
94	w *= remains d
95	else
96	w *= wremains d
97	i += c:factors d
98	else
99	if indices:d=0
100	leave try
101	eif (indices:d=c:grid:d:size-2)
102	w *= 1-remains:d
103	else
104	w *= w1remains d
105	for (var Int p) 0 c:p-1
106	point p += w*(c pos i+p)
107	for (var Int d) 0 c:n-1
108	point p += w*(shunt (k .and. 2^d)<>0 (1-remains:d)*(c in i+p d) remains:d*(c out i+p d))
109	t += w
110	for (var Int k) 0 2^c:n-1
111	var Float w := 1-t
112	var Int i := base
113	for (var Int d) 0 c:n-1
114	if (k .and. 2^d)<>0
115	w *= remains d
116	i += c:factors d
117	else
118	w *= 1-remains:d
119	for (var Int p) 0 c:p-1
120	point p += w*(c 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