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