core.py 65.9 KB
Newer Older
Sebastian Heimann's avatar
Sebastian Heimann committed
2001
2002
2003
2004
          rundir=rundir,
          status=status,
          xs_inject=xs_inject,
          **config.solver_config.get_solver_kwargs())
Sebastian Heimann's avatar
Sebastian Heimann committed
2005

Sebastian Heimann's avatar
Sebastian Heimann committed
2006
    harvest(rundir, problem)
Sebastian Heimann's avatar
Sebastian Heimann committed
2007

Sebastian Heimann's avatar
Sebastian Heimann committed
2008
2009
2010
    tstop = time.time()
    logger.info(
        'stop %i / %i (%g min)' % (ievent, nevents, (tstop - tstart)/60.))
Sebastian Heimann's avatar
Sebastian Heimann committed
2011

2012
2013
2014
    logger.info(
        'done with problem %s, rundir is %s' % (problem.name, rundir))

Sebastian Heimann's avatar
Sebastian Heimann committed
2015

Sebastian Heimann's avatar
Sebastian Heimann committed
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
class ParameterStats(Object):
    name = String.T()
    mean = Float.T()
    std = Float.T()
    best = Float.T()
    minimum = Float.T()
    percentile5 = Float.T()
    percentile16 = Float.T()
    median = Float.T()
    percentile84 = Float.T()
    percentile95 = Float.T()
    maximum = Float.T()
Sebastian Heimann's avatar
Sebastian Heimann committed
2028

Sebastian Heimann's avatar
Sebastian Heimann committed
2029
2030
2031
    def __init__(self, *args, **kwargs):
        kwargs.update(zip(self.T.propnames, args))
        Object.__init__(self, **kwargs)
Sebastian Heimann's avatar
Sebastian Heimann committed
2032
2033


Sebastian Heimann's avatar
Sebastian Heimann committed
2034
2035
2036
2037
2038
2039
class ResultStats(Object):
    problem = Problem.T()
    parameter_stats_list = List.T(ParameterStats.T())


def make_stats(problem, xs, misfits, pnames=None):
Sebastian Heimann's avatar
Sebastian Heimann committed
2040
    gms = problem.global_misfits(misfits)
Sebastian Heimann's avatar
Sebastian Heimann committed
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
    ibest = num.argmin(gms)
    rs = ResultStats(problem=problem)
    if pnames is None:
        pnames = problem.parameter_names

    for pname in pnames:
        iparam = problem.name_to_index(pname)
        vs = problem.extract(xs, iparam)
        mi, p5, p16, median, p84, p95, ma = map(float, num.percentile(
            vs, [0., 5., 16., 50., 84., 95., 100.]))

        mean = float(num.mean(vs))
        std = float(num.std(vs))
        best = float(vs[ibest])
        s = ParameterStats(
            pname, mean, std, best, mi, p5, p16, median, p84, p95, ma)

        rs.parameter_stats_list.append(s)

    return rs
Sebastian Heimann's avatar
Sebastian Heimann committed
2061
2062


Sebastian Heimann's avatar
Sebastian Heimann committed
2063
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
def format_stats(rs, fmt):
    pname_to_pindex = dict(
        (p.name, i) for (i, p) in enumerate(rs.parameter_stats_list))

    values = []
    headers = []
    for x in fmt:
        pname, qname = x.split('.')
        pindex = pname_to_pindex[pname]
        values.append(getattr(rs.parameter_stats_list[pindex], qname))
        headers.append(x)

    return ' '.join('%16.7g' % v for v in values)


def export(what, rundirs, type=None, pnames=None, filename=None):
    if pnames is not None:
        pnames_clean = [pname.split('.')[0] for pname in pnames]
        shortform = all(len(pname.split('.')) == 2 for pname in pnames)
    else:
        pnames_clean = None
        shortform = False

    if what == 'stats' and type is not None:
        raise GrondError('invalid argument combination: what=%s, type=%s' % (
            repr(what), repr(type)))

    if what != 'stats' and shortform:
        raise GrondError('invalid argument combination: what=%s, pnames=%s' % (
            repr(what), repr(pnames)))

    if what != 'stats' and type != 'vector' and pnames is not None:
        raise GrondError(
            'invalid argument combination: what=%s, type=%s, pnames=%s' % (
                repr(what), repr(type), repr(pnames)))

    if filename is None:
        out = sys.stdout
    else:
        out = open(filename, 'w')

    if type is None:
        type = 'event'

    if shortform:
        print >>out, '#', ' '.join('%16s' % x for x in pnames)

2110
    def dump(x, gm, indices):
Sebastian Heimann's avatar
Sebastian Heimann committed
2111
        if type == 'vector':
2112
2113
            print >>out, ' ', ' '.join('%16.7g' % v for v in x[indices]), \
                '%16.7g' % gm
Sebastian Heimann's avatar
Sebastian Heimann committed
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137

        elif type == 'source':
            source = problem.unpack(x)
            guts.dump(source, stream=out)

        elif type == 'event':
            ev = problem.unpack(x).pyrocko_event()
            model.dump_events([ev], stream=out)

        else:
            raise GrondError('invalid argument: type=%s' % repr(type))

    header = None
    for rundir in rundirs:
        problem, xs, misfits = load_problem_info_and_data(
            rundir, subset='harvest')

        if type == 'vector':
            pnames_take = pnames_clean or \
                problem.parameter_names[:problem.nparameters]

            indices = num.array(
                [problem.name_to_index(pname) for pname in pnames_take])

2138
2139
2140
2141
2142
2143
2144
2145
            if type == 'vector' and what in ('best', 'mean', 'ensemble'):
                extra = ['global_misfit']
            else:
                extra = []

            new_header = '# ' + ' '.join(
                '%16s' % x for x in pnames_take + extra)

Sebastian Heimann's avatar
Sebastian Heimann committed
2146
2147
2148
2149
2150
2151
2152
2153
            if type == 'vector' and header != new_header:
                print >>out, new_header

            header = new_header
        else:
            indices = None

        if what == 'best':
2154
2155
            x_best, gm_best = get_best_x_and_gm(problem, xs, misfits)
            dump(x_best, gm_best, indices)
Sebastian Heimann's avatar
Sebastian Heimann committed
2156
2157

        elif what == 'mean':
2158
2159
            x_mean, gm_mean = get_mean_x_and_gm(problem, xs, misfits)
            dump(x_mean, gm_mean, indices)
Sebastian Heimann's avatar
Sebastian Heimann committed
2160
2161
2162
2163
2164

        elif what == 'ensemble':
            gms = problem.global_misfits(misfits)
            isort = num.argsort(gms)
            for i in isort:
2165
                dump(xs[i], gms[i], indices)
Sebastian Heimann's avatar
Sebastian Heimann committed
2166
2167
2168
2169
2170
2171
2172
2173
2174
2175

        elif what == 'stats':
            rs = make_stats(problem, xs, misfits, pnames_clean)
            if shortform:
                print >>out, ' ', format_stats(rs, pnames)
            else:
                print >>out, rs

        else:
            raise GrondError('invalid argument: what=%s' % repr(what))
Sebastian Heimann's avatar
Sebastian Heimann committed
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195

    if out is not sys.stdout:
        out.close()


__all__ = '''
    GrondError
    Parameter
    ADict
    Path
    Problem
    ProblemConfig
    MisfitTarget
    MisfitResult
    Forbidden
    InnerMisfitConfig
    DatasetConfig
    TargetConfig
    SamplerDistributionChoice
    SolverConfig
Sebastian Heimann's avatar
Sebastian Heimann committed
2196
    EngineConfig
Sebastian Heimann's avatar
Sebastian Heimann committed
2197
2198
2199
2200
2201
2202
2203
2204
    Config
    HasPaths
    TargetAnalysisResult
    load_problem_info_and_data
    read_config
    forward
    harvest
    go
Sebastian Heimann's avatar
Sebastian Heimann committed
2205
    get_event_names
Sebastian Heimann's avatar
Sebastian Heimann committed
2206
    check
Sebastian Heimann's avatar
Sebastian Heimann committed
2207
2208
    export
'''.split()