def test_get_ensemble(self):
        EPSG = 32633
        upper_left_x = 436100.0
        upper_left_y = 7417800.0
        nx = 74
        ny = 94
        dx = 1000.0
        dy = 1000.0
        # Period start
        year = 2015
        month = 7
        day = 26
        hour = 0
        n_hours = 30
        utc = api.Calendar()  # No offset gives Utc
        t0 = api.YMDhms(year, month, day, hour)
        period = api.UtcPeriod(utc.time(t0), utc.time(t0) + api.deltahours(n_hours))
        t_c = utc.time(t0) + api.deltahours(1)

        base_dir = path.join(shyftdata_dir, "netcdf", "arome")
        pattern = "fc*.nc"
        bbox = ([upper_left_x, upper_left_x + nx*dx,
                 upper_left_x + nx*dx, upper_left_x],
                [upper_left_y, upper_left_y,
                 upper_left_y - ny*dy, upper_left_y - ny*dy])
        try:
            repos = AromeDataRepository(EPSG, base_dir, filename=pattern, bounding_box=bbox)
            data_names = ("temperature", "wind_speed", "relative_humidity")
            ensemble = repos.get_forecast_ensemble(data_names, period, t_c, None)
            self.assertTrue(isinstance(ensemble, list))
            self.assertEqual(len(ensemble), 10)
        except AromeDataRepositoryError as adre:
            self.skipTest("(test inconclusive- missing arome-data {0})".format(adre))

    def test_get_ensemble(self):
        """
        Simple ensemble regression test of OpenDAP data repository.
        """
        epsg, bbox = self.epsg_bbox

        dem_file = path.join(shyftdata_dir, "netcdf", "etopo180.nc")

        # Period start
        (year, month, day), hour = self.start_date, 9
        n_hours = 30
        utc = api.Calendar()  # No offset gives Utc
        t0 = api.YMDhms(year, month, day, hour)
        period = api.UtcPeriod(utc.time(t0), utc.time(t0) + api.deltahours(n_hours))
        t_c = utc.time(t0) + api.deltahours(7)

        repos = GFSDataRepository(epsg, dem_file, bounding_box=bbox)
        data_names = ("temperature", "wind_speed", "precipitation",
                      "relative_humidity", "radiation")
        ensembles = repos.get_forecast_ensemble(data_names, period, t_c, None)
        for sources in ensembles:
            self.assertEqual(set(data_names), set(sources.keys()))
            self.assertEqual(len(sources["temperature"]), 6)
            data1 = sources["temperature"][0]
            data2 = sources["temperature"][1]
            self.assertNotEqual(data1.mid_point().x, data2.mid_point().x)
            self.assertNotEqual(data1.mid_point().y, data2.mid_point().y)
            self.assertNotEqual(data1.mid_point().z, data2.mid_point().z)
            h_dt = (data1.ts.time(1) - data1.ts.time(0))/3600
            self.assertEqual(data1.ts.size(), 30//h_dt)

    def test_get_ensemble_forecast_collection(self):
        EPSG = 32633
        upper_left_x = 436100.0
        upper_left_y = 7417800.0
        nx = 74
        ny = 94
        dx = 1000.0
        dy = 1000.0
        t0 = api.YMDhms(2015, 7, 26, 0)
        n_hours = 30
        utc = api.Calendar()  # No offset gives Utc
        period = api.UtcPeriod(utc.time(t0), utc.time(t0) + api.deltahours(n_hours))
        t_c = utc.time(t0) + api.deltahours(1)

        base_dir = path.join(shyftdata_dir, "netcdf", "arome")
        pattern = "fc2015072600.nc"
        bbox = ([upper_left_x, upper_left_x + nx * dx,
                 upper_left_x + nx * dx, upper_left_x],
                [upper_left_y, upper_left_y,
                 upper_left_y - ny * dy, upper_left_y - ny * dy])
        try:
            ar1 = AromeDataRepository(EPSG, base_dir, filename=pattern, bounding_box=bbox)
            ar2 = AromeDataRepository(EPSG, base_dir, filename=pattern, bounding_box=bbox)
            repos = GeoTsRepositoryCollection([ar1, ar2])
            data_names = ("temperature", "wind_speed", "relative_humidity")
            ensemble = repos.get_forecast_ensemble(data_names, period, t_c, None)
            self.assertTrue(isinstance(ensemble, list))
            self.assertEqual(len(ensemble), 10)
            with self.assertRaises(GeoTsRepositoryCollectionError) as context:
                repos = GeoTsRepositoryCollection([ar1, ar2], reduce_type="add")
                repos.get_forecast_ensemble(data_names, period, t_c, None)
            self.assertEqual("Only replace is supported yet", context.exception.args[0])
        except AromeDataRepositoryError as adre:
            self.skipTest("(test inconclusive- missing arome-data {0})".format(adre))

    def test_get_forecast_collection(self):
        n_hours = 30
        dt = api.deltahours(1)
        utc = api.Calendar()  # No offset gives Utc
        tc = api.YMDhms(2015, 8, 24, 6)
        t0 = utc.time(tc)
        period = api.UtcPeriod(t0, t0 + api.deltahours(n_hours))
        date_str = "{}{:02}{:02}_{:02}".format(tc.year, tc.month, tc.day, tc.hour)

        epsg, bbox = self.arome_epsg_bbox

        base_dir = path.join(shyftdata_dir, "repository", "arome_data_repository")
        f1 = "arome_metcoop_red_default2_5km_{}.nc".format(date_str)
        f2 = "arome_metcoop_red_test2_5km_{}.nc".format(date_str)

        ar1 = AromeDataRepository(epsg, base_dir, filename=f1, allow_subset=True)
        ar2 = AromeDataRepository(epsg, base_dir, filename=f2, elevation_file=f1, allow_subset=True)

        geo_ts_repository = GeoTsRepositoryCollection([ar1, ar2])
        source_names = ("temperature", "radiation")
        sources = geo_ts_repository.get_forecast(source_names, period, t0,
                                                 geo_location_criteria=bbox)
        self.assertTrue(all([x in source_names for x in sources]))

        geo_ts_repository = GeoTsRepositoryCollection([ar1, ar2], reduce_type="add")
        with self.assertRaises(GeoTsRepositoryCollectionError) as context:
            sources = geo_ts_repository.get_forecast(("temperature", "radiation"),
                                                     period, t0, geo_location_criteria=bbox)

    def test_get_forecast(self):
        # Period start
        year = 2015
        month = 8
        day = 24
        hour = 6
        n_hours = 65
        utc = api.Calendar()  # No offset gives Utc
        t0 = api.YMDhms(year, month, day, hour)
        period = api.UtcPeriod(utc.time(t0), utc.time(t0) + api.deltahours(n_hours))
        t_c1 = utc.time(t0) + api.deltahours(1)
        t_c2 = utc.time(t0) + api.deltahours(7)

        base_dir = path.join(shyftdata_dir, "repository", "arome_data_repository")
        pattern = "arome_metcoop*default2_5km_*.nc"
        EPSG, bbox = self.arome_epsg_bbox

        repos = AromeDataRepository(EPSG, base_dir, filename=pattern, bounding_box=bbox)
        data_names = ("temperature", "wind_speed", "precipitation", "relative_humidity")
        tc1_sources = repos.get_forecast(data_names, period, t_c1, None)
        tc2_sources = repos.get_forecast(data_names, period, t_c2, None)

        self.assertTrue(len(tc1_sources) == len(tc2_sources))
        self.assertTrue(set(tc1_sources) == set(data_names))
        self.assertTrue(tc1_sources["temperature"][0].ts.size() == n_hours + 1)

        tc1_precip = tc1_sources["precipitation"][0].ts
        tc2_precip = tc2_sources["precipitation"][0].ts

        self.assertEqual(tc1_precip.size(), n_hours)
        self.assertTrue(tc1_precip.time(0) != tc2_precip.time(0))

    def test_bias_predictor(self):
        """
        Verify that if we feed forecast[n] and observation into the bias-predictor
        it will create the estimated bias offsets
        """
        f = api.KalmanFilter()
        bp = api.KalmanBiasPredictor(f)
        self.assertIsNotNone(bp)
        self.assertEqual(bp.filter.parameter.n_daily_observations, 8)

        n_fc = 8
        utc = api.Calendar()
        t0 = utc.time(2016, 1, 1)
        dt = api.deltahours(1)
        n_fc_steps = 36  # e.g. like arome 36 hours
        fc_dt = api.deltahours(6)
        fc_fx = lambda time_axis: self._create_fc_values(time_axis, 2.0)  # just return a constant 2.0 deg C for now
        fc_set = self._create_geo_forecast_set(n_fc, t0, dt, n_fc_steps, fc_dt, fc_fx)
        n_obs = 24
        obs_ta = api.Timeaxis2(t0, dt, n_obs)
        obs_ts = api.Timeseries(obs_ta, fill_value=0.0)
        kalman_dt = api.deltahours(3)  # suitable average for prediction temperature
        kalman_ta = api.Timeaxis2(t0, kalman_dt, 8)
        bp.update_with_forecast(fc_set, obs_ts, kalman_ta)  # here we feed in forecast-set and observation into kalman
        fc_setv = self._create_forecast_set(n_fc, t0, dt, n_fc_steps, fc_dt, fc_fx)
        bp.update_with_forecast(fc_setv, obs_ts, kalman_ta)  # also verify we can feed in a pure TsVector
        bias_pattern = bp.state.x  # the bp.state.x is now the best estimates fo the bias between fc and observation
        self.assertEqual(len(bias_pattern), 8)
        for i in range(len(bias_pattern)):
            self.assertLess(abs(bias_pattern[i] - 2.0), 0.2)  # bias should iterate to approx 2.0 degC now.

def continuous_calibration():
    utc = Calendar()
    t_start = utc.time(YMDhms(2011, 9, 1))
    t_fc_start = utc.time(YMDhms(2015, 10, 1))
    dt = deltahours(1)
    n_obs = int(round((t_fc_start - t_start)/dt))
    obs_time_axis = Timeaxis(t_start, dt, n_obs + 1)
    q_obs_m3s_ts = observed_tistel_discharge(obs_time_axis.total_period())

    ptgsk = create_tistel_simulator(PTGSKOptModel, tistel.geo_ts_repository(tistel.grid_spec.epsg()))
    initial_state = burn_in_state(ptgsk, t_start, utc.time(YMDhms(2012, 9, 1)), q_obs_m3s_ts)

    num_opt_days = 30
    # Step forward num_opt_days days and store the state for each day:
    recal_start = t_start + deltahours(num_opt_days*24)
    t = t_start
    state = initial_state
    opt_states = {t: state}
    while t < recal_start:
        ptgsk.run(Timeaxis(t, dt, 24), state)
        t += deltahours(24)
        state = ptgsk.reg_model_state
        opt_states[t] = state

    recal_stop = utc.time(YMDhms(2011, 10, 30))
    recal_stop = utc.time(YMDhms(2012, 5, 30))
    curr_time = recal_start
    q_obs_avg = TsTransform().to_average(t_start, dt, n_obs + 1, q_obs_m3s_ts)
    target_spec = TargetSpecificationPts(q_obs_avg, IntVector([0]), 1.0, KLING_GUPTA)
    target_spec_vec = TargetSpecificationVector([target_spec])
    i = 0
    times = []
    values = []
    p, p_min, p_max = construct_calibration_parameters(ptgsk)
    while curr_time < recal_stop:
        print(i)
        i += 1
        opt_start = curr_time - deltahours(24*num_opt_days)
        opt_state = opt_states.pop(opt_start)
        p = ptgsk.region_model.get_region_parameter()
        p_opt = ptgsk.optimize(Timeaxis(opt_start, dt, 24*num_opt_days), opt_state, target_spec_vec,
                               p, p_min, p_max, tr_stop=1.0e-5)
        ptgsk.region_model.set_region_parameter(p_opt)
        corr_state = adjust_simulator_state(ptgsk, curr_time, q_obs_m3s_ts)
        ptgsk.run(Timeaxis(curr_time, dt, 24), corr_state)
        curr_time += deltahours(24)
        opt_states[curr_time] = ptgsk.reg_model_state
        discharge = ptgsk.region_model.statistics.discharge([0])
        times.extend(discharge.time(i) for i in range(discharge.size()))
        values.extend(list(np.array(discharge.v)))
    plt.plot(utc_to_greg(times), values)
    plot_results(None, q_obs=observed_tistel_discharge(UtcPeriod(recal_start, recal_stop)))
    set_calendar_formatter(Calendar())
    #plt.interactive(1)
    plt.title("Continuously recalibrated discharge vs observed")
    plt.xlabel("Time in UTC")
    plt.ylabel(r"Discharge in $\mathbf{m^3s^{-1}}$", verticalalignment="top", rotation="horizontal")
    plt.gca().yaxis.set_label_coords(0, 1.1)

 def test_extract_conversion_factors_from_string(self):
     u = utime('hours since 1970-01-01 00:00:00')
     t_origin = api.Calendar(u.tzoffset).time(
         api.YMDhms(u.origin.year, u.origin.month, u.origin.day, u.origin.hour, u.origin.minute, u.origin.second))
     delta_t_dic = {'days': api.deltahours(24), 'hours': api.deltahours(1), 'minutes': api.deltaminutes(1)}
     delta_t = delta_t_dic[u.units]
     self.assertIsNotNone(u)
     self.assertEqual(delta_t, api.deltahours(1))
     self.assertEqual(t_origin, 0)

 def _predict_bias(self, obs_set, fc_set):
     # Return a set of bias_ts per observation geo_point
     bias_set = api.TemperatureSourceVector()
     kf = api.KalmanFilter()
     kbp = api.KalmanBiasPredictor(kf)
     kta = api.Timeaxis2(self.t0, api.deltahours(3), 8)
     for obs in obs_set:
         kbp.update_with_forecast(fc_set, obs.ts, kta)
         pattern = api.KalmanState.get_x(kbp.state)
         # a_ts = api.Timeseries(pattern, api.deltahours(3), self.ta)  # can do using ct of Timeseries, or:
         b_ts = api.create_periodic_pattern_ts(pattern, api.deltahours(3), self.ta.time(0), self.ta)  # function
         bias_set.append(api.TemperatureSource(obs.mid_point(), b_ts))
     return bias_set

 def test_periodic_pattern_ts(self):
     c = api.Calendar()
     t0 = c.time(2016, 1, 1)
     dt = api.deltahours(1)
     n = 240
     ta = api.Timeaxis2(t0, dt, n)
     pattern_values = api.DoubleVector.from_numpy(np.arange(8))
     pattern_dt = api.deltahours(3)
     pattern_t0 = c.time(2015,6,1)
     pattern_ts = api.create_periodic_pattern_ts(pattern_values, pattern_dt, pattern_t0, ta)  # this is how to create a periodic pattern ts (used in gridpp/kalman bias handling)
     self.assertAlmostEqual(pattern_ts.value(0), 0.0)
     self.assertAlmostEqual(pattern_ts.value(1), 0.0)
     self.assertAlmostEqual(pattern_ts.value(2), 0.0)
     self.assertAlmostEqual(pattern_ts.value(3), 1.0)  # next step in pattern starts here
     self.assertAlmostEqual(pattern_ts.value(24), 0.0)  # next day repeats the pattern

    def test_subsets(self):
        EPSG, bbox = self.arome_epsg_bbox
        # Period start
        year = 2015
        month = 8
        day = 24
        hour = 6
        n_hours = 30
        date_str = "{}{:02}{:02}_{:02}".format(year, month, day, hour)
        utc = api.Calendar()  # No offset gives Utc
        t0 = api.YMDhms(year, month, day, hour)
        period = api.UtcPeriod(utc.time(t0), utc.time(t0) + api.deltahours(n_hours))

        base_dir = path.join(shyftdata_dir, "repository", "arome_data_repository")
        filename = "arome_metcoop_red_default2_5km_{}.nc".format(date_str)

        data_names = ("temperature", "wind_speed", "precipitation", "relative_humidity", "radiation")
        allow_subset = False
        reader = AromeDataRepository(EPSG, base_dir, filename=filename,
                                     bounding_box=bbox, allow_subset=allow_subset)
        with self.assertRaises(AromeDataRepositoryError) as context:
            reader.get_timeseries(data_names, period, None)
        self.assertEqual("Could not find all data fields", context.exception.args[0])
        allow_subset = True
        reader = AromeDataRepository(EPSG, base_dir, filename=filename,
                                     bounding_box=bbox, allow_subset=allow_subset)
        try:
            sources = reader.get_timeseries(data_names, period, None)
        except AromeDataRepositoryError as e:
            self.fail("AromeDataRepository.get_timeseries(data_names, period, None) "
                      "raised AromeDataRepositoryError unexpectedly.")
        self.assertEqual(len(sources), len(data_names) - 1)

    def fetch_sources(self, input_source_types, data, params, period):
        """Method for fetching the sources in NetCDF files.

        Parameters
        ----------
        input_source_types : dict
            A map between the data to be extracted and the data containers in shyft.api.
        data : dict
            An geo-located time series shyft.api container.
        params : dict
            Additional parameters for locating the datasets.
        period : tuple
            A (start_time, stop_time) tuple that species the simulation period.

        """
        self.__dict__.update(params)
        # Fill the data with actual values
        for input_source, source_api in input_source_types.iteritems():
            ts = self._fetch_station_tseries(input_source, params['types'], period)
            assert type(ts) is list
            tsf = api.TsFactory()
            for station in ts:
                times = station['time']
                assert type(times) is list
                dt = times[1] - times[0] if len(times) > 1 else api.deltahours(1)
                total_period = api.UtcPeriod(times[0], times[-1] + dt)
                time_points = api.UtcTimeVector(times)
                time_points.push_back(total_period.end)
                values = station['values']
                value_points = api.DoubleVector.FromNdArray(values)
                api_ts = tsf.create_time_point_ts(total_period, time_points, value_points)
                data_source = source_api(api.GeoPoint(*station['location']), api_ts)
                data[input_source].append(data_source)
        return data

    def test_kling_gupta_and_nash_sutcliffe(self):
        """
        Test/verify exposure of the kling_gupta and nash_sutcliffe correlation functions

        """

        def np_nash_sutcliffe(o, p):
            return 1 - (np.sum((o - p) ** 2)) / (np.sum((o - np.mean(o)) ** 2))

        c = api.Calendar()
        t0 = c.time(2016, 1, 1)
        dt = api.deltahours(1)
        n = 240
        ta = api.Timeaxis2(t0, dt, n)
        from math import sin, pi
        rad_max = 10 * 2 * pi
        obs_values = api.DoubleVector.from_numpy(np.array([sin(i * rad_max / n) for i in range(n)]))
        mod_values = api.DoubleVector.from_numpy(np.array([0.1 + sin(pi / 10.0 + i * rad_max / n) for i in range(n)]))
        obs_ts = api.Timeseries(ta=ta, values=obs_values, point_fx=api.point_interpretation_policy.POINT_AVERAGE_VALUE)
        mod_ts = api.Timeseries(ta=ta, values=mod_values, point_fx=api.point_interpretation_policy.POINT_AVERAGE_VALUE)

        self.assertAlmostEqual(api.kling_gupta(obs_ts, obs_ts, ta, 1.0, 1.0, 1.0), 1.0, None, "1.0 for perfect match")
        self.assertAlmostEqual(api.nash_sutcliffe(obs_ts, obs_ts, ta), 1.0, None, "1.0 for perfect match")
        # verify some non trivial cases, and compare to numpy version of ns
        mod_inv = obs_ts * -1.0
        kge_inv = obs_ts.kling_gupta(mod_inv)  # also show how to use time-series.method itself to ease use
        ns_inv = obs_ts.nash_sutcliffe(mod_inv)  # similar for nash_sutcliffe, you can reach it directly from a ts
        ns_inv2 = np_nash_sutcliffe(obs_ts.values.to_numpy(), mod_inv.values.to_numpy())
        self.assertLessEqual(kge_inv, 1.0, "should be less than 1")
        self.assertLessEqual(ns_inv, 1.0, "should be less than 1")
        self.assertAlmostEqual(ns_inv, ns_inv2, 4, "should equal numpy calculated value")
        kge_obs_mod = api.kling_gupta(obs_ts, mod_ts, ta, 1.0, 1.0, 1.0)
        self.assertLessEqual(kge_obs_mod, 1.0)
        self.assertAlmostEqual(obs_ts.nash_sutcliffe( mod_ts), np_nash_sutcliffe(obs_ts.values.to_numpy(), mod_ts.values.to_numpy()))

    def test_tiny_bbox(self):
        EPSG, _ = self.arome_epsg_bbox

        x0 = 436250.0   # lower left
        y0 = 6823250.0  # lower right
        nx = 1
        ny = 1
        dx = 5.0
        dy = 5.0
        bbox = ([x0, x0 + nx*dx, x0 + nx*dx, x0], [y0, y0, y0 + ny*dy, y0 + ny*dy])
        print(bbox)
        
        # Period start
        year = 2015
        month = 8
        day = 24
        hour = 6
        n_hours = 30
        date_str = "{}{:02}{:02}_{:02}".format(year, month, day, hour)
        utc = api.Calendar()  # No offset gives Utc
        t0 = api.YMDhms(year, month, day, hour)
        period = api.UtcPeriod(utc.time(t0), utc.time(t0) + api.deltahours(n_hours))

        base_dir = path.join(shyftdata_dir, "repository", "arome_data_repository")
        filename = "arome_metcoop_red_default2_5km_{}.nc".format(date_str)
        reader = AromeDataRepository(EPSG, base_dir, filename=filename, 
                                     bounding_box=bbox, x_padding=0, y_padding=0)
        data_names = ("temperature", "wind_speed", "precipitation", "relative_humidity")
        try:
            tss = reader.get_timeseries(data_names, period, None)
        except AromeDataRepository as err:
            self.fail("reader.get_timeseries raised AromeDataRepositoryError('{}') "
                      "unexpectedly.".format(err.args[0]))

 def test_create_region_environment(self):
     cal = api.Calendar()
     time_axis = api.Timeaxis(cal.time(api.YMDhms(2015, 1, 1, 0, 0, 0)), api.deltahours(1), 240)
     re = self.create_dummy_region_environment(time_axis, api.GeoPoint(1000, 1000, 100))
     self.assertIsNotNone(re)
     self.assertEqual(len(re.radiation), 1)
     self.assertAlmostEqual(re.radiation[0].ts.value(0), 300.0)

    def test_percentiles(self):
        c = api.Calendar()
        t0 = c.time(2016, 1, 1)
        dt = api.deltahours(1)
        n = 240
        ta = api.Timeaxis(t0, dt, n)
        timeseries = api.TsVector()

        for i in range(10):
            timeseries.append(
                api.Timeseries(ta=ta, fill_value=i, point_fx=api.point_interpretation_policy.POINT_AVERAGE_VALUE))

        wanted_percentiles = api.IntVector([0, 10, 50, -1, 70, 100])
        ta_day = api.Timeaxis(t0, dt * 24, n // 24)
        ta_day2 = api.Timeaxis2(t0, dt * 24, n // 24)
        percentiles = api.percentiles(timeseries, ta_day, wanted_percentiles)
        percentiles2 = timeseries.percentiles(ta_day2, wanted_percentiles)  # just to verify it works with alt. syntax

        self.assertEqual(len(percentiles2), len(percentiles))

        for i in range(len(ta_day)):
            self.assertAlmostEqual(0.0, percentiles[0].value(i), 3, "  0-percentile")
            self.assertAlmostEqual(0.9, percentiles[1].value(i), 3, " 10-percentile")
            self.assertAlmostEqual(4.5, percentiles[2].value(i), 3, " 50-percentile")
            self.assertAlmostEqual(4.5, percentiles[3].value(i), 3, "   -average")
            self.assertAlmostEqual(6.3, percentiles[4].value(i), 3, " 70-percentile")
            self.assertAlmostEqual(9.0, percentiles[5].value(i), 3, "100-percentile")

 def prec_acc_conv(p):
     indx = np.nonzero([self.cal.calendar_units(ti).hour in self.analysis_hours for ti in time])[0]
     f = 1000.*api.deltahours(1)/(time[1] - time[0]) # conversion from m/delta_t to mm/1hour
     dp = (p[1:] - p[:-1])*f
     dp[indx] = p[indx+1]*f
     #return np.clip(p[1:] - p[:-1], 0.0, 1000.0)
     return dp

 def test_calendar_add_3h_during_dst(self):
     osl = api.Calendar("Europe/Oslo")
     t0 = osl.time(2016, 3, 27)  # dst change during spring
     t1 = osl.add(t0, api.Calendar.DAY,  1)
     dt3h=api.deltahours(3)
     d3h= osl.diff_units(t0,t1,dt3h)
     self.assertEqual(8, d3h)

 def setUp(self):
     self.c=api.Calendar()
     self.d=api.deltahours(1)
     self.n=24
     #self.t= self.c.trim(api.utctime_now(),self.d)
     self.t= self.c.trim(self.c.time(api.YMDhms(1969,12,31,0,0,0)),self.d)
     self.ta=api.Timeaxis(self.t,self.d,self.n)

def ensemble_demo():
    utc = Calendar()
    t_start = utc.time(YMDhms(2011, 9, 1))
    t_fc_ens_start = utc.time(YMDhms(2015, 7, 26))
    disp_start = utc.time(YMDhms(2015, 7, 20))
    dt = deltahours(1)
    n_obs = int(round((t_fc_ens_start - t_start)/dt))
    n_fc_ens = 30
    n_disp = int(round(t_fc_ens_start - disp_start)/dt) + n_fc_ens + 24*7

    obs_time_axis = Timeaxis(t_start, dt, n_obs + 1)
    fc_ens_time_axis = Timeaxis(t_fc_ens_start, dt, n_fc_ens)
    display_time_axis = Timeaxis(disp_start, dt, n_disp)

    q_obs_m3s_ts = observed_tistel_discharge(obs_time_axis.total_period())
    ptgsk = create_tistel_simulator(PTGSKOptModel, tistel.geo_ts_repository(tistel.grid_spec.epsg()))
    initial_state = burn_in_state(ptgsk, t_start, utc.time(YMDhms(2012, 9, 1)), q_obs_m3s_ts)

    ptgsk.run(obs_time_axis, initial_state)
    current_state = adjust_simulator_state(ptgsk, t_fc_ens_start, q_obs_m3s_ts)
    q_obs_m3s_ts = observed_tistel_discharge(display_time_axis.total_period())
    ens_repos = tistel.arome_ensemble_repository(tistel.grid_spec)
    ptgsk_fc_ens = create_tistel_simulator(PTGSKModel, ens_repos)
    sims = ptgsk_fc_ens.create_ensembles(fc_ens_time_axis, t_fc_ens_start, current_state)
    for sim in sims:
        sim.simulate()
    plt.hold(1)
    percentiles = [10, 25, 50, 75, 90]
    plot_percentiles(sims, percentiles, obs=q_obs_m3s_ts)
    #plt.interactive(1)
    plt.show()