1 | package problems |
---|
2 | |
---|
3 | import ( |
---|
4 | "bufio" |
---|
5 | "bytes" |
---|
6 | "fmt" |
---|
7 | "math" |
---|
8 | "math/rand" |
---|
9 | "os" |
---|
10 | "C" |
---|
11 | "unsafe" |
---|
12 | ) |
---|
13 | |
---|
14 | type Point struct { |
---|
15 | indep []float64 |
---|
16 | depnd []float64 |
---|
17 | } |
---|
18 | |
---|
19 | func (d *Point) NumIndep() int { return len(d.indep) } |
---|
20 | func (d *Point) SetNumIndep(sz int) { d.indep = make([]float64, sz) } |
---|
21 | func (d *Point) Indep(p int) float64 { return d.indep[p] } |
---|
22 | func (d *Point) SetIndep(p int, v float64) { d.indep[p] = v } |
---|
23 | func (d *Point) Indeps() []float64 { return d.indep } |
---|
24 | func (d *Point) SetIndeps(v []float64) { d.indep = v } |
---|
25 | func (d *Point) NumDepnd() int { return len(d.depnd) } |
---|
26 | func (d *Point) SetNumDepnd(sz int) { d.depnd = make([]float64, sz) } |
---|
27 | func (d *Point) Depnd(p int) float64 { return d.depnd[p] } |
---|
28 | func (d *Point) SetDepnd(p int, v float64) { d.depnd[p] = v } |
---|
29 | func (d *Point) Depnds() []float64 { return d.depnd } |
---|
30 | func (d *Point) SetDepnds(v []float64) { d.depnd = v } |
---|
31 | |
---|
32 | type PointSet struct { |
---|
33 | filename string |
---|
34 | id int |
---|
35 | |
---|
36 | numDim int |
---|
37 | indepNames []string |
---|
38 | depndNames []string |
---|
39 | sysNames []string |
---|
40 | |
---|
41 | dataPoints []Point |
---|
42 | sysVals []float64 |
---|
43 | } |
---|
44 | |
---|
45 | func (d *PointSet) FN() string { return d.filename } |
---|
46 | func (d *PointSet) SetFN(fn string) { d.filename = fn } |
---|
47 | func (d *PointSet) ID() int { return d.id } |
---|
48 | func (d *PointSet) SetID(id int) { d.id = id } |
---|
49 | func (d *PointSet) SetNumPoints(cnt int) { d.dataPoints = make([]Point, cnt) } |
---|
50 | |
---|
51 | func (d *PointSet) NumIndep() int { return len(d.indepNames) } |
---|
52 | func (d *PointSet) NumDepnd() int { return len(d.depndNames) } |
---|
53 | func (d *PointSet) NumDim() int { return d.numDim } // TODO check to see if TIME is a variable |
---|
54 | func (d *PointSet) SetNumDim(dim int) { d.numDim = dim } // TODO check to see if TIME is a variable |
---|
55 | func (d *PointSet) NumSys() int { return len(d.sysNames) } |
---|
56 | func (d *PointSet) NumPoints() int { return len(d.dataPoints) } |
---|
57 | |
---|
58 | func (d *PointSet) IndepName(xi int) string { return d.indepNames[xi] } |
---|
59 | func (d *PointSet) GetIndepNames() []string { return d.indepNames } |
---|
60 | func (d *PointSet) SetIndepNames(names []string) { d.indepNames = names } |
---|
61 | func (d *PointSet) DepndName(xi int) string { return d.depndNames[xi] } |
---|
62 | func (d *PointSet) GetDepndNames() []string { return d.depndNames } |
---|
63 | func (d *PointSet) SetDepndNames(names []string) { d.depndNames = names } |
---|
64 | |
---|
65 | func (d *PointSet) SysName(si int) string { return d.sysNames[si] } |
---|
66 | func (d *PointSet) GetSysNames() []string { return d.sysNames } |
---|
67 | func (d *PointSet) SetSysNames(names []string) { d.sysNames = names } |
---|
68 | func (d *PointSet) SetSysVals(sv []float64) { d.sysVals = sv } |
---|
69 | |
---|
70 | func (d *PointSet) Point(p int) *Point { return &(d.dataPoints[p]) } |
---|
71 | func (d *PointSet) Points() []Point { return d.dataPoints } |
---|
72 | func (d *PointSet) SetPoints(pts []Point) { d.dataPoints = pts } |
---|
73 | func (d *PointSet) SysVal(p int) float64 { return d.sysVals[p] } |
---|
74 | func (d *PointSet) SysVals() []float64 { return d.sysVals } |
---|
75 | |
---|
76 | // read function at end of file [ func (d *PointSet) Read(filename string) ] |
---|
77 | |
---|
78 | type PntSubset struct { |
---|
79 | ds *PointSet |
---|
80 | |
---|
81 | index []int |
---|
82 | input []Point |
---|
83 | output []Point |
---|
84 | sysVals []float64 |
---|
85 | } |
---|
86 | |
---|
87 | func (s *PntSubset) ID() int { return s.ds.id } |
---|
88 | func (s *PntSubset) DS() *PointSet { return s.ds } |
---|
89 | func (s *PntSubset) SetDS(ds *PointSet) { s.ds = ds } |
---|
90 | |
---|
91 | func (s *PntSubset) NumIndep() int { return s.ds.NumIndep() } |
---|
92 | func (s *PntSubset) NumDepnd() int { return s.ds.NumDepnd() } |
---|
93 | func (s *PntSubset) NumSys() int { return s.ds.NumSys() } |
---|
94 | func (s *PntSubset) NumPoints() int { return len(s.index) } |
---|
95 | |
---|
96 | func (s *PntSubset) SysVals() []float64 { return s.sysVals } |
---|
97 | func (s *PntSubset) SetSysVals(svls []float64) { s.sysVals = svls } |
---|
98 | |
---|
99 | func (s *PntSubset) Index(p int) int { return s.index[p] } |
---|
100 | func (s *PntSubset) Indexes() []int { return s.index } |
---|
101 | func (s *PntSubset) SetIndexes(idxs []int) { s.index = idxs } |
---|
102 | func (s *PntSubset) Input(p int) *Point { return &s.input[p] } |
---|
103 | func (s *PntSubset) Output(p int) *Point { return &s.output[p] } |
---|
104 | |
---|
105 | func (s *PntSubset) AddPoint(p int, input, output *Point) { |
---|
106 | s.index = append(s.index, p) |
---|
107 | s.input = append(s.input, *input) |
---|
108 | s.output = append(s.output, *output) |
---|
109 | } |
---|
110 | |
---|
111 | // using indexes, update the input/output data |
---|
112 | func (s *PntSubset) Refresh() { |
---|
113 | L := len(s.index) |
---|
114 | if s.input == nil { |
---|
115 | s.input = make([]Point, L) |
---|
116 | } |
---|
117 | if s.output == nil { |
---|
118 | s.output = make([]Point, L) |
---|
119 | } |
---|
120 | |
---|
121 | for i := 0; i < L; i++ { |
---|
122 | s.input[i] = *s.ds.Point(s.index[i]) |
---|
123 | if s.index[i]+1 >= s.ds.NumPoints() { |
---|
124 | continue |
---|
125 | } |
---|
126 | s.output[i] = *s.ds.Point(s.index[i] + 1) |
---|
127 | } |
---|
128 | } |
---|
129 | |
---|
130 | func (trainData *PointSet) InitTrainDataF(indepNames []string, depndNames []string, matrix unsafe.Pointer, nEntries int) { |
---|
131 | trainData.ReadPointSet("C:/Users/Hansi/Desktop/Data.trn") |
---|
132 | } |
---|
133 | |
---|
134 | func (trainData *PointSet) InitTestDataF(indepNames []string, depndNames []string, matrix unsafe.Pointer, nEntries int) { |
---|
135 | trainData.ReadPointSet("C:/Users/Hansi/Desktop/Data.tst") |
---|
136 | } |
---|
137 | |
---|
138 | func (trainData *PointSet) InitTrainData(indepNames []string, depndNames []string, matrix unsafe.Pointer, nEntries int) { |
---|
139 | trainData.indepNames = indepNames //First line in .trn file |
---|
140 | trainData.numDim = len(trainData.indepNames) |
---|
141 | trainData.depndNames = depndNames //Second line in .trn file |
---|
142 | |
---|
143 | fmt.Printf("Var Names = %v | %v\n", trainData.depndNames, trainData.indepNames) |
---|
144 | var nClum int = len(indepNames) + len(depndNames) |
---|
145 | |
---|
146 | for i := 0; i < nEntries; i++ { |
---|
147 | var pnt Point |
---|
148 | for j := 0; j < len(trainData.indepNames); j++ { |
---|
149 | elemNr := i * nClum + j |
---|
150 | item := (*float64) (unsafe.Pointer( uintptr(unsafe.Pointer(matrix)) + unsafe.Sizeof(float64(0)) * uintptr(elemNr) )) |
---|
151 | pnt.indep = append(pnt.indep, *item) |
---|
152 | |
---|
153 | if os.Getenv("PGEDEBUG") == "1" { |
---|
154 | fmt.Printf("Adding independend Test/Train Data (%v,%v): %v\n", i, j, *item) |
---|
155 | } |
---|
156 | } |
---|
157 | |
---|
158 | for j := len(trainData.indepNames); j < (len(trainData.indepNames) + len(trainData.depndNames)); j++ { |
---|
159 | elemNr := i * nClum + j |
---|
160 | item := (*float64) (unsafe.Pointer( uintptr(unsafe.Pointer(matrix)) + unsafe.Sizeof(float64(0)) * uintptr(elemNr) )) |
---|
161 | pnt.depnd = append(pnt.depnd, *item ) |
---|
162 | |
---|
163 | if os.Getenv("PGEDEBUG") == "1" { |
---|
164 | fmt.Printf("Adding dependend (right side) Test/Train Data (%v, %v): %v\n", i, j, *item) |
---|
165 | } |
---|
166 | } |
---|
167 | |
---|
168 | if len(pnt.indep) > 0 { |
---|
169 | trainData.dataPoints = append(trainData.dataPoints, pnt) |
---|
170 | } |
---|
171 | } |
---|
172 | fmt.Printf("Num Points: %v\n", len(trainData.dataPoints)) |
---|
173 | |
---|
174 | //if os.Getenv("PGEDEBUG") == "2" { |
---|
175 | // DebugPrint(trainData) |
---|
176 | //} |
---|
177 | } |
---|
178 | |
---|
179 | //func (trainData *PointSet) DebugPrint() { |
---|
180 | // for index, elem in range |
---|
181 | //} |
---|
182 | |
---|
183 | func (d *PointSet) ReadPointSet(filename string) { |
---|
184 | ftotal, err := os.OpenFile(filename, os.O_RDONLY, 0) |
---|
185 | if err != nil { |
---|
186 | fmt.Printf("err: %v\n", err) |
---|
187 | return |
---|
188 | } |
---|
189 | defer ftotal.Close() |
---|
190 | file := bufio.NewReader(ftotal) |
---|
191 | |
---|
192 | var word string |
---|
193 | |
---|
194 | // get independent variables (x_i...) |
---|
195 | for i := 0; ; i++ { |
---|
196 | _, err := fmt.Fscanf(file, "%s", &word) |
---|
197 | if err != nil { |
---|
198 | break |
---|
199 | } |
---|
200 | d.indepNames = append(d.indepNames, word) |
---|
201 | } |
---|
202 | d.numDim = len(d.indepNames) |
---|
203 | |
---|
204 | // get dependent variables (y_j...) |
---|
205 | for i := 0; ; i++ { |
---|
206 | _, err := fmt.Fscanf(file, "%s", &word) |
---|
207 | if err != nil { |
---|
208 | break |
---|
209 | } |
---|
210 | d.depndNames = append(d.depndNames, word) |
---|
211 | } |
---|
212 | |
---|
213 | fmt.Printf("Var Names = %v | %v\n", d.depndNames, d.indepNames) |
---|
214 | |
---|
215 | for i := 0; ; i++ { |
---|
216 | var pnt Point |
---|
217 | var dval, ival float64 |
---|
218 | if err != nil { |
---|
219 | break |
---|
220 | } |
---|
221 | |
---|
222 | for j := 0; j < len(d.indepNames); j++ { |
---|
223 | _, err = fmt.Fscanf(file, "%f", &ival) |
---|
224 | if err != nil { |
---|
225 | break |
---|
226 | } |
---|
227 | pnt.indep = append(pnt.indep, ival) |
---|
228 | } |
---|
229 | |
---|
230 | for j := 0; j < len(d.depndNames); j++ { |
---|
231 | _, err = fmt.Fscanf(file, "%f\n", &dval) |
---|
232 | if err != nil { |
---|
233 | break |
---|
234 | } |
---|
235 | pnt.depnd = append(pnt.depnd, dval) |
---|
236 | } |
---|
237 | |
---|
238 | if len(pnt.indep) > 0 { |
---|
239 | d.dataPoints = append(d.dataPoints, pnt) |
---|
240 | } |
---|
241 | if i%100 == 0 { |
---|
242 | fmt.Println("Point(%d): %v\n", i, pnt) |
---|
243 | } |
---|
244 | } |
---|
245 | fmt.Printf("Num Points: %v\n", len(d.dataPoints)) |
---|
246 | } |
---|
247 | |
---|
248 | func ReadBytesPointSet(ftotal []byte) (d *PointSet) { |
---|
249 | d = new(PointSet) |
---|
250 | // ftotal, err := os.OpenFile(filename, os.O_RDONLY, 0) |
---|
251 | // if err != nil { |
---|
252 | // fmt.Printf("err: %v\n", err) |
---|
253 | // return |
---|
254 | // } |
---|
255 | // defer ftotal.Close() |
---|
256 | var err error |
---|
257 | file := bytes.NewReader(ftotal) |
---|
258 | |
---|
259 | var word string |
---|
260 | |
---|
261 | // get independent variables (x_i...) |
---|
262 | for i := 0; ; i++ { |
---|
263 | _, err := fmt.Fscanf(file, "%s", &word) |
---|
264 | if err != nil { |
---|
265 | break |
---|
266 | } |
---|
267 | d.indepNames = append(d.indepNames, word) |
---|
268 | } |
---|
269 | d.numDim = len(d.indepNames) |
---|
270 | |
---|
271 | // get dependent variables (y_j...) |
---|
272 | for i := 0; ; i++ { |
---|
273 | _, err := fmt.Fscanf(file, "%s", &word) |
---|
274 | if err != nil { |
---|
275 | break |
---|
276 | } |
---|
277 | d.depndNames = append(d.depndNames, word) |
---|
278 | } |
---|
279 | |
---|
280 | fmt.Printf("Var Names = %v | %v\n", d.depndNames, d.indepNames) |
---|
281 | |
---|
282 | for i := 0; ; i++ { |
---|
283 | var pnt Point |
---|
284 | var dval, ival float64 |
---|
285 | if err != nil { |
---|
286 | break |
---|
287 | } |
---|
288 | |
---|
289 | for j := 0; j < len(d.indepNames); j++ { |
---|
290 | _, err = fmt.Fscanf(file, "%f", &ival) |
---|
291 | if err != nil { |
---|
292 | break |
---|
293 | } |
---|
294 | pnt.indep = append(pnt.indep, ival) |
---|
295 | } |
---|
296 | |
---|
297 | for j := 0; j < len(d.depndNames); j++ { |
---|
298 | _, err = fmt.Fscanf(file, "%f\n", &dval) |
---|
299 | if err != nil { |
---|
300 | break |
---|
301 | } |
---|
302 | pnt.depnd = append(pnt.depnd, dval) |
---|
303 | } |
---|
304 | |
---|
305 | if len(pnt.indep) > 0 { |
---|
306 | d.dataPoints = append(d.dataPoints, pnt) |
---|
307 | } |
---|
308 | if i%100 == 0 { |
---|
309 | fmt.Println("Point(%d): %v\n", i, pnt) |
---|
310 | } |
---|
311 | } |
---|
312 | fmt.Printf("Num Points: %v\n", len(d.dataPoints)) |
---|
313 | return d |
---|
314 | } |
---|
315 | |
---|
316 | func (d *PointSet) WritePointSet(filename string) { |
---|
317 | fmt.Printf("Writing file: %s\n", filename) |
---|
318 | ftotal, err := os.Create(filename) |
---|
319 | if err != nil { |
---|
320 | fmt.Printf("Error creating file: %s %v\n", filename, err) |
---|
321 | return |
---|
322 | } |
---|
323 | defer ftotal.Close() |
---|
324 | file := bufio.NewWriter(ftotal) |
---|
325 | defer file.Flush() |
---|
326 | |
---|
327 | // write independent variable names (x_i...) |
---|
328 | for i := 0; i < d.NumIndep(); i++ { |
---|
329 | if err != nil { |
---|
330 | fmt.Errorf("error writing pointset to file: %v\n", err) |
---|
331 | break |
---|
332 | } |
---|
333 | } |
---|
334 | fmt.Fprintln(file) |
---|
335 | |
---|
336 | // trite dependent variable names (y_j...) |
---|
337 | for i := 0; i < d.NumDepnd(); i++ { |
---|
338 | if err != nil { |
---|
339 | break |
---|
340 | } |
---|
341 | } |
---|
342 | fmt.Fprintln(file) |
---|
343 | |
---|
344 | // write points |
---|
345 | points := d.Points() |
---|
346 | for i := 0; i < d.NumPoints(); i++ { |
---|
347 | indep := points[i].Indeps() |
---|
348 | depnd := points[i].Depnds() |
---|
349 | for j := 0; j < len(indep); j++ { |
---|
350 | if err != nil { |
---|
351 | break |
---|
352 | } |
---|
353 | } |
---|
354 | for j := 0; j < len(depnd); j++ { |
---|
355 | if err != nil { |
---|
356 | break |
---|
357 | } |
---|
358 | } |
---|
359 | fmt.Fprintln(file) |
---|
360 | } |
---|
361 | } |
---|
362 | |
---|
363 | func SplitPointSetTrainTest(pnts *PointSet, pcnt_train float64, seed int) (train, test *PointSet) { |
---|
364 | |
---|
365 | train = new(PointSet) |
---|
366 | test = new(PointSet) |
---|
367 | |
---|
368 | train.filename, test.filename = pnts.filename, pnts.filename |
---|
369 | train.id, test.id = pnts.id, pnts.id |
---|
370 | train.indepNames, test.indepNames = pnts.indepNames, pnts.indepNames |
---|
371 | train.depndNames, test.depndNames = pnts.depndNames, pnts.depndNames |
---|
372 | train.sysNames, test.sysNames = pnts.sysNames, pnts.sysNames |
---|
373 | train.sysVals, test.sysVals = pnts.sysVals, pnts.sysVals |
---|
374 | |
---|
375 | L := len(pnts.dataPoints) |
---|
376 | Tst := int(float64(L) * (1.0 - pcnt_train)) |
---|
377 | |
---|
378 | tmp := make([]Point, L) |
---|
379 | copy(tmp, pnts.dataPoints) |
---|
380 | |
---|
381 | rand.Seed(int64(seed)) |
---|
382 | |
---|
383 | for i := 0; i < Tst; i++ { |
---|
384 | p := rand.Intn(L - i) |
---|
385 | tmp[i], tmp[p] = tmp[p], tmp[i] |
---|
386 | } |
---|
387 | |
---|
388 | test.dataPoints = tmp[:Tst] |
---|
389 | train.dataPoints = tmp[Tst:] |
---|
390 | |
---|
391 | return |
---|
392 | } |
---|
393 | |
---|
394 | func (d *PointSet) ReadLakeFile(filename string) { |
---|
395 | ftotal, err := os.OpenFile(filename, os.O_RDONLY, 0) |
---|
396 | if err != nil { |
---|
397 | fmt.Printf("err: %v\n", err) |
---|
398 | return |
---|
399 | } |
---|
400 | defer ftotal.Close() |
---|
401 | file := bufio.NewReader(ftotal) |
---|
402 | |
---|
403 | var word string |
---|
404 | |
---|
405 | // get independent variables (x_i...) |
---|
406 | for i := 0; ; i++ { |
---|
407 | _, err := fmt.Fscanf(file, "%s", &word) |
---|
408 | if err != nil { |
---|
409 | break |
---|
410 | } |
---|
411 | d.indepNames = append(d.indepNames, word) |
---|
412 | } |
---|
413 | d.numDim = len(d.indepNames) |
---|
414 | |
---|
415 | // // get dependent variables (y_j...) |
---|
416 | // for i := 0; ; i++ { |
---|
417 | // _, err := fmt.Fscanf(file, "%s", &word) |
---|
418 | // if err != nil { |
---|
419 | // break |
---|
420 | // } |
---|
421 | // d.depndNames = append(d.depndNames, word) |
---|
422 | // } |
---|
423 | |
---|
424 | // remove time names from indepNames |
---|
425 | d.indepNames = d.indepNames[2:] |
---|
426 | |
---|
427 | // fmt.Printf("Var Names = %v | %v\n", d.depndNames, d.indepNames) |
---|
428 | |
---|
429 | for i := 0; ; i++ { |
---|
430 | var pnt Point |
---|
431 | var ival float64 |
---|
432 | if err != nil { |
---|
433 | break |
---|
434 | } |
---|
435 | |
---|
436 | // read time values and disgaurd |
---|
437 | var dummy string |
---|
438 | for j := 0; j < 2; j++ { |
---|
439 | _, err = fmt.Fscanf(file, "%s", &dummy) |
---|
440 | // fmt.Println(i, dummy) |
---|
441 | if err != nil { |
---|
442 | // fmt.Println("err:", err) |
---|
443 | break |
---|
444 | } |
---|
445 | } |
---|
446 | |
---|
447 | // append dummy time value |
---|
448 | pnt.indep = append(pnt.indep, -0.1) |
---|
449 | |
---|
450 | // read real data |
---|
451 | for j := 0; j < len(d.indepNames); j++ { |
---|
452 | _, err = fmt.Fscanf(file, "%f\n", &ival) |
---|
453 | if err != nil { |
---|
454 | // fmt.Println("err:", err) |
---|
455 | break |
---|
456 | } |
---|
457 | // hack for PAR |
---|
458 | // if j == 1 { |
---|
459 | // ival /= 2000.0 * (math.Pi / 2.0) |
---|
460 | // } |
---|
461 | |
---|
462 | pnt.indep = append(pnt.indep, ival) |
---|
463 | } |
---|
464 | |
---|
465 | // fmt.Println(i, pnt) |
---|
466 | |
---|
467 | // for j := 0; j < len(d.depndNames); j++ { |
---|
468 | // _, err = fmt.Fscanf(file, "%f\n", &dval) |
---|
469 | // if err != nil { |
---|
470 | // break |
---|
471 | // } |
---|
472 | // pnt.depnd = append(pnt.depnd, dval) |
---|
473 | // } |
---|
474 | |
---|
475 | if len(pnt.indep) > 1 { |
---|
476 | for p := 0; p < len(pnt.indep); p++ { |
---|
477 | if math.IsNaN(pnt.indep[p]) { |
---|
478 | // fmt.Println("NaN @ ", i, p) |
---|
479 | goto skip |
---|
480 | } |
---|
481 | } |
---|
482 | d.dataPoints = append(d.dataPoints, pnt) |
---|
483 | } |
---|
484 | skip: |
---|
485 | } |
---|
486 | |
---|
487 | // calculate numerical derivatives |
---|
488 | calcDerivs(d.dataPoints) |
---|
489 | |
---|
490 | for i := 0; i < len(d.dataPoints); i++ { |
---|
491 | p := d.dataPoints[i] |
---|
492 | if len(p.Indeps()) == 0 || len(p.Depnds()) == 0 { |
---|
493 | fmt.Println("Bad Point @", i) |
---|
494 | } |
---|
495 | // if i%100 == 0 { |
---|
496 | // fmt.Printf("Point(%d): %v\n", i, d.dataPoints[i]) |
---|
497 | // } |
---|
498 | } |
---|
499 | |
---|
500 | // fmt.Printf("Num Points: %v\n", len(d.dataPoints)) |
---|
501 | } |
---|
502 | |
---|
503 | /* Calculate the first derivative of four points with: h = 0.25 |
---|
504 | * (from: http://www.trentfguidry.net/post/2010/09/04/Numerical-differentiation-formulas.aspx) |
---|
505 | * |
---|
506 | * xF0 = ( -3.0*xF4 + 16.0*xF3 - 36.0*xF2 + 48.0*xF1 - 25.0*xF0) / (12.0*h) |
---|
507 | * xF1 = ( xF4 - 6.0-xF3 + 18.0*xF2 - 10.0*xF1 - 3.0*xF0) / (12.0*h) |
---|
508 | * xF2 = ( -xF4 + 8.0*xF3 - 8.0*xF1 + xF0) / (12.0*h) |
---|
509 | * xF3 = ( 3.0*xF4 + 10.0*xF3 - 18.0*xF2 + 6.0*xF1 - xF0) / (12.0*h) |
---|
510 | * xF4 = ( 25.0*xF4 - 48.0*xF3 + 36.0*xF2 - 16.0*xF1 + 3.0*xF0) / (12.0*h) |
---|
511 | */ |
---|
512 | func calcDerivs(pts []Point) { |
---|
513 | h := 24.0 // / (24.0 * 60.0) |
---|
514 | |
---|
515 | NP := len(pts) |
---|
516 | ND := pts[0].NumIndep() |
---|
517 | |
---|
518 | // for summing and averaging on a point/variable~wise basis |
---|
519 | cnts := make([][]int, NP) |
---|
520 | vals := make([][]float64, NP) |
---|
521 | for p := 0; p < NP; p++ { |
---|
522 | cnts[p] = make([]int, ND) |
---|
523 | vals[p] = make([]float64, ND) |
---|
524 | } |
---|
525 | |
---|
526 | for p := 0; p < NP-5; p++ { |
---|
527 | for i := 0; i < ND; i++ { |
---|
528 | var F [5]float64 |
---|
529 | var dF [5]float64 |
---|
530 | |
---|
531 | for j := 0; j < 5; j++ { |
---|
532 | F[j] = pts[p+j].Indep(i) |
---|
533 | } |
---|
534 | dF[0] = (-3.0*F[4] + 16.0*F[3] - 36.0*F[2] + 48.0*F[1] - 25.0*F[0]) / (12.0 * h) |
---|
535 | dF[1] = (F[4] - 6.0 - F[3] + 18.0*F[2] - 10.0*F[1] - 3.0*F[0]) / (12.0 * h) |
---|
536 | dF[2] = (-F[4] + 8.0*F[3] - 8.0*F[1] + F[0]) / (12.0 * h) |
---|
537 | dF[3] = (3.0*F[4] + 10.0*F[3] - 18.0*F[2] + 6.0*F[1] - F[0]) / (12.0 * h) |
---|
538 | dF[4] = (25.0*F[4] - 48.0*F[3] + 36.0*F[2] - 16.0*F[1] + 3.0*F[0]) / (12.0 * h) |
---|
539 | |
---|
540 | for j := 0; j < 5; j++ { |
---|
541 | |
---|
542 | vals[p+j][i] += dF[j] |
---|
543 | cnts[p+j][i]++ |
---|
544 | } |
---|
545 | } |
---|
546 | } |
---|
547 | |
---|
548 | for p := 0; p < NP; p++ { |
---|
549 | depnds := make([]float64, ND) |
---|
550 | for i := 0; i < ND; i++ { |
---|
551 | depnds[i] = vals[p][i] / float64(cnts[p][i]) |
---|
552 | } |
---|
553 | pts[p].SetDepnds(depnds) |
---|
554 | } |
---|
555 | } |
---|