Python scipy.stats.gmean() Examples

def geometric_mean_var(z):
    for row in np.eye(z.shape[1]):
        if not np.any(np.all(row == z, axis=1)):
            z = np.row_stack([z, row])
    n_points, n_dim = z.shape

    D = vectorized_cdist(z, z)
    np.fill_diagonal(D, np.inf)

    k = n_dim - 1
    I = D.argsort(axis=1)[:, :k]

    first = np.column_stack([np.arange(n_points) for _ in range(k)])

    val = gmean(D[first, I], axis=1)

    return val.var() 

def get_statistical_property_scores_per_input_per_impl(self, func: StatisticalPropertyFunc, input_num: int,
                                                           reduce: ReduceFunc = stats.gmean) -> t.Dict[str, t.List[float]]:
        """
        Assumptions:
            - Most programs have the same number of input (known as max input number)
            - The input number n takes roughly the same amount of time for every program category
        """
        cats = self._get_categories_for_number_of_inputs(self.get_max_input_num())
        scores_per_impl = InsertionTimeOrderedDict()
        for cat in cats:
            scores = cat.get_statistical_property_scores_per_input_per_impl(func, cat.get_input_strs()[input_num])
            for impl in scores:
                if impl not in scores_per_impl:
                    scores_per_impl[impl] = []
                scores_per_impl[impl].append(reduce(scores[impl]))
        return scores_per_impl 

def rl_reward(env):

    delay = np.asarray(env.env_D)
    mask = delay == np.inf
    delay[mask] = len(delay)*np.max(delay[~mask])

    if env.PRAEMIUM == 'AVG':
        reward = -np.mean(matrix_to_rl(delay))
    elif env.PRAEMIUM == 'MAX':
        reward = -np.max(matrix_to_rl(delay))
    elif env.PRAEMIUM == 'AXM':
        reward = -(np.mean(matrix_to_rl(delay)) + np.max(matrix_to_rl(delay)))/2
    elif env.PRAEMIUM == 'GEO':
        reward = -stats.gmean(matrix_to_rl(delay))
    elif env.PRAEMIUM == 'LOST':
        reward = -env.env_L
    return reward


# WRAPPER ITSELF 

def main() -> None:
    throughputs_hase = []
    throughputs_original = []
    for i in range(args.n):
        with open(f"{args.name}_{i}.out") as file:
            benchmarks, throughput = parse(file)
            throughputs_original.append(throughput)
        with open(f"{args.name}_hase_{i}.out") as file:
            benchmarks, throughput = parse(file)
            throughputs_hase.append(throughput)

    throughputs_hase = np.array(throughputs_hase)
    throughputs_original = np.array(throughputs_original)

    ratios = aggregate(throughputs_hase) / aggregate(throughputs_original)

    for i in range(len(benchmarks)):
        print(f"{benchmarks[i]}\t{ratios[i]:.4f}")

    print("GeoMean\t" + str(gmean(ratios))) 

def main() -> None:
    throughputs_hase = []
    throughputs_original = []
    for i in range(args.n):
        with open(f"{args.outdir}/{args.name}_{i}.out") as file:
            benchmarks, throughput = parse(file)
            throughputs_original.append(throughput)
        with open(f"{args.outdir}/{args.name}_hase_{i}.out") as file:
            benchmarks, throughput = parse(file)
            throughputs_hase.append(throughput)

    throughputs_hase = np.array(throughputs_hase)
    throughputs_original = np.array(throughputs_original)

    ratios = aggregate(throughputs_hase) / aggregate(throughputs_original)

    for i in range(len(benchmarks)):
        print(f"{benchmarks[i]}\t{ratios[i]:.4f}")

    print("GeoMean\t" + str(gmean(ratios))) 

def main() -> None:
    throughputs_hase = []
    throughputs_original = []
    for i in range(args.n):
        with open(f"{args.name}_{i}.out") as file:
            benchmarks, throughput = parse(file)
            throughputs_original.append(throughput)
        with open(f"{args.name}_hase_{i}.out") as file:
            benchmarks, throughput = parse(file)
            throughputs_hase.append(throughput)

    throughputs_hase = np.array(throughputs_hase)
    throughputs_original = np.array(throughputs_original)

    ratios = aggregate(throughputs_hase) / aggregate(throughputs_original)

    for i in range(len(benchmarks)):
        print(f"{benchmarks[i]}\t{ratios[i]:.4f}")

    # print("GeoMean\t" + str(gmean(ratios))) 

def main() -> None:
    throughputs_hase = []
    throughputs_original = []
    for i in range(args.n):
        with open(f"{args.outdir}/{args.name}_{i}.out") as file:
            benchmarks, throughput = parse(file)
            throughputs_original.append(throughput)
        with open(f"{args.outdir}/{args.name}_hase_{i}.out") as file:
            benchmarks, throughput = parse(file)
            throughputs_hase.append(throughput)

    throughputs_hase = np.array(throughputs_hase)
    throughputs_original = np.array(throughputs_original)

    ratios = aggregate(throughputs_hase) / aggregate(throughputs_original)

    for i in range(len(benchmarks)):
        print(f"{benchmarks[i]}\t{ratios[i]:.4f}")

    print("GeoMean\t" + str(gmean(ratios))) 

def test_1D(self):
        a = (1,2,3,4)
        actual = mstats.gmean(a)
        desired = np.power(1*2*3*4,1./4.)
        assert_almost_equal(actual, desired, decimal=14)

        desired1 = mstats.gmean(a,axis=-1)
        assert_almost_equal(actual, desired1, decimal=14)
        assert_(not isinstance(desired1, ma.MaskedArray))

        a = ma.array((1,2,3,4),mask=(0,0,0,1))
        actual = mstats.gmean(a)
        desired = np.power(1*2*3,1./3.)
        assert_almost_equal(actual, desired,decimal=14)

        desired1 = mstats.gmean(a,axis=-1)
        assert_almost_equal(actual, desired1, decimal=14) 

def _aggregate_test_prediction(out_of_fold_test_predictions):
    agg_methods = {'mean': np.mean,
                   'gmean': gmean}
    prediction_column = [col for col in out_of_fold_test_predictions.columns if '_prediction' in col][0]
    if params.aggregation_method == 'rank_mean':
        rank_column = prediction_column.replace('_prediction', '_rank')
        test_predictions_with_ranks = []
        for fold_id, fold_df in out_of_fold_test_predictions.groupby('fold_id'):
            fold_df[rank_column] = calculate_rank(fold_df[prediction_column])
            test_predictions_with_ranks.append(fold_df)
        test_predictions_with_ranks = pd.concat(test_predictions_with_ranks, axis=0)

        test_prediction_aggregated = test_predictions_with_ranks.groupby(cfg.ID_COLUMNS)[rank_column].apply(
            np.mean).reset_index()
    else:
        test_prediction_aggregated = out_of_fold_test_predictions.groupby(cfg.ID_COLUMNS)[prediction_column].apply(
            agg_methods[params.aggregation_method]).reset_index()

    test_prediction_aggregated.columns = [cfg.ID_COLUMNS + cfg.TARGET_COLUMNS]
    return test_prediction_aggregated 

def geom_std(values: t.List[float]) -> float:
    """
    Calculates the geometric standard deviation for the passed values.
    Source: https://en.wikipedia.org/wiki/Geometric_standard_deviation
    """
    import scipy.stats as stats
    import scipy as sp
    gmean = stats.gmean(values)
    return sp.exp(sp.sqrt(sp.sum([sp.log(x / gmean) ** 2 for x in values]) / len(values))) 

def gmean(X):
    """ geometric mean for each variable """
    return stats.gmean(X, axis=1) 

def agg_method(self):
        methods = {'mean': np.mean,
                   'max': np.max,
                   'min': np.min,
                   'gmean': gmean
                   }
        return partial(methods[self.method], axis=-1) 

def returns_gmean(returns):
    """
    Calculates geometric average returns from a given returns series.
    """
    if isinstance(returns, pd.DataFrame) or isinstance(returns, pd.Series):
        returns = returns.fillna(0)
    else:
        returns = np.nan_to_num(returns)
    return ss.gmean(1 + returns, axis=0) - 1 

def log_to_pct_return(log_returns):
    return np.exp(log_returns) - 1


# def validate_mean_method(method):
#     if method not in {'gmean', 'simple'}:
#         raise AssertionError('mean_method can only be {"gmean", "simple"}')


# def validate_return_type(return_type):
#     if return_type not in {'pct', 'log'}:
#         raise AssertionError('mean_method can only be {"pct", "log"}') 

def agg_method(self):
        methods = {'mean': np.mean,
                   'max': np.max,
                   'min': np.min,
                   'gmean': gmean
                   }
        return partial(methods[self.method], axis=-1) 

def easystep(self, action):
        self.counter += 1

        self.upd_env_R_from_R(action)

        # write to file input for Omnet: Routing
        vector_to_file(matrix_to_omnet_v(self.env_R), self.folder + OMROUTING, 'w')
        # VERIFY FILE POSITION AND FORMAT (separator, matrix/vector) np.savetxt("tmp.txt", routing, fmt="%d")

        # execute omnet
        omnet_wrapper(self)

        # read Omnet's output: Delay and Lost
        om_output = file_to_csv(self.folder + OMDELAY)
        self.upd_env_D(csv_to_matrix(om_output, self.ACTIVE_NODES))
        self.upd_env_L(csv_to_lost(om_output))

        reward = rl_reward(self)

        # log everything to file
        vector_to_file([-reward], self.folder + REWARDLOG, 'a')
        cur_state = rl_state(self)
        log = np.concatenate(([self.counter], [self.env_L], [np.mean(matrix_to_rl(self.env_D))], [np.max(matrix_to_rl(self.env_D))], [(np.mean(matrix_to_rl(self.env_D)) + np.max(matrix_to_rl(self.env_D)))/2], [stats.gmean(matrix_to_rl(self.env_D))]))
        vector_to_file(log, self.folder + WHOLELOG, 'a')

        # generate traffic for next iteration
        self.upd_env_T(self.tgen.generate())
        # write to file input for Omnet: Traffic, or do nothing if static
        if self.TRAFFIC.split(':')[0] not in ('STAT', 'STATEQ', 'FILE', 'DIR'):
            vector_to_file(matrix_to_omnet_v(self.env_T), self.folder + OMTRAFFIC, 'w')

        new_state = rl_state(self)
        # return new status and reward
        return new_state, reward, 0 

def used_summarize_mean(values: t.List[float]) -> float:
    if CALC_MODE in [Mode.geom_mean_rel_to_best, Mode.mean_rel_to_one]:
        return stats.gmean(values)
    elif CALC_MODE in [Mode.mean_rel_to_first]:
        return sp.mean(values)
    assert False 

def get_geom_over_rel_means(self) -> t.Dict[str, float]:
        return self.get_reduced_x_per_impl(used_rel_mean_property, stats.gmean) 

def get_geom_over_rel_stds(self) -> t.Dict[str, float]:
        return self.get_reduced_x_per_impl(rel_std_property, stats.gmean) 

def get_html(self, base_file_name: str, h_level: int) -> str:
        html = """
            <h{}>Program: {!r} ({} lines, {} entropy)</h{}>
            The following plot shows the mean score per input distribution for every implementation.
        """.format(h_level, self.name, self.line_number, self.entropy, h_level)
        html += self.get_box_plot_html(base_file_name)
        scores = self.get_impl_mean_scores()
        std_devs = self.get_statistical_property_scores(rel_std_dev_func)
        html += """
            <table class="table">
                <tr><th>implementation</th><th>geom mean over means relative to best (per input) aka mean score</th>
                    <th>... std dev rel. to the best mean</th>
                </tr>
        """
        for impl in scores.keys():
            html += """
                <tr><td>{}</td><td>{:5.2%}</td><td>{:5.2%}</td></tr>
            """.format(impl, stats.gmean(scores[impl]), stats.gmean(std_devs[impl]))
        html += "</table>"
        impl_names = list(scores.keys())
        for (i, input) in enumerate(self.prog_inputs.keys()):
            app = html_escape_property(input)
            if len(app) > 20:
                app = str(i)
            html += self.prog_inputs[input].get_html(base_file_name + "_" + app, h_level + 1)
        return html 

def get_statistical_property_scores_per_impl(self, func: StatisticalPropertyFunc,
                                                 reduce: ReduceFunc = stats.gmean) -> t.Dict[str, float]:
        impl_scores = InsertionTimeOrderedDict()
        for prog in self.programs:
            scores = self.programs[prog].get_statistical_property_scores(func)
            for impl in scores:
                if impl not in impl_scores:
                    impl_scores[impl] = []
                impl_scores[impl].append(reduce(scores[impl]))
        return impl_scores 

def get_statistical_property_scores(self, func: StatisticalPropertyFunc,
                                        reduce: ReduceFunc = stats.gmean) -> t.Dict[str, float]:
        ret = InsertionTimeOrderedDict()
        scores_per_impl = self.get_statistical_property_scores_per_impl(func)
        for impl in scores_per_impl:
            ret[impl] = reduce(scores_per_impl[impl])
        return ret 

def get_statistical_property_scores(self, func: StatisticalPropertyFunc,
                                        reduce: ReduceFunc = stats.gmean) -> t.Dict[str, float]:
        ret = InsertionTimeOrderedDict()
        scores_per_impl = self.get_statistical_property_scores_per_impl(func)
        for impl in scores_per_impl:
            ret[impl] = reduce(scores_per_impl[impl])
        return ret 

def _report_block_with_baseline(self, block: RunData, print_func: t.Callable[[str], None], baseline: RunData,
                                    descr_size: int):
        print_func("{descr:<20} ({num:>5}) with baseline {descr2:<20} ({num2:>5})"
                   .format(descr=block.description(), num=len(block.data[block.properties[0]]),
                           descr2=block.description(), num2=len(block.data[block.properties[0]])))
        combined_props = set(block.properties) & set(baseline.properties)
        tester = TesterRegistry.get_tester()
        for prop in sorted(combined_props):
            mean = np.mean(block[prop])
            std = np.std(block[prop])
            base_mean = baseline.get_single_properties()[prop].mean()
            base_std = baseline.get_single_properties()[prop].std()
            mean_str = str(FNumber(mean / base_mean, abs_deviation=std / base_mean, is_percent=True))
            dev = "{:>5.5%}".format(std / mean) if mean != 0 else "{:>5.5}".format(std)
            print_func("\t {{prop:<{}}} mean = {{mean:>15s}}, confidence = {{conf:>5.0%}}, dev = {{dev:>11s}}, "
                       "{{dbase:>11s}}".format(descr_size)
                .format(
                    prop=prop,
                    mean=mean_str,
                    dev=dev,
                    conf=tester.test(block[prop], baseline[prop]),
                    dbase="{:>5.5%}".format(base_std / base_mean) if base_mean != 0 else "{:>5.5}".format(base_std)))
        rels = [(block.get_single_properties()[prop].mean() / baseline.get_single_properties()[prop].mean())
                            for prop in combined_props]
        gmean = stats.gmean(rels)
        gstd = util.geom_std(rels)
        print_func("geometric mean of relative mean = {:>15}, std dev = {:>15}"
                   .format(FNumber(gmean, is_percent=True).format(), FNumber(gstd, is_percent=True).format())) 

def test_2D(self):
        a = ma.array(((1, 2, 3, 4), (1, 2, 3, 4), (1, 2, 3, 4)),
                     mask=((0, 0, 0, 0), (1, 0, 0, 1), (0, 1, 1, 0)))
        actual = mstats.gmean(a)
        desired = np.array((1,2,3,4))
        assert_array_almost_equal(actual, desired, decimal=14)

        desired1 = mstats.gmean(a,axis=0)
        assert_array_almost_equal(actual, desired1, decimal=14)

        actual = mstats.gmean(a, -1)
        desired = ma.array((np.power(1*2*3*4,1./4.),
                            np.power(2*3,1./2.),
                            np.power(1*4,1./2.)))
        assert_array_almost_equal(actual, desired, decimal=14) 

def agg_method(self):
        methods = {'mean': np.mean,
                   'max': np.max,
                   'min': np.min,
                   'gmean': gmean
                   }
        return partial(methods[self.method], axis=-1) 

def agg_method(self):
        methods = {'mean': np.mean,
                   'max': np.max,
                   'min': np.min,
                   'gmean': gmean
                   }
        return partial(methods[self.method], axis=-1) 

def gmean(self):
        """Returns the gmean of the models predictions.

        Returns
        -------
        `PipeApply`
        """
        return self.apply(lambda x: gmean(x, axis=0)) 

def test_1D_list(self):
        a = (1,2,3,4)
        actual = stats.gmean(a)
        desired = power(1*2*3*4,1./4.)
        assert_almost_equal(actual, desired,decimal=14)

        desired1 = stats.gmean(a,axis=-1)
        assert_almost_equal(actual, desired1, decimal=14) 

def test_1D_array(self):
        a = array((1,2,3,4), float32)
        actual = stats.gmean(a)
        desired = power(1*2*3*4,1./4.)
        assert_almost_equal(actual, desired, decimal=7)

        desired1 = stats.gmean(a,axis=-1)
        assert_almost_equal(actual, desired1, decimal=7)