
# -*- coding: utf-8 -*-
import sys
import os
# sys.path.append('/home/puudeli/0Pilvi/webPy/pylib/bvp_solver')
webpy = os.environ.get('WEBPY')
sys.path.append(webpy + 'pylib/bvp_solver')
# import pylib.bvp_solver


from solver import solve
from ProblemDefinition import ProblemDefinition
from Solution import Solution
from template_generator import get_template


import numpy
import time
from pylab import ioff
# from math import sin, cos
from timeit import default_timer as timer
from opt_dfuns import odefun, bcfun, init_guess  # , gradodefun
from opt_params import par
from OptPiirto import opt_nosturi_piirto, tallenna_ohjaus


# guess palauttaa tilavektorin arvon ajanhetkellä t
def guess(t):
    (tmax, dim) = guess_arr.shape
    it = int((tmax - 1.0) / par.T_f * t)
    return guess_arr[it, :]


def ratkaisu(i, s_paiva):
    global guess_arr
    start = timer()
    solution = solve(
        problem,
        solution_guess=guess,
        tolerance=1.0e-5,
        method=6,
        trace=2,
        max_subintervals=500)

    # x: aika-akseli, jolla on plotdim pistettä alku- ja loppuhetken välillä
    x = numpy.linspace(problem.boundary_points[0], problem.boundary_points[1],
                       plotdim)
    # y: ratkaisun mukaiset arvot tilamuuttujille kussakin aikapisteessä
    y = solution(x)
    # ratkaisusta uusi alkuarvaus
    guess_arr = y.transpose()
    end = timer()

    # piirretään muuttujien kuvaajat ajan funktiona
    opt_nosturi_piirto(s_paiva, i, x, y, plotdim)
    # välitulostus, ratkaisun iterointiin kulunut aika ja
    # muutaman mielenkiintoisen parametrin arvo
    s1 = '\niteration: {0:2d} time: {1:.4g} '.format(i, end - start)
    s2 = ' C_W {0:.4g}; C_F1 {1:.4g}; C_F2 {2:.4g}'
    s2 = s2.format(par.C_W, par.C_F1, par.C_F2)
    print(s1 + s2)
    return x, y


# Interaktiivinen grafiikka pois käytöstä. 
ioff()


plotdim = int(48.0 * par.T_f)


guess_arr = numpy.array(
    [init_guess(float(i) * par.T_f / 100.0) for i in range(0, 100)])


# kahden pisteen reuna-arvotehtävän määrittely
problem = ProblemDefinition(
    num_ODE=4 * par.sysdim,
    num_parameters=0,
    num_left_boundary_conditions=2 * par.sysdim,
    boundary_points=(0, par.T_f),
    function=odefun,
    # function_derivative=gradodefun,
    boundary_conditions=bcfun)


paiva = time.localtime(time.time())
s_paiva = str(paiva.tm_year) + str(paiva.tm_mon) + str(paiva.tm_mday) + str(
    paiva.tm_hour) + str(paiva.tm_min)

maxiter = 50
last = 0
# Muutetaan optimointikriteerin parametreja ja lasketaan uusi ratkaisu
# käyttäen edellistä ratkaisua alkuarvauksena
# Alla lisätään iteraatio iteraatiolta energian käytön painotusta.
# Voisi olla mielenkiintoista muutella ilmanvastusta ja kitkoja.
try:
    for i in range(0, maxiter):
        par.C_F1 = 0.8 * par.C_F1
        par.C_F2 = 0.8 * par.C_F2
        x, y = ratkaisu(i, s_paiva)
        last = i
        tallenna_ohjaus(x, y, plotdim)

finally:
    print("se siitä")
    
# try: .. finally: rakenne
# Jos optimiratkaisua ei löydy, `bvp_solver` ilmoittaa suoritusvirheestä, 
# jonka seurauksena suoritus siirtyy kohtaan finally:

iteration:  0 time: 10.25  C_W 1; C_F1 8; C_F2 1.6

iteration:  1 time: 10.03  C_W 1; C_F1 6.4; C_F2 1.28

iteration:  2 time: 6.639  C_W 1; C_F1 5.12; C_F2 1.024

iteration:  3 time: 9.691  C_W 1; C_F1 4.096; C_F2 0.8192

iteration:  4 time: 9.912  C_W 1; C_F1 3.277; C_F2 0.6554

iteration:  5 time: 7.844  C_W 1; C_F1 2.621; C_F2 0.5243

iteration:  6 time: 7.98  C_W 1; C_F1 2.097; C_F2 0.4194

iteration:  7 time: 6.208  C_W 1; C_F1 1.678; C_F2 0.3355

iteration:  8 time: 6.131  C_W 1; C_F1 1.342; C_F2 0.2684

iteration:  9 time: 7.279  C_W 1; C_F1 1.074; C_F2 0.2147

Sisällysluettelo