<commit_before>from __future__ import division

import warnings

from six.moves import zip

warnings.filterwarnings("once", category=DeprecationWarning)

import logging
logging.basicConfig(format='%(asctime)s: %(name)s: %(message)s',
                    datefmt='%Y-%m-%d %H:%M:%S', level=logging.DEBUG)

import unittest
import os
import copy

import shapely.geometry as shpg
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Local imports
from oggm.tests.test_graphics import init_hef
from oggm.core.models import massbalance, flowline
from oggm.core.models.massbalance import ConstantBalanceModel
from oggm.tests import is_slow, requires_working_conda
from oggm.tests import HAS_NEW_GDAL

from oggm import utils, cfg
from oggm.cfg import N, SEC_IN_DAY, SEC_IN_MONTH, SEC_IN_YEAR

# Globals
current_dir = os.path.dirname(os.path.abspath(__file__))

# test directory
testdir = os.path.join(current_dir, 'tmp')

do_plot = False

DOM_BORDER = 80

def dummy_constant_bed(hmax=3000., hmin=1000., nx=200):

    map_dx = 100.
    dx = 1.

    surface_h = np.linspace(hmax, hmin, nx)
    bed_h = surface_h
    widths = surface_h * 0. + 3.

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.VerticalWallFlowline(line, dx, map_dx, surface_h,
                                          bed_h, widths)]

def dummy_constant_bed_cliff(hmax=3000., hmin=1000., nx=200):
    """
    I introduce a cliff in the bed to test the mass conservation of the models
    Such a cliff could be real or a DEM error/artifact
    """
    
    map_dx = 100.
    dx = 1.

    surface_h = np.linspace(hmax, hmin, nx)
    
    cliff_height = 250.0
    surface_h[50:] = surface_h[50:] - cliff_height
    
    bed_h = surface_h
    widths = surface_h * 0. + 1.

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.VerticalWallFlowline(line, dx, map_dx, surface_h,
                                          bed_h, widths)]

def dummy_bumpy_bed():

    map_dx = 100.
    dx = 1.
    nx = 200

    coords = np.arange(0, nx-0.5, 1)
    surface_h = np.linspace(3000, 1000, nx)
    surface_h += 170. * np.exp(-((coords-30)/5)**2)

    bed_h = surface_h
    widths = surface_h * 0. + 3.
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.VerticalWallFlowline(line, dx, map_dx, surface_h,
                                          bed_h, widths)]


def dummy_noisy_bed():

    map_dx = 100.
    dx = 1.
    nx = 200
    np.random.seed(42)
    coords = np.arange(0, nx-0.5, 1)
    surface_h = np.linspace(3000, 1000, nx)
    surface_h += 100 * np.random.rand(nx) - 50.

    bed_h = surface_h
    widths = surface_h * 0. + 3.
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.VerticalWallFlowline(line, dx, map_dx, surface_h,
                                          bed_h, widths)]

def dummy_parabolic_bed():

    map_dx = 100.
    dx = 1.
    nx = 200

    surface_h = np.linspace(3000, 1000, nx)
    bed_h = surface_h
    shape = surface_h * 0. + 5.e-03

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.ParabolicFlowline(line, dx, map_dx, surface_h,
                                       bed_h, shape)]

def dummy_mixed_bed():

    map_dx = 100.
    dx = 1.
    nx = 200

    surface_h = np.linspace(3000, 1000, nx)
    bed_h = surface_h
    shape = surface_h * 0. + 3.e-03
    shape[10:20] = 0.00001

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.MixedFlowline(line, dx, map_dx, surface_h,
                                   bed_h, shape, lambdas=3.5)]


def dummy_trapezoidal_bed(hmax=3000., hmin=1000., nx=200):

    map_dx = 100.
    dx = 1.

    surface_h = np.linspace(hmax, hmin, nx)
    bed_h = surface_h
    widths = surface_h * 0. + 1.6

    lambdas = surface_h * 0. + 2

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)

    return [flowline.TrapezoidalFlowline(line, dx, map_dx, surface_h,
                                         bed_h, widths, lambdas)]


def dummy_width_bed():
    """This bed has a width of 6 during the first 20 points and then 3"""

    map_dx = 100.
    dx = 1.
    nx = 200

    surface_h = np.linspace(3000, 1000, nx)
    bed_h = surface_h
    widths = surface_h * 0. + 3.
    widths[0:20] = 6.

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.VerticalWallFlowline(line, dx, map_dx, surface_h,
                                          bed_h, widths)]

def dummy_width_bed_tributary():

    map_dx = 100.
    dx = 1.
    nx = 200

    surface_h = np.linspace(3000, 1000, nx)
    bed_h = surface_h
    widths = surface_h * 0. + 3.
    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)

    fl_0 = flowline.VerticalWallFlowline(line, dx, map_dx, surface_h, bed_h,
                                         widths)
    coords = np.arange(0, 19.1, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.+1]).T)
    fl_1 = flowline.VerticalWallFlowline(line, dx, map_dx, surface_h[0:20],
                                         bed_h[0:20], widths[0:20])
    fl_1.set_flows_to(fl_0)
    return [fl_1, fl_0]

class TestInitFlowline(unittest.TestCase):

    def setUp(self):
        pass

    def tearDown(self):
        pass

    def test_init_present_time_glacier(self):

        gdir = init_hef(border=DOM_BORDER)
        flowline.init_present_time_glacier(gdir)

        fls = gdir.read_pickle('model_flowlines')

        lens = [len(gdir.read_pickle('centerlines', div_id=i)) for i in [1,
                                                                         2, 3]]
        ofl = gdir.read_pickle('inversion_flowlines',
                               div_id=np.argmax(lens)+1)[-1]

        self.assertTrue(gdir.rgi_date.year == 2003)
        self.assertTrue(len(fls) == 5)

        vol = 0.
        area = 0.
        for fl in fls:
            refo = 1 if fl is fls[-1] else 0
            self.assertTrue(fl.order == refo)
            ref = np.arange(len(fl.surface_h)) * fl.dx
            np.testing.assert_allclose(ref, fl.dis_on_line,
                                       rtol=0.001,
                                       atol=0.01)
            self.assertTrue(len(fl.surface_h) ==
                            len(fl.bed_h) ==
                            len(fl.bed_shape) ==
                            len(fl.dis_on_line) ==
                            len(fl.widths))

            self.assertTrue(np.all(fl.bed_shape >= 0))
            self.assertTrue(np.all(fl.widths >= 0))
            vol += fl.volume_km3
            area += fl.area_km2

            if refo == 1:
                rtol = 0.07 if HAS_NEW_GDAL else 0.05
                np.testing.assert_allclose(ofl.widths * gdir.grid.dx,
                                           fl.widths_m[0:len(ofl.widths)],
                                           rtol=rtol)

        np.testing.assert_allclose(0.573, vol, rtol=0.001)
        np.testing.assert_allclose(7350.0, fls[-1].length_m)

        rtol = 0.01 if HAS_NEW_GDAL else 0.001
        np.testing.assert_allclose(gdir.rgi_area_km2, area, rtol=rtol)

        if do_plot:  # pragma: no cover
            plt.plot(fls[-1].bed_h)
            plt.plot(fls[-1].surface_h)
            plt.show()

    def test_present_time_glacier_massbalance(self):

        gdir = init_hef(border=DOM_BORDER)
        flowline.init_present_time_glacier(gdir)

        mb_mod = massbalance.HistalpMassBalanceModel(gdir)

        fls = gdir.read_pickle('model_flowlines')
        glacier = flowline.FlowlineModel(fls)

        hef_file = utils.get_demo_file('mbdata_RGI40-11.00897.csv')
        mbdf = pd.read_csv(hef_file).set_index('YEAR')

        hgts = np.array([])
        widths = np.array([])
        for fl in glacier.fls:
            hgts = np.concatenate((hgts, fl.surface_h))
            widths = np.concatenate((widths, fl.widths_m))
        tot_mb = []
        refmb = []
        grads = hgts * 0
        for yr, mb in mbdf.iterrows():
            refmb.append(mb['ANNUAL_BALANCE'])
            mbh = mb_mod.get_mb(hgts, yr) * SEC_IN_YEAR * cfg.RHO
            grads += mbh
            tot_mb.append(np.average(mbh, weights=widths))
        grads /= len(tot_mb)

        # Bias
        self.assertTrue(np.abs(utils.md(tot_mb, refmb)) < 50)

        # Gradient
        dfg = pd.read_csv(utils.get_demo_file('mbgrads_RGI40-11.00897.csv'),
                          index_col='ALTITUDE').mean(axis=1)

        # Take the altitudes below 3100 and fit a line
        dfg = dfg[dfg.index < 3100]
        pok = np.where(hgts < 3100)
        from scipy.stats import linregress
        slope_obs, _, _, _, _ = linregress(dfg.index, dfg.values)
        slope_our, _, _, _, _ = linregress(hgts[pok], grads[pok])
        np.testing.assert_allclose(slope_obs, slope_our, rtol=0.1)


class TestMassBalance(unittest.TestCase):
    def setUp(self):
        pass

    def tearDown(self):
        pass

    def test_tstar_mb(self):

        gdir = init_hef(border=DOM_BORDER)
        flowline.init_present_time_glacier(gdir)

        mb_mod = massbalance.TstarMassBalanceModel(gdir)

        ofl = gdir.read_pickle('inversion_flowlines', div_id=0)
        h = np.array([])
        w = np.array([])
        for fl in ofl:
            h = np.append(h, fl.surface_h)
            w = np.append(w, fl.widths)

        ombh = mb_mod.get_mb(h, None)
        otmb = np.sum(ombh * w)
        np.testing.assert_allclose(0., otmb, atol=0.26)

        mb_mod = massbalance.TodayMassBalanceModel(gdir)

        ofl = gdir.read_pickle('inversion_flowlines', div_id=0)
        h = np.array([])
        w = np.array([])
        for fl in ofl:
            h = np.append(h, fl.surface_h)
            w = np.append(w, fl.widths)

        mbh = mb_mod.get_mb(h, None)
        tmb = np.sum(mbh * w)
        self.assertTrue(tmb < otmb)

        if do_plot:  # pragma: no cover
            plt.plot(h, ombh, 'o')
            plt.plot(h, mbh, 'o')
            plt.show()

    def test_backwards_mb(self):

        gdir = init_hef(border=DOM_BORDER)
        flowline.init_present_time_glacier(gdir)

        mb_mod_ref = massbalance.TstarMassBalanceModel(gdir)
        mb_mod = massbalance.BackwardsMassBalanceModel(gdir, use_tstar=True)

        ofl = gdir.read_pickle('inversion_flowlines', div_id=0)
        h = np.array([])
        w = np.array([])
        for fl in ofl:
            h = np.append(h, fl.surface_h)
            w = np.append(w, fl.widths)

        ombh = mb_mod_ref.get_mb(h, None) * SEC_IN_YEAR
        mbh = mb_mod.get_mb(h, None) * SEC_IN_YEAR
        mb_mod.set_bias(100.)
        mbhb = mb_mod.get_mb(h, None) * SEC_IN_YEAR

        np.testing.assert_allclose(ombh, mbh, rtol=0.001)
        self.assertTrue(np.mean(mbhb) > np.mean(mbh))

        if do_plot:  # pragma: no cover
            plt.plot(h, ombh, 'o')
            plt.plot(h, mbh, 'x')
            plt.plot(h, mbhb, 'x')
            plt.show()

        mb_mod_ref = massbalance.TodayMassBalanceModel(gdir)
        mb_mod = massbalance.BackwardsMassBalanceModel(gdir)

        ofl = gdir.read_pickle('inversion_flowlines', div_id=0)
        h = np.array([])
        w = np.array([])
        for fl in ofl:
            h = np.append(h, fl.surface_h)
            w = np.append(w, fl.widths)

        ombh = mb_mod_ref.get_mb(h, None) * SEC_IN_YEAR
        mbh = mb_mod.get_mb(h, None) * SEC_IN_YEAR
        mb_mod.set_bias(-100.)
        mbhb = mb_mod.get_mb(h, None) * SEC_IN_YEAR

        np.testing.assert_allclose(ombh, mbh, rtol=0.005)
        self.assertTrue(np.mean(mbhb) < np.mean(mbh))

        if do_plot:  # pragma: no cover
            plt.plot(h, ombh, 'o')
            plt.plot(h, mbh, 'x')
            plt.plot(h, mbhb, 'x')
            plt.show()


    def test_histalp_mb(self):

        gdir = init_hef(border=DOM_BORDER)
        flowline.init_present_time_glacier(gdir)

        mb_mod = massbalance.TodayMassBalanceModel(gdir)

        ofl = gdir.read_pickle('inversion_flowlines', div_id=0)
        h = np.array([])
        w = np.array([])
        for fl in ofl:
            h = np.append(h, fl.surface_h)
            w = np.append(w, fl.widths)

        ref_mbh = mb_mod.get_mb(h, None) * 365 * 24 * 3600

        mb_mod = massbalance.HistalpMassBalanceModel(gdir)

        # Climate period
        yrs = np.arange(1983, 2003.1, 1)

        my_mb = ref_mbh * 0.
        for yr in yrs:
            my_mb += mb_mod.get_mb(h, yr)
        my_mb = my_mb / len(yrs) * 365 * 24 * 3600

        np.testing.assert_allclose(ref_mbh, my_mb, rtol=0.01)


class TestIdealisedCases(unittest.TestCase):

    def setUp(self):
        self.glen_a = 2.4e-24    # Modern style Glen parameter A
        self.aglen_old = (N+2) * 1.9e-24 / 2. # outdated value
        self.fd = 2. * self.glen_a / (N + 2.)  # equivalent to glen_a
        self.fs = 0             # set slidin
        self.fs_old = 5.7e-20  # outdated value
        
    def tearDown(self):
        pass

    def test_constant_bed(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]

        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 700, 2)
        for model in models:
            fls = dummy_constant_bed()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, y0=0., glen_a=self.glen_a,
                          fs=self.fs, fixed_dt=10*SEC_IN_DAY)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        if do_plot:  # pragma: no cover
            plt.figure()
            plt.plot(yrs, lens[0], 'r')
            plt.plot(yrs, lens[1], 'b')
            plt.plot(yrs, lens[2], 'g')
            plt.title('Compare Length')
            plt.xlabel('years')
            plt.ylabel('[m]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(yrs, volume[0], 'r')
            plt.plot(yrs, volume[1], 'b')
            plt.plot(yrs, volume[2], 'g')
            plt.title('Compare Volume')
            plt.xlabel('years')
            plt.ylabel('[km^3]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(fls[-1].bed_h, 'k')
            plt.plot(surface_h[0], 'r')
            plt.plot(surface_h[1], 'b')
            plt.plot(surface_h[2], 'g')
            plt.title('Compare Shape')
            plt.xlabel('[m]')
            plt.ylabel('Elevation [m]')
            plt.legend(['Bed','Karthaus','Flux','MUSCL-SuperBee'],loc=3)
            plt.show()

        np.testing.assert_almost_equal(lens[0][-1], lens[1][-1])
        np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=3e-3)
        np.testing.assert_allclose(volume[1][-1], volume[2][-1], atol=3e-3)

        self.assertTrue(utils.rmsd(lens[0], lens[2])<50.)
        self.assertTrue(utils.rmsd(lens[1], lens[2])<50.)
        self.assertTrue(utils.rmsd(volume[0], volume[2])<2e-3)
        self.assertTrue(utils.rmsd(volume[1], volume[2])<2e-3)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[2])<1.0)
        self.assertTrue(utils.rmsd(surface_h[1], surface_h[2])<1.0)

    def test_constant_bed_cliff(self):
        """ a test case for mass conservation in the flowline models
            the idea is to introduce a cliff in the sloping bed and see
            what the models do when the cliff height is changed
        """
        
        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]

        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 700, 2)
        for model in models:
            fls = dummy_constant_bed_cliff()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, y0=0., glen_a=self.glen_a,
                          fs=self.fs, fixed_dt=2*SEC_IN_DAY)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        if do_plot:  # pragma: no cover
            plt.figure()
            plt.plot(yrs, lens[0], 'r')
            plt.plot(yrs, lens[1], 'b')
            plt.plot(yrs, lens[2], 'g')
            plt.title('Compare Length')
            plt.xlabel('years')
            plt.ylabel('[m]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(yrs, volume[0], 'r')
            plt.plot(yrs, volume[1], 'b')
            plt.plot(yrs, volume[2], 'g')
            plt.title('Compare Volume')
            plt.xlabel('years')
            plt.ylabel('[km^3]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(fls[-1].bed_h, 'k')
            plt.plot(surface_h[0], 'r')
            plt.plot(surface_h[1], 'b')
            plt.plot(surface_h[2], 'g')
            plt.title('Compare Shape')
            plt.xlabel('[m]')
            plt.ylabel('Elevation [m]')
            plt.legend(['Bed','Karthaus','Flux','MUSCL-SuperBee'],loc=3)
            plt.show()

        # OK, so basically, Alex's tests below show that the other models
        # are wrong and produce too much mass. There is also another more
        # more trivial issue with the computation of the length, I added a
        # "to do" in the code.

        # Unit-testing perspective:
        # verify that indeed the models are wrong of more than 50%
        self.assertTrue(volume[1][-1] > volume[2][-1] * 1.5)
        # Karthaus is even worse
        self.assertTrue(volume[0][-1] > volume[1][-1])

        if False:
            # TODO: this will always fail so ignore it for now
            np.testing.assert_almost_equal(lens[0][-1], lens[1][-1])
            np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=2e-3)
            np.testing.assert_allclose(volume[1][-1], volume[2][-1], atol=2e-3)

            self.assertTrue(utils.rmsd(lens[0], lens[2])<50.)
            self.assertTrue(utils.rmsd(lens[1], lens[2])<50.)
            self.assertTrue(utils.rmsd(volume[0], volume[2])<1e-3)
            self.assertTrue(utils.rmsd(volume[1], volume[2])<1e-3)
            self.assertTrue(utils.rmsd(surface_h[0], surface_h[2])<1.0)
            self.assertTrue(utils.rmsd(surface_h[1], surface_h[2])<1.0)
            
    @requires_working_conda
    def test_equilibrium(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel]

        vols = []
        for model in models:
            fls = dummy_constant_bed()
            mb = ConstantBalanceModel(2600.)
            model = model(fls, mb_model=mb, glen_a=self.glen_a,
                          fixed_dt=10 * SEC_IN_DAY)

            model.run_until_equilibrium()
            vols.append(model.volume_km3)

        ref_vols = []
        for model in models:
            fls = dummy_constant_bed()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, glen_a=self.glen_a,
                          fixed_dt=10 * SEC_IN_DAY)

            model.run_until(600)
            ref_vols.append(model.volume_km3)

        np.testing.assert_allclose(ref_vols, vols, atol=0.005)

    def test_adaptive_ts(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]
        steps = [SEC_IN_MONTH, None, None]
        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 500, 2)
        for model, step in zip(models, steps):
            fls = dummy_constant_bed()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, glen_a=self.glen_a, fixed_dt=step)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        np.testing.assert_almost_equal(lens[0][-1], lens[1][-1])
        np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2)
        np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2)

        self.assertTrue(utils.rmsd(lens[0], lens[1])<50.)
        self.assertTrue(utils.rmsd(volume[0], volume[1])<1e-3)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[1]) < 5)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[2]) < 5)

    def test_bumpy_bed(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]
        steps = [15 * SEC_IN_DAY, None, None]
        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 500, 2)
        for model, step in zip(models, steps):
            fls = dummy_bumpy_bed()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, glen_a=self.glen_a, fixed_dt=step)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        if do_plot:  # pragma: no cover
            plt.figure()
            plt.plot(yrs, lens[0], 'r')
            plt.plot(yrs, lens[1], 'b')
            plt.plot(yrs, lens[2], 'g')
            plt.title('Compare Length')
            plt.xlabel('years')
            plt.ylabel('[m]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(yrs, volume[0], 'r')
            plt.plot(yrs, volume[1], 'b')
            plt.plot(yrs, volume[2], 'g')
            plt.title('Compare Volume')
            plt.xlabel('years')
            plt.ylabel('[km^3]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(fls[-1].bed_h, 'k')
            plt.plot(surface_h[0], 'r')
            plt.plot(surface_h[1], 'b')
            plt.plot(surface_h[2], 'g')
            plt.title('Compare Shape')
            plt.xlabel('[m]')
            plt.ylabel('Elevation [m]')
            plt.legend(['Bed','Karthaus','Flux','MUSCL-SuperBee'],loc=3)
            plt.show()

        np.testing.assert_almost_equal(lens[0][-1], lens[1][-1])
        np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2)
        np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2)

        self.assertTrue(utils.rmsd(lens[0], lens[1])<50.)
        self.assertTrue(utils.rmsd(volume[0], volume[1])<1e-2)
        self.assertTrue(utils.rmsd(volume[0], volume[2])<1e-2)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[1])<5)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[2])<5)

    def test_noisy_bed(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]
        steps = [15 * SEC_IN_DAY, None, None]
        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 500, 2)
        fls_orig = dummy_noisy_bed()
        for model, step in zip(models, steps):
            fls = copy.deepcopy(fls_orig)
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, glen_a=self.glen_a, fixed_dt=step)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        if do_plot:  # pragma: no cover
            plt.figure()
            plt.plot(yrs, lens[0], 'r')
            plt.plot(yrs, lens[1], 'b')
            plt.plot(yrs, lens[2], 'g')
            plt.title('Compare Length')
            plt.xlabel('years')
            plt.ylabel('[m]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(yrs, volume[0], 'r')
            plt.plot(yrs, volume[1], 'b')
            plt.plot(yrs, volume[2], 'g')
            plt.title('Compare Volume')
            plt.xlabel('years')
            plt.ylabel('[km^3]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(fls[-1].bed_h, 'k')
            plt.plot(surface_h[0], 'r')
            plt.plot(surface_h[1], 'b')
            plt.plot(surface_h[2], 'g')
            plt.title('Compare Shape')
            plt.xlabel('[m]')
            plt.ylabel('Elevation [m]')
            plt.legend(['Bed','Karthaus','Flux','MUSCL-SuperBee'],loc=3)
            plt.show()

        np.testing.assert_almost_equal(lens[0][-1], lens[1][-1])
        np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2)
        np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2)

        self.assertTrue(utils.rmsd(lens[0], lens[1])<100.)
        self.assertTrue(utils.rmsd(volume[0], volume[1])<1e-1)
        self.assertTrue(utils.rmsd(volume[0], volume[2])<1e-1)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[1])<10)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[2])<10)

    def test_varying_width(self):
        """This test is for a flowline glacier of variying width, i.e with an
         accumulation area twice as wide as the tongue."""

        # TODO: @alexjarosch here we should have a look at MUSCLSuperBeeModel
        # set do_plot = True to see the plots

        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]
        steps = [15 * SEC_IN_DAY, None, None]
        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 500, 2)
        for model, step in zip(models, steps):
            fls = dummy_width_bed()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, glen_a=self.glen_a, fixed_dt=step)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        if do_plot:  # pragma: no cover
            plt.figure()
            plt.plot(yrs, lens[0], 'r')
            plt.plot(yrs, lens[1], 'b')
            plt.plot(yrs, lens[2], 'g')
            plt.title('Compare Length')
            plt.xlabel('years')
            plt.ylabel('[m]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(yrs, volume[0], 'r')
            plt.plot(yrs, volume[1], 'b')
            plt.plot(yrs, volume[2], 'g')
            plt.title('Compare Volume')
            plt.xlabel('years')
            plt.ylabel('[km^3]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(fls[-1].bed_h, 'k')
            plt.plot(surface_h[0], 'r')
            plt.plot(surface_h[1], 'b')
            plt.plot(surface_h[2], 'g')
            plt.title('Compare Shape')
            plt.xlabel('[m]')
            plt.ylabel('Elevation [m]')
            plt.legend(['Bed','Karthaus','Flux','MUSCL-SuperBee'],loc=3)
            plt.show()

        np.testing.assert_almost_equal(lens[0][-1], lens[1][-1])
        np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2)

        np.testing.assert_allclose(utils.rmsd(lens[0], lens[1]), 0., atol=50)
        np.testing.assert_allclose(utils.rmsd(volume[0], volume[1]), 0.,
                                   atol=3e-3)
        np.testing.assert_allclose(utils.rmsd(surface_h[0], surface_h[1]), 0.,
                                   atol=2)

    def test_tributary(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel]
        steps = [15 * SEC_IN_DAY, None]
        flss = [dummy_width_bed(), dummy_width_bed_tributary()]
        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 500, 2)
        for model, step, fls in zip(models, steps, flss):
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, fs=self.fs_old,
                          glen_a=self.aglen_old,
                          fixed_dt=step)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = np.sum([f.volume_km3 for f in fls])
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())<commit_msg>constant mass-balance units changed, and so did the tests<commit_after>from __future__ import division

import warnings

from six.moves import zip

warnings.filterwarnings("once", category=DeprecationWarning)

import logging
logging.basicConfig(format='%(asctime)s: %(name)s: %(message)s',
                    datefmt='%Y-%m-%d %H:%M:%S', level=logging.DEBUG)

import unittest
import os
import copy

import shapely.geometry as shpg
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Local imports
from oggm.tests.test_graphics import init_hef
from oggm.core.models import massbalance, flowline
from oggm.core.models.massbalance import ConstantBalanceModel
from oggm.tests import is_slow, requires_working_conda
from oggm.tests import HAS_NEW_GDAL

from oggm import utils, cfg
from oggm.cfg import N, SEC_IN_DAY, SEC_IN_MONTH, SEC_IN_YEAR

# Globals
current_dir = os.path.dirname(os.path.abspath(__file__))

# test directory
testdir = os.path.join(current_dir, 'tmp')

do_plot = False

DOM_BORDER = 80

def dummy_constant_bed(hmax=3000., hmin=1000., nx=200):

    map_dx = 100.
    dx = 1.

    surface_h = np.linspace(hmax, hmin, nx)
    bed_h = surface_h
    widths = surface_h * 0. + 3.

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.VerticalWallFlowline(line, dx, map_dx, surface_h,
                                          bed_h, widths)]

def dummy_constant_bed_cliff(hmax=3000., hmin=1000., nx=200):
    """
    I introduce a cliff in the bed to test the mass conservation of the models
    Such a cliff could be real or a DEM error/artifact
    """
    
    map_dx = 100.
    dx = 1.

    surface_h = np.linspace(hmax, hmin, nx)
    
    cliff_height = 250.0
    surface_h[50:] = surface_h[50:] - cliff_height
    
    bed_h = surface_h
    widths = surface_h * 0. + 1.

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.VerticalWallFlowline(line, dx, map_dx, surface_h,
                                          bed_h, widths)]

def dummy_bumpy_bed():

    map_dx = 100.
    dx = 1.
    nx = 200

    coords = np.arange(0, nx-0.5, 1)
    surface_h = np.linspace(3000, 1000, nx)
    surface_h += 170. * np.exp(-((coords-30)/5)**2)

    bed_h = surface_h
    widths = surface_h * 0. + 3.
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.VerticalWallFlowline(line, dx, map_dx, surface_h,
                                          bed_h, widths)]


def dummy_noisy_bed():

    map_dx = 100.
    dx = 1.
    nx = 200
    np.random.seed(42)
    coords = np.arange(0, nx-0.5, 1)
    surface_h = np.linspace(3000, 1000, nx)
    surface_h += 100 * np.random.rand(nx) - 50.

    bed_h = surface_h
    widths = surface_h * 0. + 3.
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.VerticalWallFlowline(line, dx, map_dx, surface_h,
                                          bed_h, widths)]

def dummy_parabolic_bed():

    map_dx = 100.
    dx = 1.
    nx = 200

    surface_h = np.linspace(3000, 1000, nx)
    bed_h = surface_h
    shape = surface_h * 0. + 5.e-03

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.ParabolicFlowline(line, dx, map_dx, surface_h,
                                       bed_h, shape)]

def dummy_mixed_bed():

    map_dx = 100.
    dx = 1.
    nx = 200

    surface_h = np.linspace(3000, 1000, nx)
    bed_h = surface_h
    shape = surface_h * 0. + 3.e-03
    shape[10:20] = 0.00001

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.MixedFlowline(line, dx, map_dx, surface_h,
                                   bed_h, shape, lambdas=3.5)]


def dummy_trapezoidal_bed(hmax=3000., hmin=1000., nx=200):

    map_dx = 100.
    dx = 1.

    surface_h = np.linspace(hmax, hmin, nx)
    bed_h = surface_h
    widths = surface_h * 0. + 1.6

    lambdas = surface_h * 0. + 2

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)

    return [flowline.TrapezoidalFlowline(line, dx, map_dx, surface_h,
                                         bed_h, widths, lambdas)]


def dummy_width_bed():
    """This bed has a width of 6 during the first 20 points and then 3"""

    map_dx = 100.
    dx = 1.
    nx = 200

    surface_h = np.linspace(3000, 1000, nx)
    bed_h = surface_h
    widths = surface_h * 0. + 3.
    widths[0:20] = 6.

    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)
    return [flowline.VerticalWallFlowline(line, dx, map_dx, surface_h,
                                          bed_h, widths)]

def dummy_width_bed_tributary():

    map_dx = 100.
    dx = 1.
    nx = 200

    surface_h = np.linspace(3000, 1000, nx)
    bed_h = surface_h
    widths = surface_h * 0. + 3.
    coords = np.arange(0, nx-0.5, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.]).T)

    fl_0 = flowline.VerticalWallFlowline(line, dx, map_dx, surface_h, bed_h,
                                         widths)
    coords = np.arange(0, 19.1, 1)
    line = shpg.LineString(np.vstack([coords, coords*0.+1]).T)
    fl_1 = flowline.VerticalWallFlowline(line, dx, map_dx, surface_h[0:20],
                                         bed_h[0:20], widths[0:20])
    fl_1.set_flows_to(fl_0)
    return [fl_1, fl_0]

class TestInitFlowline(unittest.TestCase):

    def setUp(self):
        pass

    def tearDown(self):
        pass

    def test_init_present_time_glacier(self):

        gdir = init_hef(border=DOM_BORDER)
        flowline.init_present_time_glacier(gdir)

        fls = gdir.read_pickle('model_flowlines')

        lens = [len(gdir.read_pickle('centerlines', div_id=i)) for i in [1,
                                                                         2, 3]]
        ofl = gdir.read_pickle('inversion_flowlines',
                               div_id=np.argmax(lens)+1)[-1]

        self.assertTrue(gdir.rgi_date.year == 2003)
        self.assertTrue(len(fls) == 5)

        vol = 0.
        area = 0.
        for fl in fls:
            refo = 1 if fl is fls[-1] else 0
            self.assertTrue(fl.order == refo)
            ref = np.arange(len(fl.surface_h)) * fl.dx
            np.testing.assert_allclose(ref, fl.dis_on_line,
                                       rtol=0.001,
                                       atol=0.01)
            self.assertTrue(len(fl.surface_h) ==
                            len(fl.bed_h) ==
                            len(fl.bed_shape) ==
                            len(fl.dis_on_line) ==
                            len(fl.widths))

            self.assertTrue(np.all(fl.bed_shape >= 0))
            self.assertTrue(np.all(fl.widths >= 0))
            vol += fl.volume_km3
            area += fl.area_km2

            if refo == 1:
                rtol = 0.07 if HAS_NEW_GDAL else 0.05
                np.testing.assert_allclose(ofl.widths * gdir.grid.dx,
                                           fl.widths_m[0:len(ofl.widths)],
                                           rtol=rtol)

        np.testing.assert_allclose(0.573, vol, rtol=0.001)
        np.testing.assert_allclose(7350.0, fls[-1].length_m)

        rtol = 0.01 if HAS_NEW_GDAL else 0.001
        np.testing.assert_allclose(gdir.rgi_area_km2, area, rtol=rtol)

        if do_plot:  # pragma: no cover
            plt.plot(fls[-1].bed_h)
            plt.plot(fls[-1].surface_h)
            plt.show()

    def test_present_time_glacier_massbalance(self):

        gdir = init_hef(border=DOM_BORDER)
        flowline.init_present_time_glacier(gdir)

        mb_mod = massbalance.HistalpMassBalanceModel(gdir)

        fls = gdir.read_pickle('model_flowlines')
        glacier = flowline.FlowlineModel(fls)

        hef_file = utils.get_demo_file('mbdata_RGI40-11.00897.csv')
        mbdf = pd.read_csv(hef_file).set_index('YEAR')

        hgts = np.array([])
        widths = np.array([])
        for fl in glacier.fls:
            hgts = np.concatenate((hgts, fl.surface_h))
            widths = np.concatenate((widths, fl.widths_m))
        tot_mb = []
        refmb = []
        grads = hgts * 0
        for yr, mb in mbdf.iterrows():
            refmb.append(mb['ANNUAL_BALANCE'])
            mbh = mb_mod.get_mb(hgts, yr) * SEC_IN_YEAR * cfg.RHO
            grads += mbh
            tot_mb.append(np.average(mbh, weights=widths))
        grads /= len(tot_mb)

        # Bias
        self.assertTrue(np.abs(utils.md(tot_mb, refmb)) < 50)

        # Gradient
        dfg = pd.read_csv(utils.get_demo_file('mbgrads_RGI40-11.00897.csv'),
                          index_col='ALTITUDE').mean(axis=1)

        # Take the altitudes below 3100 and fit a line
        dfg = dfg[dfg.index < 3100]
        pok = np.where(hgts < 3100)
        from scipy.stats import linregress
        slope_obs, _, _, _, _ = linregress(dfg.index, dfg.values)
        slope_our, _, _, _, _ = linregress(hgts[pok], grads[pok])
        np.testing.assert_allclose(slope_obs, slope_our, rtol=0.1)


class TestMassBalance(unittest.TestCase):
    def setUp(self):
        pass

    def tearDown(self):
        pass

    def test_tstar_mb(self):

        gdir = init_hef(border=DOM_BORDER)
        flowline.init_present_time_glacier(gdir)

        mb_mod = massbalance.TstarMassBalanceModel(gdir)

        ofl = gdir.read_pickle('inversion_flowlines', div_id=0)
        h = np.array([])
        w = np.array([])
        for fl in ofl:
            h = np.append(h, fl.surface_h)
            w = np.append(w, fl.widths)

        ombh = mb_mod.get_mb(h, None)
        otmb = np.sum(ombh * w)
        np.testing.assert_allclose(0., otmb, atol=0.26)

        mb_mod = massbalance.TodayMassBalanceModel(gdir)

        ofl = gdir.read_pickle('inversion_flowlines', div_id=0)
        h = np.array([])
        w = np.array([])
        for fl in ofl:
            h = np.append(h, fl.surface_h)
            w = np.append(w, fl.widths)

        mbh = mb_mod.get_mb(h, None)
        tmb = np.sum(mbh * w)
        self.assertTrue(tmb < otmb)

        if do_plot:  # pragma: no cover
            plt.plot(h, ombh, 'o')
            plt.plot(h, mbh, 'o')
            plt.show()

    def test_backwards_mb(self):

        gdir = init_hef(border=DOM_BORDER)
        flowline.init_present_time_glacier(gdir)

        mb_mod_ref = massbalance.TstarMassBalanceModel(gdir)
        mb_mod = massbalance.BackwardsMassBalanceModel(gdir, use_tstar=True)

        ofl = gdir.read_pickle('inversion_flowlines', div_id=0)
        h = np.array([])
        w = np.array([])
        for fl in ofl:
            h = np.append(h, fl.surface_h)
            w = np.append(w, fl.widths)

        ombh = mb_mod_ref.get_mb(h, None) * SEC_IN_YEAR
        mbh = mb_mod.get_mb(h, None) * SEC_IN_YEAR
        mb_mod.set_bias(100.)
        mbhb = mb_mod.get_mb(h, None) * SEC_IN_YEAR

        np.testing.assert_allclose(ombh, mbh, rtol=0.001)
        self.assertTrue(np.mean(mbhb) > np.mean(mbh))

        if do_plot:  # pragma: no cover
            plt.plot(h, ombh, 'o')
            plt.plot(h, mbh, 'x')
            plt.plot(h, mbhb, 'x')
            plt.show()

        mb_mod_ref = massbalance.TodayMassBalanceModel(gdir)
        mb_mod = massbalance.BackwardsMassBalanceModel(gdir)

        ofl = gdir.read_pickle('inversion_flowlines', div_id=0)
        h = np.array([])
        w = np.array([])
        for fl in ofl:
            h = np.append(h, fl.surface_h)
            w = np.append(w, fl.widths)

        ombh = mb_mod_ref.get_mb(h, None) * SEC_IN_YEAR
        mbh = mb_mod.get_mb(h, None) * SEC_IN_YEAR
        mb_mod.set_bias(-100.)
        mbhb = mb_mod.get_mb(h, None) * SEC_IN_YEAR

        np.testing.assert_allclose(ombh, mbh, rtol=0.005)
        self.assertTrue(np.mean(mbhb) < np.mean(mbh))

        if do_plot:  # pragma: no cover
            plt.plot(h, ombh, 'o')
            plt.plot(h, mbh, 'x')
            plt.plot(h, mbhb, 'x')
            plt.show()


    def test_histalp_mb(self):

        gdir = init_hef(border=DOM_BORDER)
        flowline.init_present_time_glacier(gdir)

        mb_mod = massbalance.TodayMassBalanceModel(gdir)

        ofl = gdir.read_pickle('inversion_flowlines', div_id=0)
        h = np.array([])
        w = np.array([])
        for fl in ofl:
            h = np.append(h, fl.surface_h)
            w = np.append(w, fl.widths)

        ref_mbh = mb_mod.get_mb(h, None) * 365 * 24 * 3600

        mb_mod = massbalance.HistalpMassBalanceModel(gdir)

        # Climate period
        yrs = np.arange(1983, 2003.1, 1)

        my_mb = ref_mbh * 0.
        for yr in yrs:
            my_mb += mb_mod.get_mb(h, yr)
        my_mb = my_mb / len(yrs) * 365 * 24 * 3600

        np.testing.assert_allclose(ref_mbh, my_mb, rtol=0.01)


class TestIdealisedCases(unittest.TestCase):

    def setUp(self):
        self.glen_a = 2.4e-24    # Modern style Glen parameter A
        self.aglen_old = (N+2) * 1.9e-24 / 2. # outdated value
        self.fd = 2. * self.glen_a / (N + 2.)  # equivalent to glen_a
        self.fs = 0             # set slidin
        self.fs_old = 5.7e-20  # outdated value
        
    def tearDown(self):
        pass

    def test_constant_bed(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]

        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 700, 2)
        for model in models:
            fls = dummy_constant_bed()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, y0=0., glen_a=self.glen_a,
                          fs=self.fs, fixed_dt=10*SEC_IN_DAY)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        if do_plot:  # pragma: no cover
            plt.figure()
            plt.plot(yrs, lens[0], 'r')
            plt.plot(yrs, lens[1], 'b')
            plt.plot(yrs, lens[2], 'g')
            plt.title('Compare Length')
            plt.xlabel('years')
            plt.ylabel('[m]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(yrs, volume[0], 'r')
            plt.plot(yrs, volume[1], 'b')
            plt.plot(yrs, volume[2], 'g')
            plt.title('Compare Volume')
            plt.xlabel('years')
            plt.ylabel('[km^3]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(fls[-1].bed_h, 'k')
            plt.plot(surface_h[0], 'r')
            plt.plot(surface_h[1], 'b')
            plt.plot(surface_h[2], 'g')
            plt.title('Compare Shape')
            plt.xlabel('[m]')
            plt.ylabel('Elevation [m]')
            plt.legend(['Bed','Karthaus','Flux','MUSCL-SuperBee'],loc=3)
            plt.show()

        np.testing.assert_almost_equal(lens[0][-1], lens[1][-1])
        np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=3e-3)
        np.testing.assert_allclose(volume[1][-1], volume[2][-1], atol=3e-3)

        self.assertTrue(utils.rmsd(lens[0], lens[2])<50.)
        self.assertTrue(utils.rmsd(lens[1], lens[2])<50.)
        self.assertTrue(utils.rmsd(volume[0], volume[2])<2e-3)
        self.assertTrue(utils.rmsd(volume[1], volume[2])<2e-3)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[2])<1.0)
        self.assertTrue(utils.rmsd(surface_h[1], surface_h[2])<1.0)

    def test_constant_bed_cliff(self):
        """ a test case for mass conservation in the flowline models
            the idea is to introduce a cliff in the sloping bed and see
            what the models do when the cliff height is changed
        """
        
        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]

        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 700, 2)
        for model in models:
            fls = dummy_constant_bed_cliff()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, y0=0., glen_a=self.glen_a,
                          fs=self.fs, fixed_dt=2*SEC_IN_DAY)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        if do_plot:  # pragma: no cover
            plt.figure()
            plt.plot(yrs, lens[0], 'r')
            plt.plot(yrs, lens[1], 'b')
            plt.plot(yrs, lens[2], 'g')
            plt.title('Compare Length')
            plt.xlabel('years')
            plt.ylabel('[m]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(yrs, volume[0], 'r')
            plt.plot(yrs, volume[1], 'b')
            plt.plot(yrs, volume[2], 'g')
            plt.title('Compare Volume')
            plt.xlabel('years')
            plt.ylabel('[km^3]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(fls[-1].bed_h, 'k')
            plt.plot(surface_h[0], 'r')
            plt.plot(surface_h[1], 'b')
            plt.plot(surface_h[2], 'g')
            plt.title('Compare Shape')
            plt.xlabel('[m]')
            plt.ylabel('Elevation [m]')
            plt.legend(['Bed','Karthaus','Flux','MUSCL-SuperBee'],loc=3)
            plt.show()

        # OK, so basically, Alex's tests below show that the other models
        # are wrong and produce too much mass. There is also another more
        # more trivial issue with the computation of the length, I added a
        # "to do" in the code.

        # Unit-testing perspective:
        # verify that indeed the models are wrong of more than 50%
        self.assertTrue(volume[1][-1] > volume[2][-1] * 1.5)
        # Karthaus is even worse
        self.assertTrue(volume[0][-1] > volume[1][-1])

        if False:
            # TODO: this will always fail so ignore it for now
            np.testing.assert_almost_equal(lens[0][-1], lens[1][-1])
            np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=2e-3)
            np.testing.assert_allclose(volume[1][-1], volume[2][-1], atol=2e-3)

            self.assertTrue(utils.rmsd(lens[0], lens[2])<50.)
            self.assertTrue(utils.rmsd(lens[1], lens[2])<50.)
            self.assertTrue(utils.rmsd(volume[0], volume[2])<1e-3)
            self.assertTrue(utils.rmsd(volume[1], volume[2])<1e-3)
            self.assertTrue(utils.rmsd(surface_h[0], surface_h[2])<1.0)
            self.assertTrue(utils.rmsd(surface_h[1], surface_h[2])<1.0)
            
    @requires_working_conda
    def test_equilibrium(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel]

        vols = []
        for model in models:
            fls = dummy_constant_bed()
            mb = ConstantBalanceModel(2600.)
            model = model(fls, mb_model=mb, glen_a=self.glen_a,
                          fixed_dt=10 * SEC_IN_DAY)

            model.run_until_equilibrium()
            vols.append(model.volume_km3)

        ref_vols = []
        for model in models:
            fls = dummy_constant_bed()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, glen_a=self.glen_a,
                          fixed_dt=10 * SEC_IN_DAY)

            model.run_until(600)
            ref_vols.append(model.volume_km3)

        np.testing.assert_allclose(ref_vols, vols, atol=0.01)

    def test_adaptive_ts(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]
        steps = [SEC_IN_MONTH, None, None]
        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 500, 2)
        for model, step in zip(models, steps):
            fls = dummy_constant_bed()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, glen_a=self.glen_a, fixed_dt=step)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        np.testing.assert_almost_equal(lens[0][-1], lens[1][-1])
        np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2)
        np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2)

        self.assertTrue(utils.rmsd(lens[0], lens[1])<50.)
        self.assertTrue(utils.rmsd(volume[0], volume[1])<1e-3)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[1]) < 5)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[2]) < 5)

    def test_bumpy_bed(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]
        steps = [15 * SEC_IN_DAY, None, None]
        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 500, 2)
        for model, step in zip(models, steps):
            fls = dummy_bumpy_bed()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, glen_a=self.glen_a, fixed_dt=step)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        if do_plot:  # pragma: no cover
            plt.figure()
            plt.plot(yrs, lens[0], 'r')
            plt.plot(yrs, lens[1], 'b')
            plt.plot(yrs, lens[2], 'g')
            plt.title('Compare Length')
            plt.xlabel('years')
            plt.ylabel('[m]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(yrs, volume[0], 'r')
            plt.plot(yrs, volume[1], 'b')
            plt.plot(yrs, volume[2], 'g')
            plt.title('Compare Volume')
            plt.xlabel('years')
            plt.ylabel('[km^3]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(fls[-1].bed_h, 'k')
            plt.plot(surface_h[0], 'r')
            plt.plot(surface_h[1], 'b')
            plt.plot(surface_h[2], 'g')
            plt.title('Compare Shape')
            plt.xlabel('[m]')
            plt.ylabel('Elevation [m]')
            plt.legend(['Bed','Karthaus','Flux','MUSCL-SuperBee'],loc=3)
            plt.show()

        np.testing.assert_almost_equal(lens[0][-1], lens[1][-1])
        np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2)
        np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2)

        self.assertTrue(utils.rmsd(lens[0], lens[1])<50.)
        self.assertTrue(utils.rmsd(volume[0], volume[1])<1e-2)
        self.assertTrue(utils.rmsd(volume[0], volume[2])<1e-2)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[1])<5)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[2])<5)

    def test_noisy_bed(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]
        steps = [15 * SEC_IN_DAY, None, None]
        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 500, 2)
        fls_orig = dummy_noisy_bed()
        for model, step in zip(models, steps):
            fls = copy.deepcopy(fls_orig)
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, glen_a=self.glen_a, fixed_dt=step)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        if do_plot:  # pragma: no cover
            plt.figure()
            plt.plot(yrs, lens[0], 'r')
            plt.plot(yrs, lens[1], 'b')
            plt.plot(yrs, lens[2], 'g')
            plt.title('Compare Length')
            plt.xlabel('years')
            plt.ylabel('[m]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(yrs, volume[0], 'r')
            plt.plot(yrs, volume[1], 'b')
            plt.plot(yrs, volume[2], 'g')
            plt.title('Compare Volume')
            plt.xlabel('years')
            plt.ylabel('[km^3]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(fls[-1].bed_h, 'k')
            plt.plot(surface_h[0], 'r')
            plt.plot(surface_h[1], 'b')
            plt.plot(surface_h[2], 'g')
            plt.title('Compare Shape')
            plt.xlabel('[m]')
            plt.ylabel('Elevation [m]')
            plt.legend(['Bed','Karthaus','Flux','MUSCL-SuperBee'],loc=3)
            plt.show()

        np.testing.assert_allclose(lens[0][-1], lens[1][-1], atol=101)
        np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2)
        np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2)

        self.assertTrue(utils.rmsd(lens[0], lens[1])<100.)
        self.assertTrue(utils.rmsd(volume[0], volume[1])<1e-1)
        self.assertTrue(utils.rmsd(volume[0], volume[2])<1e-1)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[1])<10)
        self.assertTrue(utils.rmsd(surface_h[0], surface_h[2])<10)

    def test_varying_width(self):
        """This test is for a flowline glacier of variying width, i.e with an
         accumulation area twice as wide as the tongue."""

        # TODO: @alexjarosch here we should have a look at MUSCLSuperBeeModel
        # set do_plot = True to see the plots

        models = [flowline.KarthausModel, flowline.FluxBasedModel,
                  flowline.MUSCLSuperBeeModel]
        steps = [15 * SEC_IN_DAY, None, None]
        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 500, 2)
        for model, step in zip(models, steps):
            fls = dummy_width_bed()
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, glen_a=self.glen_a, fixed_dt=step)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = fls[-1].volume_km3
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())

        if do_plot:  # pragma: no cover
            plt.figure()
            plt.plot(yrs, lens[0], 'r')
            plt.plot(yrs, lens[1], 'b')
            plt.plot(yrs, lens[2], 'g')
            plt.title('Compare Length')
            plt.xlabel('years')
            plt.ylabel('[m]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(yrs, volume[0], 'r')
            plt.plot(yrs, volume[1], 'b')
            plt.plot(yrs, volume[2], 'g')
            plt.title('Compare Volume')
            plt.xlabel('years')
            plt.ylabel('[km^3]')
            plt.legend(['Karthaus','Flux','MUSCL-SuperBee'],loc=2)

            plt.figure()
            plt.plot(fls[-1].bed_h, 'k')
            plt.plot(surface_h[0], 'r')
            plt.plot(surface_h[1], 'b')
            plt.plot(surface_h[2], 'g')
            plt.title('Compare Shape')
            plt.xlabel('[m]')
            plt.ylabel('Elevation [m]')
            plt.legend(['Bed','Karthaus','Flux','MUSCL-SuperBee'],loc=3)
            plt.show()

        np.testing.assert_almost_equal(lens[0][-1], lens[1][-1])
        np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2)

        np.testing.assert_allclose(utils.rmsd(lens[0], lens[1]), 0., atol=50)
        np.testing.assert_allclose(utils.rmsd(volume[0], volume[1]), 0.,
                                   atol=3e-3)
        np.testing.assert_allclose(utils.rmsd(surface_h[0], surface_h[1]), 0.,
                                   atol=5)

    def test_tributary(self):

        models = [flowline.KarthausModel, flowline.FluxBasedModel]
        steps = [15 * SEC_IN_DAY, None]
        flss = [dummy_width_bed(), dummy_width_bed_tributary()]
        lens = []
        surface_h = []
        volume = []
        yrs = np.arange(1, 500, 2)
        for model, step, fls in zip(models, steps, flss):
            mb = ConstantBalanceModel(2600.)

            model = model(fls, mb_model=mb, fs=self.fs_old,
                          glen_a=self.aglen_old,
                          fixed_dt=step)

            length = yrs * 0.
            vol = yrs * 0.
            for i, y in enumerate(yrs):
                model.run_until(y)
                length[i] = fls[-1].length_m
                vol[i] = np.sum([f.volume_km3 for f in fls])
            lens.append(length)
            volume.append(vol)
            surface_h.append(fls[-1].surface_h.copy())
{'fraction_non_alphanumeric': 0.13934220015765728, 'fraction_numerical': 0.029015648323633497, 'mean_word_length': 1.8970409601088383, 'pattern_counts': {'":': 0, '<': 65, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 58, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4195217', 'n_tokens_mistral': 22290, 'n_tokens_neox': 20913, 'n_words': 3712}