text:
"""
Code based on Alex Lang and Oscar Beijbom, 2018.
Licensed under MIT License [see LICENSE].
"""

import time
from enum import Enum
from functools import reduce

import numpy as np
import sparseconvnet as scn
import torch
from torch import nn
from torch.nn import functional as F
import libs
from libs.tools import metrics

from models.bones.pillars import (PillarFeatureNet, PointPillarsScatter,
                                  PillarFeatureNet_trt, PointPillarsScatter_trt)
from models.bones.rpn import RPN,RPN_trt

from core.losses import (WeightedSigmoidClassificationLoss,
                        WeightedSmoothL1LocalizationLoss,
                        WeightedSoftmaxClassificationLoss)
from libs.ops import box_np_ops,box_torch_ops

class PointPillars(nn.Module):
    def __init__(self,
                output_shape,
                model_cfg,
                target_assigner):
        super().__init__()

        self.name =model_cfg.NAME
        self._pc_range = model_cfg.POINT_CLOUD_RANGE
        self._voxel_size = model_cfg.GRID_SIZE
        self._num_class = model_cfg.NUM_CLASS
        self._use_bev = model_cfg.BACKBONE.use_bev
        self._total_forward_time = 0.0
        self._total_postprocess_time = 0.0
        self._total_inference_count = 0

        #for prepare loss weights
        self._pos_cls_weight = model_cfg.pos_cls_weight
        self._neg_cls_weight = model_cfg.neg_cls_weight
        self._loss_norm_type = model_cfg.loss_norm_type
        #for create loss
        self._loc_loss_ftor = model_cfg.loc_loss_ftor
        self._cls_loss_ftor = model_cfg.cls_loss_ftor
        self._dir_loss_ftor = WeightedSoftmaxClassificationLoss()
        self._direction_loss_weight = model_cfg.LOSS.direction_loss_weight
        self._encode_rad_error_by_sin = model_cfg.ENCODE_RAD_ERROR_BY_SIN
        #
        self._cls_loss_weight = model_cfg.cls_weight
        self._loc_loss_weight = model_cfg.loc_weight
        # for direction classifier
        self._use_direction_classifier = model_cfg.BACKBONE.use_direction_classifier
        # for predict
        self._use_sigmoid_score = model_cfg.POST_PROCESSING.use_sigmoid_score
        self._box_coder = target_assigner.box_coder
        self.target_assigner = target_assigner
        #for nms
        self._multiclass_nms = model_cfg.PREDICT.multiclass_nms
        self._use_rotate_nms = model_cfg.PREDICT.use_rotate_nms

        self._nms_score_threshold = model_cfg.POST_PROCESSING.nms_score_threshold
        self._nms_pre_max_size = model_cfg.POST_PROCESSING.nms_pre_max_size
        self._nms_post_max_size = model_cfg.POST_PROCESSING.nms_post_max_size
        self._nms_iou_threshold = model_cfg.POST_PROCESSING.nms_iou_threshold
        self._use_sigmoid_score = model_cfg.POST_PROCESSING.use_sigmoid_score
        # self._use_sigmoid_score = use_sigmoid_score
        self._encode_background_as_zeros=model_cfg.BACKBONE.encode_background_as_zeros
        #1.PFN
        self.pfn = PillarFeatureNet(model_cfg.num_input_features,
                                    model_cfg.PILLAR_FEATURE_EXTRACTOR.use_norm,
                                    num_filters=model_cfg.pfn_num_filters,
                                    with_distance=model_cfg.PILLAR_FEATURE_EXTRACTOR.with_distance,
                                    voxel_size = self._voxel_size,
                                    pc_range = self._pc_range)

        #2.sparse middle
        self.mfe = PointPillarsScatter(output_shape=output_shape,
                                       num_input_features=model_cfg.pfn_num_filters[-1])
        num_rpn_input_filters=self.mfe.nchannels
        #3.rpn
        self.rpn =  RPN(
                        use_norm=model_cfg.BACKBONE.use_norm,
                        num_class=self._num_class,
                        layer_nums=model_cfg.BACKBONE.layer_nums,
                        layer_strides=model_cfg.BACKBONE.layer_strides,
                        num_filters=model_cfg.BACKBONE.num_filters,
                        upsample_strides=model_cfg.BACKBONE.upsample_strides,
                        num_upsample_filters=model_cfg.BACKBONE.num_upsample_filters,
                        num_input_filters=num_rpn_input_filters,
                        num_anchor_per_loc=target_assigner.num_anchors_per_location,
                        encode_background_as_zeros=self._encode_background_as_zeros,
                        use_direction_classifier=model_cfg.BACKBONE.use_direction_classifier,
                        use_bev=model_cfg.BACKBONE.use_bev,
                        use_groupnorm=model_cfg.BACKBONE.use_groupnorm,
                        num_groups=model_cfg.BACKBONE.num_groups,
                        box_code_size=target_assigner.box_coder.code_size
                        )
        self.rpn_acc = metrics.Accuracy(
            dim=-1, encode_background_as_zeros=self._encode_background_as_zeros)
        self.rpn_precision = metrics.Precision(dim=-1)
        self.rpn_recall = metrics.Recall(dim=-1)
        self.rpn_metrics = metrics.PrecisionRecall(
            dim=-1,
            thresholds=[0.1, 0.3, 0.5, 0.7, 0.8, 0.9, 0.95],
            use_sigmoid_score=self._use_sigmoid_score,
            encode_background_as_zeros=self._encode_background_as_zeros)

        self.rpn_cls_loss =metrics.Scalar()
        self.rpn_loc_loss =metrics.Scalar()
        self.rpn_total_loss = metrics.Scalar()
        self.register_buffer("global_step", torch.LongTensor(1).zero_())

    def update_global_step(self):
        self.global_step += 1

    def get_global_step(self):
        return int(self.global_step.cpu().numpy()[0])

    def forward(self,example):
        '''nutonomy/second.pytorch'''
        voxels = example['voxels']
        num_points = example['num_points']
        coors = example["coordinates"]
        batch_anchors = example["anchors"]
        batch_size_dev = batch_anchors.shape[0]
        t = time.time()
        # features: [num_voxels, max_num_points_per_voxel, 7]
        # num_points: [num_voxels]
        # coors: [num_voxels, 4]
        voxel_features = self.pfn(voxels, num_points, coors)

        spatial_features = self.mfe(voxel_features,coors,batch_size_dev)
        if self._use_bev:
            preds_dict = self.rpn(spatial_features,example['bev_map'])
        else:
            preds_dict = self.rpn(spatial_features)

        box_preds = preds_dict["box_preds"]
        cls_preds = preds_dict["cls_preds"]
        self._total_forward_time += time.time() - t
        if self.training:
            labels = example['labels']
            reg_targets = example['reg_targets']

            cls_weights, reg_weights, cared = prepare_loss_weights(
                labels,
                pos_cls_weight=self._pos_cls_weight,
                neg_cls_weight=self._neg_cls_weight,
                loss_norm_type=self._loss_norm_type,
                dtype=voxels.dtype)
            cls_targets = labels * cared.type_as(labels)
            cls_targets = cls_targets.unsqueeze(-1)

            loc_loss,cls_loss = create_loss(
                self._loc_loss_ftor,
                self._cls_loss_ftor,
                box_preds=box_preds,
                cls_preds=cls_preds,
                cls_targets=cls_targets,
                cls_weights=cls_weights,
                reg_targets=reg_targets,
                reg_weights=reg_weights,
                num_class=self._num_class,
                encode_rad_error_by_sin=self._encode_rad_error_by_sin,
                encode_background_as_zeros=self._encode_background_as_zeros,
                box_code_size=self._box_coder.code_size,)

            loc_loss_reduced = loc_loss.sum() / batch_size_dev
            loc_loss_reduced *= self._loc_loss_weight
            cls_pos_loss, cls_neg_loss = _get_pos_neg_loss(cls_loss, labels)
            cls_pos_loss /= self._pos_cls_weight
            cls_neg_loss /= self._neg_cls_weight
            cls_loss_reduced = cls_loss.sum() / batch_size_dev
            cls_loss_reduced *= self._cls_loss_weight
            loss = loc_loss_reduced + cls_loss_reduced
            if self._use_direction_classifier:
                dir_targets =  get_direction_target(example['anchors'],
                                                   reg_targets)
                dir_logits = preds_dict["dir_cls_preds"].view(
                    batch_size_dev, -1, 2)
                weights = (labels > 0).type_as(dir_logits)
                weights /= torch.clamp(weights.sum(-1, keepdim=True), min=1.0)
                dir_loss = self._dir_loss_ftor(
                    dir_logits, dir_targets, weights=weights)
                dir_loss = dir_loss.sum() / batch_size_dev
                loss += dir_loss * self._direction_loss_weight

            return {
                "loss": loss,
                "cls_loss": cls_loss,
                "loc_loss": loc_loss,
                "cls_pos_loss": cls_pos_loss,
                "cls_neg_loss": cls_neg_loss,
                "cls_preds": cls_preds,
                "dir_loss_reduced": dir_loss,
                "cls_loss_reduced": cls_loss_reduced,
                "loc_loss_reduced": loc_loss_reduced,
                "cared": cared,
            }
        else:
            return self.predict(example, preds_dict)

    def predict(self, example, preds_dict):
        t = time.time()
        batch_size = example['anchors'].shape[0]
        batch_anchors = example["anchors"].view(batch_size, -1, 7)

        self._total_inference_count += batch_size
        if "anchors_mask" not in example:
            batch_anchors_mask = [None] * batch_size
        else:
            batch_anchors_mask = example["anchors_mask"].view(batch_size, -1)
        batch_imgidx = example['image_idx']

        self._total_forward_time += time.time() - t
        t = time.time()

        batch_box_preds = preds_dict["box_preds"]
        batch_cls_preds = preds_dict["cls_preds"]
        batch_box_preds = batch_box_preds.view(batch_size, -1,
                                               self._box_coder.code_size)

        num_class_with_bg = self._num_class
        if not self._encode_background_as_zeros:
            num_class_with_bg = self._num_class + 1

        batch_cls_preds = batch_cls_preds.view(batch_size, -1,
                                               num_class_with_bg)
        batch_box_preds = self._box_coder.decode_torch(batch_box_preds,
                                                       batch_anchors)
        if self._use_direction_classifier:
            batch_dir_preds = preds_dict["dir_cls_preds"]
            batch_dir_preds = batch_dir_preds.view(batch_size, -1, 2)
        else:
            batch_dir_preds = [None] * batch_size

        predictions_dicts = []

        for box_preds, cls_preds, dir_preds, img_idx, a_mask in zip(
            batch_box_preds, batch_cls_preds, batch_dir_preds,batch_imgidx, batch_anchors_mask):

            if a_mask is not None:
                box_preds = box_preds[a_mask]
                cls_preds = cls_preds[a_mask]
            if self._use_direction_classifier:
                if a_mask is not None:
                    dir_preds = dir_preds[a_mask]
                dir_labels = torch.max(dir_preds,dim = -1)[1]
            if self._encode_background_as_zeros:
                # this don't support softmax
                assert self._use_sigmoid_score is True
                total_scores = torch.sigmoid(cls_preds)
            else:
                # encode background as first element in one-hot vector
                if self._use_sigmoid_score:
                    total_scores = torch.sigmoid(cls_preds)[..., 1:]
                else:
                    total_scores = F.softmax(cls_preds, dim=-1)[..., 1:]

            #apply nms in birdeye view
            if self._use_rotate_nms:
                nms_func = box_torch_ops.rotate_nms
            else:
                nms_func = box_torch_ops.nms
            selected_boxes = None
            selected_labels = None
            selected_scores = None
            selected_dir_labels = None

            if self._multiclass_nms:
                #curently only support class-agnostic boxes.
                boxes_for_nms = box_preds[:, [0, 1, 3, 4, 6]]
                if not self._use_rotate_nms:
                    box_preds_corners = box_torch_ops.center_to_corner_box2d(
                        boxes_for_nms[:,:2],boxes_for_nms[:,2:4],
                        boxes_for_nms[:,4]
                    )
                    box_for_nms = box_torch_ops.corner_to_standup_nd(box_preds_corners)
                boxes_for_mcnms = boxes_for_nms.unsqueeze(1)
                selected_per_class = box_torch_ops.multiclass_nms(
                    nms_func=nms_func,
                    boxes=boxes_for_mcnms,
                    scores=total_scores,
                    num_class=self._num_class,
                    pre_max_size=self._nms_pre_max_size,
                    post_max_size=self._nms_post_max_size,
                    iou_threshold=self._nms_iou_threshold,
                    score_thresh=self._nms_score_threshold,
                )
                selected_boxes, selected_labels, selected_scores = [] , [], []
                selected_dir_labels = []
                for i , selected in enumerate(selected_per_class):
                    if selected is not None:
                        num_dets = selected.shape[0]
                        selected_boxes.append(box_preds[selected])
                        selected_labels.append(
                            torch.full([num_dets],i,dtype=torch.int64)
                        )
                        if len(selected_boxes) > 0:
                            selected_boxes = torch.cat(selected_boxes, dim= 0)
                            selected_labels = torch.cat(selected_labels,dim=0)
                            selected_scores = torch.cat(selected_scores,dim = 0)
                            if self._use_direction_classifier:
                                selected_dir_labels = torch.cat(
                                    selected_dir_labels,dim= 0
                                )
                else:
                    selected_boxes = None
                    selected_labels = None
                    selected_scores = None
                    selected_dir_labels = None
            else:
                # get highest score per prediction, than apply nms
                # to remove overlapped box
                if num_class_with_bg ==1:
                    top_scores = total_scores.squeeze(-1)
                    top_labels = torch.zeros(total_scores.shape[0],device=total_scores.device,dtype=torch.long)
                else:
                    top_scores, top_labels = torch.max(total_scores,dim=-1)

                if self._nms_score_threshold > 0.0:
                    thresh = torch.tensor(
                        [self._nms_score_threshold],
                        device = total_scores.device
                    ).type_as(total_scores)
                    top_scores_keep = (top_scores >= thresh)
                    top_scores = top_scores.masked_select(top_scores_keep)
                if top_scores.shape[0] != 0:
                    if self._nms_score_threshold > 0.0:
                        box_preds = box_preds[top_scores_keep]
                        if self._use_direction_classifier:
                            dir_labels = dir_labels[top_scores_keep]
                        top_labels = top_labels[top_scores_keep]
                    boxes_for_nms = box_preds[:,[0,1,3,4,6]]
                    if not self._use_rotate_nms:
                        box_preds_corners = box_torch_ops.center_to_corner_box2d(
                            boxes_for_nms[:, :2], boxes_for_nms[:, 2:4],
                            boxes_for_nms[:, 4])
                        boxes_for_nms = box_torch_ops.corner_to_standup_nd(
                            box_preds_corners)
                    # the nms in 3d detection just remove overlap boxes.
                    selected = nms_func(
                        boxes_for_nms,
                        top_scores,
                        pre_max_size=self._nms_pre_max_size,
                        post_max_size=self._nms_post_max_size,
                        iou_threshold=self._nms_iou_threshold,
                    )
                else:
                    selected = None
                if selected is not None:
                    selected_boxes = box_preds[selected]
                    if self._use_direction_classifier:
                        selected_dir_labels = dir_labels[selected]
                    selected_labels = top_labels[selected]
                    selected_scores = top_scores[selected]

            #finally generate predictions.
            if selected_boxes is not None:
                box_preds = selected_boxes
                scores = selected_scores
                label_preds = selected_labels
                if self._use_direction_classifier:
                    dir_labels = selected_dir_labels
                    opp_labels = (box_preds[..., -1] > 0) ^ dir_labels.bool()
                    box_preds[...,-1] += torch.where(
                        opp_labels,
                        torch.tensor(np.pi).type_as(box_preds),
                        torch.tensor(0.0).type_as(box_preds)
                    )
                final_box_preds = box_preds
                final_scores = scores
                final_labels = label_preds
                # predictions
                predictions_dict = {
                    "bbox": None,
                    "box3d_camera": None ,
                    "box3d_lidar": final_box_preds,
                    "scores": final_scores,
                    "label_preds": final_labels,
                    "image_idx": img_idx,
                }
            else:
                predictions_dict = {
                    "bbox": None,
                    "box3d_camera": None,
                    "box3d_lidar": None,
                    "scores": None,
                    "label_preds": None,
                    "image_idx": img_idx,
                }
            predictions_dicts.append(predictions_dict)
        self._total_postprocess_time += time.time() - t
        return predictions_dicts

    @property
    def avg_forward_time(self):
        return self._total_forward_time / self._total_inference_count

    @property
    def avg_postprocess_time(self):
        return self._total_postprocess_time / self._total_inference_count

    def clear_time_metrics(self):
        self._total_forward_time = 0.0
        self._total_postprocess_time = 0.0
        self._total_inference_count = 0

    def metrics_to_float(self):
        self.rpn_acc.float()
        self.rpn_metrics.float()
        self.rpn_cls_loss.float()
        self.rpn_loc_loss.float()
        self.rpn_total_loss.float()

    def update_metrics(self,
                       cls_loss,
                       loc_loss,
                       cls_preds,
                       labels,
                       sampled):
        batch_size = cls_preds.shape[0]
        num_class = self._num_class
        if not self._encode_background_as_zeros:
            num_class += 1
        cls_preds = cls_preds.view(batch_size, -1, num_class)
        rpn_acc = self.rpn_acc(labels, cls_preds, sampled).numpy()[0]
        prec, recall = self.rpn_metrics(labels, cls_preds, sampled)
        prec = prec.numpy()
        recall = recall.numpy()
        rpn_cls_loss = self.rpn_cls_loss(cls_loss).numpy()[0]
        rpn_loc_loss = self.rpn_loc_loss(loc_loss).numpy()[0]
        ret = {
            "cls_loss": float(rpn_cls_loss),
            "cls_loss_rt": float(cls_loss.data.cpu().numpy()),
            'loc_loss': float(rpn_loc_loss),
            "loc_loss_rt": float(loc_loss.data.cpu().numpy()),
            "rpn_acc": float(rpn_acc),
        }
        for i, thresh in enumerate(self.rpn_metrics.thresholds):
            ret[f"prec@{int(thresh*100)}"] = float(prec[i])
            ret[f"rec@{int(thresh*100)}"] = float(recall[i])
        return ret

    def clear_metrics(self):
        self.rpn_acc.clear()
        self.rpn_metrics.clear()
        self.rpn_cls_loss.clear()
        self.rpn_loc_loss.clear()
        self.rpn_total_loss.clear()

    @staticmethod
    def convert_norm_to_float(net):
        '''
        BatchNorm layers to have parameters in single precision.
        Find all layers and convert them back to float. This can't
        be done with built in .apply as that function will apply
        fn to all modules, parameters, and buffers. Thus we wouldn't
        be able to guard the float conversion based on the module type.
        '''
        if isinstance(net, torch.nn.modules.batchnorm._BatchNorm):
            net.float()
        for child in net.children():
            PointPillars.convert_norm_to_float(net)
        return net

class PointPillars_trt(nn.Module):
    def __init__(self,
                output_shape,
                model_cfg,
                target_assigner):
        super().__init__()

        self.name =model_cfg.NAME
        self._pc_range = model_cfg.POINT_CLOUD_RANGE
        self._voxel_size = model_cfg.GRID_SIZE
        self._num_class = model_cfg.NUM_CLASS
        self._use_bev = model_cfg.BACKBONE.use_bev
        self._total_forward_time = 0.0
        self._total_postprocess_time = 0.0
        self._total_inference_count = 0

        # #for prepare loss weights
        # self._pos_cls_weight = model_cfg.pos_cls_weight
        # self._neg_cls_weight = model_cfg.neg_cls_weight
        # self._loss_norm_type = model_cfg.loss_norm_type
        # #for create loss
        # self._loc_loss_ftor = model_cfg.loc_loss_ftor
        # self._cls_loss_ftor = model_cfg.cls_loss_ftor
        # self._dir_loss_ftor = WeightedSoftmaxClassificationLoss()
        # self._direction_loss_weight = model_cfg.LOSS.direction_loss_weight
        # self._encode_rad_error_by_sin = model_cfg.ENCODE_RAD_ERROR_BY_SIN
        # #
        # self._cls_loss_weight = model_cfg.cls_weight
        # self._loc_loss_weight = model_cfg.loc_weight
        # for direction classifier
        self._use_direction_classifier = model_cfg.BACKBONE.use_direction_classifier
        # for predict
        self._use_sigmoid_score = model_cfg.POST_PROCESSING.use_sigmoid_score
        self._box_coder = target_assigner.box_coder
        self.target_assigner = target_assigner
        #for nms
        self._multiclass_nms = model_cfg.PREDICT.multiclass_nms
        self._use_rotate_nms = model_cfg.PREDICT.use_rotate_nms

        self._nms_score_threshold = model_cfg.POST_PROCESSING.nms_score_threshold
        self._nms_pre_max_size = model_cfg.POST_PROCESSING.nms_pre_max_size
        self._nms_post_max_size = model_cfg.POST_PROCESSING.nms_post_max_size
        self._nms_iou_threshold = model_cfg.POST_PROCESSING.nms_iou_threshold
        self._use_sigmoid_score = model_cfg.POST_PROCESSING.use_sigmoid_score
        # self._use_sigmoid_score = use_sigmoid_score
        self._encode_background_as_zeros=model_cfg.BACKBONE.encode_background_as_zeros
        #1.PFN
        self.pfn = PillarFeatureNet_trt(model_cfg.num_input_features,
                                    model_cfg.PILLAR_FEATURE_EXTRACTOR.use_norm,
                                    num_filters=model_cfg.pfn_num_filters,
                                    with_distance=model_cfg.PILLAR_FEATURE_EXTRACTOR.with_distance,
                                    voxel_size = self._voxel_size,
                                    pc_range = self._pc_range)

        #2.sparse middle
        self.mfe = PointPillarsScatter_trt(output_shape=output_shape,
                                       num_input_features=model_cfg.pfn_num_filters[-1])
        num_rpn_input_filters=self.mfe.nchannels
        #3.rpn
        self.rpn =  RPN_trt(
                        use_norm=model_cfg.BACKBONE.use_norm,
                        num_class=self._num_class,
                        layer_nums=model_cfg.BACKBONE.layer_nums,
                        layer_strides=model_cfg.BACKBONE.layer_strides,
                        num_filters=model_cfg.BACKBONE.num_filters,
                        upsample_strides=model_cfg.BACKBONE.upsample_strides,
                        num_upsample_filters=model_cfg.BACKBONE.num_upsample_filters,
                        num_input_filters=num_rpn_input_filters,
                        num_anchor_per_loc=target_assigner.num_anchors_per_location,
                        encode_background_as_zeros=self._encode_background_as_zeros,
                        use_direction_classifier=model_cfg.BACKBONE.use_direction_classifier,
                        use_bev=model_cfg.BACKBONE.use_bev,
                        use_groupnorm=model_cfg.BACKBONE.use_groupnorm,
                        num_groups=model_cfg.BACKBONE.num_groups,
                        box_code_size=target_assigner.box_coder.code_size
                        )

    def forward(self,example):
        '''nutonomy/second.pytorch'''
        # training input [0:pillar_x, 1:pillar_y, 2:pillar_z, 3:pillar_i,
        #                 4:num_points_per_pillar, 5:x_sub_shaped, 6:y_sub_shaped, 7:mask, 8:coors
        #                 9:voxel_mask, 10:anchors, #11:anchors_mask]
        #print(len(example))
        pillar_x = example[0]
        pillar_y = example[1]
        pillar_z = example[2]
        pillar_i = example[3]
        num_points = example[4]
        x_sub_shaped = example[5]
        y_sub_shaped = example[6]
        mask = example[7]
        #print(pillar_x.shape)

        # features: [num_voxels, max_num_points_per_voxel, 7]
        # num_points: [num_voxels]
        # coors: [num_voxels, 4]
        voxel_features = self.pfn(pillar_x, pillar_y, pillar_z, pillar_i,
                                  num_points, x_sub_shaped, y_sub_shaped, mask)
        voxel_features = voxel_features.squeeze()
        voxel_features = voxel_features.permute(1, 0)
        #print(voxel_features.shape)
        coors = example[8]
        spatial_features = self.mfe(voxel_features, coors,example[9])
        # spatial_features input size is : [1, 64, 496, 432]

        preds_dict = self.rpn(spatial_features)

        return self.predict(example, preds_dict)

    def predict(self, example, preds_dict):
        batch_size = 1
        batch_anchors = example[10].view(batch_size, -1, 7)

        self._total_inference_count += batch_size
        # batch_rect = example[11]
        # batch_Trv2c = example[12]
        # batch_P2 = example[13]

        batch_anchors_mask = [None] * batch_size

        # if "anchors_mask" not in example:
        #     batch_anchors_mask = [None] * batch_size
        # else:
        #     batch_anchors_mask = example["anchors_mask"].view(batch_size, -1)
        # assert 15==len(example), "somthing write with example size!"
        #batch_anchors_mask = anchors_mask.view(batch_size, -1)
        # batch_imgidx = example['image_idx']
        batch_imgidx = example[11]

        # self._total_forward_time += time.time() - t
        # t = time.time()
        batch_box_preds = preds_dict[0]
        batch_cls_preds = preds_dict[1]
        batch_box_preds = batch_box_preds.view(batch_size, -1,
                                               self._box_coder.code_size)
        num_class_with_bg = self._num_class
        if not self._encode_background_as_zeros:
            num_class_with_bg = self._num_class + 1

        batch_cls_preds = batch_cls_preds.view(batch_size, -1,
                                               num_class_with_bg)
        batch_box_preds = self._box_coder.decode_torch(batch_box_preds,
                                                       batch_anchors)
        if self._use_direction_classifier:
            batch_dir_preds = preds_dict[2]
            batch_dir_preds = batch_dir_preds.view(batch_size, -1, 2)
        else:
            batch_dir_preds = [None] * batch_size

        # predictions_dicts = []
        predictions_dicts = ()
        for box_preds, cls_preds, dir_preds,  img_idx, a_mask in zip(
                batch_box_preds, batch_cls_preds, batch_dir_preds, batch_imgidx, batch_anchors_mask):
            if a_mask is not None:
                box_preds = box_preds[a_mask]
                cls_preds = cls_preds[a_mask]
            if self._use_direction_classifier:
                if a_mask is not None:
                    dir_preds = dir_preds[a_mask]
                # print(dir_preds.shape)
                dir_labels = torch.max(dir_preds, dim=-1)[1]
            if self._encode_background_as_zeros:
                # this don't support softmax
                assert self._use_sigmoid_score is True
                total_scores = torch.sigmoid(cls_preds)
            else:
                # encode background as first element in one-hot vector
                if self._use_sigmoid_score:
                    total_scores = torch.sigmoid(cls_preds)[..., 1:]
                else:
                    total_scores = F.softmax(cls_preds, dim=-1)[..., 1:]
            # Apply NMS in birdeye view
            if self._use_rotate_nms:
                nms_func = box_torch_ops.rotate_nms
            else:
                nms_func = box_torch_ops.nms
            selected_boxes = None
            selected_labels = None
            selected_scores = None
            selected_dir_labels = None

            if self._multiclass_nms:
                # curently only support class-agnostic boxes.
                boxes_for_nms = box_preds[:, [0, 1, 3, 4, 6]]
                if not self._use_rotate_nms:
                    box_preds_corners = box_torch_ops.center_to_corner_box2d(
                        boxes_for_nms[:, :2], boxes_for_nms[:, 2:4],
                        boxes_for_nms[:, 4])
                    boxes_for_nms = box_torch_ops.corner_to_standup_nd(
                        box_preds_corners)
                boxes_for_mcnms = boxes_for_nms.unsqueeze(1)
                selected_per_class = box_torch_ops.multiclass_nms(
                    nms_func=nms_func,
                    boxes=boxes_for_mcnms,
                    scores=total_scores,
                    num_class=self._num_class,
                    pre_max_size=self._nms_pre_max_size,
                    post_max_size=self._nms_post_max_size,
                    iou_threshold=self._nms_iou_threshold,
                    score_thresh=self._nms_score_threshold,
                )
                selected_boxes, selected_labels, selected_scores = [], [], []
                selected_dir_labels = []
                for i, selected in enumerate(selected_per_class):
                    if selected is not None:
                        num_dets = selected.shape[0]
                        selected_boxes.append(box_preds[selected])
                        selected_labels.append(
                            torch.full([num_dets], i, dtype=torch.int64))
                        if self._use_direction_classifier:
                            selected_dir_labels.append(dir_labels[selected])
                        selected_scores.append(total_scores[selected, i])
                if len(selected_boxes) > 0:
                    selected_boxes = torch.cat(selected_boxes, dim=0)
                    selected_labels = torch.cat(selected_labels, dim=0)
                    selected_scores = torch.cat(selected_scores, dim=0)
                    if self._use_direction_classifier:
                        selected_dir_labels = torch.cat(
                            selected_dir_labels, dim=0)
                else:
                    selected_boxes = None
                    selected_labels = None
                    selected_scores = None
                    selected_dir_labels = None
            else:
                # get highest score per prediction, than apply nms
                # to remove overlapped box.
                if num_class_with_bg == 1:
                    top_scores = total_scores.squeeze(-1)
                    top_labels = torch.zeros(
                        total_scores.shape[0],
                        device=total_scores.device,
                        dtype=torch.long)
                else:
                    top_scores, top_labels = torch.max(total_scores, dim=-1)

                if self._nms_score_threshold > 0.0:
                    thresh = torch.tensor(
                        [self._nms_score_threshold],
                        device=total_scores.device).type_as(total_scores)
                    top_scores_keep = (top_scores >= thresh)
                    top_scores = top_scores.masked_select(top_scores_keep)
                if top_scores.shape[0] != 0:
                    if self._nms_score_threshold > 0.0:
                        box_preds = box_preds[top_scores_keep]
                        if self._use_direction_classifier:
                            dir_labels = dir_labels[top_scores_keep]
                        top_labels = top_labels[top_scores_keep]
                    boxes_for_nms = box_preds[:, [0, 1, 3, 4, 6]]
                    if not self._use_rotate_nms:
                        box_preds_corners = box_torch_ops.center_to_corner_box2d(
                            boxes_for_nms[:, :2], boxes_for_nms[:, 2:4],
                            boxes_for_nms[:, 4])
                        boxes_for_nms = box_torch_ops.corner_to_standup_nd(
                            box_preds_corners)
                    # the nms in 3d detection just remove overlap boxes.
                    selected = nms_func(
                        boxes_for_nms,
                        top_scores,
                        pre_max_size=self._nms_pre_max_size,
                        post_max_size=self._nms_post_max_size,
                        iou_threshold=self._nms_iou_threshold,
                    )
                else:
                    selected = None
                if selected is not None:
                    selected_boxes = box_preds[selected]
                    if self._use_direction_classifier:
                        selected_dir_labels = dir_labels[selected]
                    selected_labels = top_labels[selected]
                    selected_scores = top_scores[selected]
            # finally generate predictions.

            if selected_boxes is not None:
                box_preds = selected_boxes
                scores = selected_scores
                label_preds = selected_labels
                if self._use_direction_classifier:
                    dir_labels = selected_dir_labels
                    opp_labels = (box_preds[..., -1] > 0) ^ (dir_labels.byte() > 0)
                    box_preds[..., -1] += torch.where(
                        opp_labels,
                        torch.tensor(np.pi).type_as(box_preds),
                        torch.tensor(0.0).type_as(box_preds))
                    # box_preds[..., -1] += (
                    #     ~(dir_labels.byte())).type_as(box_preds) * np.pi
                final_box_preds = box_preds
                final_scores = scores
                final_labels = label_preds

                predictions_dict = (final_box_preds, final_scores, label_preds, img_idx)
            else:
                predictions_dict = (None, None, None, None, None, img_idx)
            # predictions_dicts.append(predictions_dict)
            predictions_dicts += (predictions_dict, )
        #self._total_postprocess_time += time.time() - t
        return predictions_dicts

def prepare_loss_weights(labels,
                         pos_cls_weight= 1.0,
                         neg_cls_weight= 1.0,
                         loss_norm_type= 'NormByNumPositives',
                         dtype= torch.float32):

    cared = labels >= 0
    #cared : [N,num_anchors]
    positives = labels >0
    negatives = labels ==0
    negatives_cls_weights = negatives.type(dtype) * neg_cls_weight
    cls_weights = neg_cls_weight + pos_cls_weight * positives.type(dtype)
    reg_weights = positives.type(dtype)
    if loss_norm_type == 'NormByNumExamples':
        num_examples = cared.type(dtype).sum(1, keepdim=True)
        num_examples = torch.clamp(num_examples, min=1.0)
        cls_weights /= num_examples
        bbox_normalizer = positives.sum(1, keepdim=True).type(dtype)
        reg_weights /= torch.clamp(bbox_normalizer, min=1.0)
    elif loss_norm_type == 'NormByNumPositives':  # for focal loss
        pos_normalizer = positives.sum(1, keepdim=True).type(dtype)
        reg_weights /= torch.clamp(pos_normalizer, min=1.0)
        cls_weights /= torch.clamp(pos_normalizer, min=1.0)
    elif loss_norm_type == 'NormByNumPosNeg':
        pos_neg = torch.stack([positives, negatives], dim=-1).type(dtype)
        normalizer = pos_neg.sum(1, keepdim=True)  # [N, 1, 2]
        cls_normalizer = (pos_neg * normalizer).sum(-1)  # [N, M]
        cls_normalizer = torch.clamp(cls_normalizer, min=1.0)
        # cls_normalizer will be pos_or_neg_weight/num_pos_or_neg
        normalizer = torch.clamp(normalizer, min=1.0)
        reg_weights /= normalizer[:, 0:1, 0]
        cls_weights /= cls_normalizer
    else:
        raise ValueError(
            f"unknown loss norm type. available: {list(LossNormType)}")
    return cls_weights, reg_weights, cared

def create_loss(loc_loss_ftor,
                cls_loss_ftor,
                box_preds,
                cls_preds,
                cls_targets,
                cls_weights,
                reg_targets,
                reg_weights,
                num_class,
                encode_background_as_zeros=True,
                encode_rad_error_by_sin=True,
                box_code_size = 7):
    batch_size = int(box_preds.shape[0])
    box_preds = box_preds.view(batch_size, -1, box_code_size)
    if encode_background_as_zeros:
        cls_preds = cls_preds.view(batch_size,-1,num_class)
    else:
        cls_preds = cls_preds.view(batch_size,-1,num_class+1)
    cls_targets = cls_targets.squeeze(-1)
    one_hot_targets = libs.tools.one_hot(
        cls_targets, depth=num_class+1,dtype=box_preds.dtype
    )
    if encode_background_as_zeros:
        one_hot_targets = one_hot_targets[..., 1:]
    if encode_rad_error_by_sin:
        # sin(a-b) = sina*cosb-cosa*sinb
        box_preds, reg_targets = add_sin_difference(box_preds, reg_targets)
    loc_losses = loc_loss_ftor(
        box_preds, reg_targets, weights=reg_weights)  # [N, M]
    cls_losses = cls_loss_ftor(
        cls_preds, one_hot_targets, weights=cls_weights)  # [N, M]
    return loc_losses, cls_losses

def add_sin_difference(boxes1, boxes2):
    rad_pred_encoding = torch.sin(boxes1[..., -1:]) * torch.cos(
        boxes2[..., -1:])
    rad_tg_encoding = torch.cos(boxes1[..., -1:]) * torch.sin(boxes2[..., -1:])
    boxes1 = torch.cat([boxes1[..., :-1], rad_pred_encoding], dim=-1)
    boxes2 = torch.cat([boxes2[..., :-1], rad_tg_encoding], dim=-1)
    return boxes1, boxes2

def _get_pos_neg_loss(cls_loss, labels):
    # cls_loss: [N, num_anchors, num_class]
    # labels: [N, num_anchors]
    batch_size = cls_loss.shape[0]
    if cls_loss.shape[-1] == 1 or len(cls_loss.shape) == 2:
        cls_pos_loss = (labels > 0).type_as(cls_loss) * cls_loss.view(
            batch_size, -1)
        cls_neg_loss = (labels == 0).type_as(cls_loss) * cls_loss.view(
            batch_size, -1)
        cls_pos_loss = cls_pos_loss.sum() / batch_size
        cls_neg_loss = cls_neg_loss.sum() / batch_size
    else:
        cls_pos_loss = cls_loss[..., 1:].sum() / batch_size
        cls_neg_loss = cls_loss[..., 0].sum() / batch_size
    return cls_pos_loss, cls_neg_loss

def get_direction_target(anchors, reg_targets, one_hot=True):
    batch_size = reg_targets.shape[0]
    anchors = anchors.view(batch_size, -1, 7)
    rot_gt = reg_targets[..., -1] + anchors[..., -1]
    dir_cls_targets = (rot_gt > 0).long()
    if one_hot:
        dir_cls_targets = libs.tools.one_hot(
            dir_cls_targets, 2, dtype=anchors.dtype)
    return dir_cls_targets


text:
import numpy as np
from si.util.util import euclidian_distance

__all__ = ['Kmeans']

class Kmeans:

	def __init__(self, k: int, iterations=100):
		self.k = k
		self.max_iter = iterations
		self.centroids = None
		self.distance = euclidian_distance

	def fit(self, dataset):
		x = dataset.X
		self._min = np.min(x, axis=0)
		self._max = np.max(x, axis=0)

	def generate_random_centroid(self, dataset):
		x = dataset.X
		self.centroids = np.array([np.random.uniform(low=self._min[i], high=self._max[i], size=(self.k,)) for i in range(x.shape[1])]).T

	def closest_centroid(self, x):
		dist = self.distance(x, self.centroids)
		closest_centroid = np.argmin(dist, axis=0)
		return closest_centroid

	def transform(self, dataset):
		self.generate_random_centroid(dataset)
		x = dataset.X
		changed = True
		count = 0
		old_idxs = np.zeros(x.shape[0])
		while changed and count < self.max_iter:
			self.idxs = np.apply_along_axis(self.closest_centroid, axis=0, arr=x.T)
			cent = list()
			for i in range(self.k):
				cent.append(np.mean(x[self.idxs == i], axis=0))
			self.centroids = np.array(cent)
			changed = np.any(old_idxs != self.idxs)
			old_idxs = self.idxs
			count += 1
		return self.centroids, self.idxs

	def fit_transform(self, dataset):
		self.fit(dataset)
		return self.transform(dataset)

ext: py

reason: import numpy


text:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
import sklearn
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression, Lasso, Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.datasets import load_boston, load_iris, load_wine, load_digits, \
    load_breast_cancer, load_diabetes, fetch_mldata
import statsmodels.api as sm
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
import itertools

np.set_printoptions(precision=3)

class PolyModel:
    def __init__(self, order=2, model_type=LinearRegression(normalize=True)):
        self.order = order
        self.model_type = model_type

    def predict(self, X):
        return self.model.predict(X[:,self.vars])

    def fit(self, X, y, vars):
        """
        Get nth order polynomial with interaction terms; X is full X,
        vars is tuple, list, or set like (1,5) to indicate columns of X.
        """
        self.vars = list(sorted(vars))
        X = X[:,self.vars]
        poly = PolynomialFeatures(self.order)
        model = make_pipeline(poly, sklearn.base.clone(self.model_type))
        model.fit(X, y)
        lr = model.named_steps[self.model_type.__class__.__name__.lower()]
        terms = poly.get_feature_names()
        terms = [t.replace('x1','x'+str(self.vars[1])) for t in terms]
        terms = [t.replace('x0','x'+str(self.vars[0])) for t in terms]
        terms[0] = 'c' # replace the '1' with 'c'
        terms = [f'{c:.2f}{t}' for c, t in zip(lr.coef_, terms)]
        self.eqn = ' + '.join(terms)
        self.model = model

    def __str__(self):
        return self.eqn
    def __repr__(self):
        return self.__str__()

class GreedyLayer:
    def __init__(self, k, order=2, metric=r2_score, model_type=LinearRegression(normalize=True)):
        self.k = k
        self.order = order
        self.metric = metric
        self.model_type = model_type
        self.models = []

    def predict(self, X): # X is full matrix, each model knows its var pair
        X_output = np.empty(shape=(X.shape[0],len(self.models)), dtype=float)
        for j,m in enumerate(self.models):
            X_output[:,j] = m.predict(X)
        return X_output

    def fit(self, X_train, y_train):
        p = X_train.shape[1]
        allvars = range(p)
        # if self.k == p:
        #     selected_vars = allvars
        # else: # pick subset of vars
        #     selected_vars = np.random.choice(allvars, min(p,self.k), replace=False)
        pairs = list(itertools.combinations(allvars, 2))
        for j, pair in enumerate(pairs):
            m = PolyModel(order=self.order, model_type=self.model_type)
            m.fit(X_train, y_train, vars=pair)
            self.models.append(m)

    def cull(self, idx):
        self.models = [self.models[i] for i in idx]

    def __str__(self):
        return '\n'.join([str(m) for m in self.models])

class GreedyNet:
    def __init__(self, n_layers=3, k=None, order=2, metric=r2_score, model_type=LinearRegression(normalize=True)):
        self.n_layers = n_layers
        self.k = k # None implies all n choose k combos of input parameters
        self.order = order
        self.metric = metric
        self.model_type=model_type
        self.layers = [] # list of GreedyLayer

    def predict(self, X):
        X_input = X
        for layer in self.layers:
            X_output = layer.predict(X_input)
            X_input = X_output
        return X_output

    def fit(self, X_train, y_train):
        X_input = X_train
        k = X_train.shape[1] if self.k is None else self.k
        for layer in range(self.n_layers):
            gl = GreedyLayer(k=k, order=self.order, model_type=self.model_type)
            gl.fit(X_input, y_train)
            self.layers.append(gl)
            Y_pred = gl.predict(X_input)
            # Y_test_pred = gl.predict(X_test) # move test cases along through the layer

            # print(f"LAYER {layer}\n",gl,"\nY_pred\n",Y_pred,"\nY_test_pred\n",X_test)
            train_metrics = []
            # test_metrics = []
            for mi,m in enumerate(gl.models):
                train_score = self.metric(y_train, Y_pred[:,mi])
                # test_score = self.metric(y_test, Y_test_pred[:,mi])
                train_metrics.append(train_score)
                # test_metrics.append(test_score)
                # print(f"pair {m.vars}", train_score, test_score)

            train_metrics = np.array(train_metrics)
            best_idx = np.array(list(reversed(np.argsort(train_metrics))))
            best_idx = best_idx[:k]
            print("BEST train", train_metrics[best_idx][0:5])

            gl.cull(best_idx)  # cull to top k per training metric
            X_input = Y_pred[:,best_idx]
            X_input = np.hstack([X_train,X_input]) # pack original vars in as well
            # X_test = Y_test_pred[:,best_idx]

            # test_metrics = np.array(test_metrics)
            # best_idx = np.array(list(reversed(np.argsort(test_metrics))))
            # best_idx = best_idx[:k]
            # print("BEST test", test_metrics[best_idx])

            # print(f"LAYER {layer}\n",gl)

def synthetic():
    np.random.seed(1)
    n = 20
    p = 3

    X = np.empty(shape=(n,p), dtype=float)
    for j in range(p):
        X[:,j] = (np.random.random(size=n)*10).astype(int)

    y = X[:,0] + X[:,1] + X[:,2]
    print("X\n",X)
    print("y\n", y.reshape(-1,1))

    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)

    net = GreedyNet(n_layers=5, model_type=Lasso(alpha=.01, normalize=True))
    net.fit(X_train, X_test, y_train, y_test)
    X_output = net.layers[0].predict(X_train)

def rent():
    np.random.seed(1)
    df = pd.read_csv("rent10k.csv")

    X = df.drop('price', axis=1).values
    # X = X[:,0:4].copy()
    y = df['price'].values
    # print("X\n",X)
    # print("y\n", y.reshape(-1,1))
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

    # net = GreedyNet(order=2, n_layers=100, k=25, model_type=Lasso(alpha=1, normalize=True))
    net = GreedyNet(order=3, n_layers=1500, k=15, model_type=LinearRegression(normalize=True))
    # net.fit(X_train, X_test, y_train, y_test)
    net.fit(X_train, y_train)
    return

    X_output = net.predict(X_train)
    y_pred = X_output[:, 0]
    print("BEST OUTPUT:", y_pred[:10]) # should be best
    print("y_train:", y_train[:10])
    print(r2_score(y_train, y_pred), mean_absolute_error(y_train, y_pred))

rent()

max_stars_repo_name: parrt/selfnet

max_issues_repo_name: parrt/selfnet

max_forks_repo_name: parrt/selfnet

reason: import numpy,import statsmodels


text:
"""
Wavesloshing Problem
"""
import numpy as np
from math import cos
from proteus import (Domain, Context,
                     FemTools as ft,
                     MeshTools as mt,
                     WaveTools as wt)
from proteus.mprans import SpatialTools as st
from proteus.Profiling import logEvent
from proteus.mprans.SpatialTools import Tank2D

# predefined options
opts=Context.Options([
    # water
    ("water_depth_fraction", 0.5, "Nondimensional water level normalised to the tank height"),
    ("amplitude", 0.005, "Nondimensional amplitude normalised to tank height"),
    ("nu",0.,"Kinematic viscosity m/sec^2. Set to 0 when comparing to analytical results for inviscid flow"),
    # tank
    ("tank_dim", (0.1 , 0.1), "Dimensions of the tank (length, height) in m"),
    #gravity
    ("g",(0,-9.81,0), "Gravity vector"),
    # gauges
    ("gauge_output", True, "Produce gauge data. A free-surface gauge is located at the far left boundary"),
    # mesh refinement and boundary condition
    ("he", 0.001,"Characteristic element size in m"),
    ("cfl", 0.5 ,"Target cfl"),
    ("opentop", True,"Enabling atmosphere boundary at the top"),
    # run time options
    ("T", 0.1 ,"Simulation time in s"),
    ("nsave", 0., "number of saves per second"),
    ("dt_init", 0.001 ,"Maximum initial time step in s"),
    ("dt_fixed", None, "fixed time step for proteus (scale with period) in s"),
    ("structured", False, "Use a (triangular) structured mesh"),
    ("useHex", False, "Use a hexahedral structured mesh"),
    ("gen_mesh", True ,"Generate new mesh"),
    ])

# ----- CONTEXT ------ #

# tank
tank_dim = opts.tank_dim

m = 1
km = m*np.pi/tank_dim[0]
water_amplitude = opts.amplitude
water_depth = h = tank_dim[1]*opts.water_depth_fraction

wavelength = tank_dim[0]*2.
k = 2*np.pi/wavelength
h = opts.water_depth_fraction*opts.tank_dim[1]
eps = k*water_amplitude

##########################################
#     Discretization Input Options       #
##########################################

#[temp] temporary location
backgroundDiffusionFactor = 0.01 #[temp]

genMesh = opts.gen_mesh
movingDomain = False
checkMass = False
applyRedistancing = True
useOldPETSc = False
useSuperlu = False
timeDiscretization = 'be'  # 'vbdf', 'be', 'flcbdf'
spaceOrder = 1
useHex = opts.useHex
structured = opts.structured
useRBLES = 0.0
useMetrics = 1.0
applyCorrection = True
useVF = 1.0
useOnlyVF = False
useRANS = 0  # 0 -- None
             # 1 -- K-Epsilon
             # 2 -- K-Omega

# ----- INPUT CHECKS ----- #
if spaceOrder not in [1,2]:
    raise ValueError("INVALID: spaceOrder(" + str(spaceOrder) + ")")

if useRBLES not in [0.0, 1.0]:
    raise ValueError("INVALID: useRBLES(" + str(useRBLES) + ")")

if useMetrics not in [0.0, 1.0]:
    raise ValueError("INVALID: useMetrics(" + str(useMetrics) + ")")

# ----- DISCRETIZATION ----- #
nd = 2
if spaceOrder == 1:
    hFactor=1.0
    if useHex:
        basis=ft.C0_AffineLinearOnCubeWithNodalBasis
        elementQuadrature = ft.CubeGaussQuadrature(nd,2)
        elementBoundaryQuadrature = ft.CubeGaussQuadrature(nd-1,2)
    else:
    	 basis=ft.C0_AffineLinearOnSimplexWithNodalBasis
         elementQuadrature = ft.SimplexGaussQuadrature(nd,3)
         elementBoundaryQuadrature = ft.SimplexGaussQuadrature(nd-1,3)
elif spaceOrder == 2:
    hFactor=0.5
    if useHex:
        basis=ft.C0_AffineLagrangeOnCubeWithNodalBasis
        elementQuadrature = ft.CubeGaussQuadrature(nd,4)
        elementBoundaryQuadrature = ft.CubeGaussQuadrature(nd-1,4)
    else:
        basis=ft.C0_AffineQuadraticOnSimplexWithNodalBasis
        elementQuadrature = ft.SimplexGaussQuadrature(nd,4)
        elementBoundaryQuadrature = ft.SimplexGaussQuadrature(nd-1,4)

##########################################
# Numerical Options and Other Parameters #
##########################################

weak_bc_penalty_constant = 100.0
nLevels = 1

# ----- PHYSICAL PROPERTIES ----- #

# Water
rho_0 = 998.2
nu_0 = opts.nu#0.#1.004e-6

# Air
rho_1 = 1.205
nu_1 = 0.#1.500e-5

# Surface Tension
sigma_01 = 0.0

# Gravity
g = opts.g

# ----- TIME STEPPING & VELOCITY----- #

T = opts.T
dt_init = opts.dt_init
runCFL = opts.cfl
nDTout = int(round(T*opts.nsave))
if nDTout > 0:
    dt_out= (T-dt_init)/nDTout
else:
    dt_out = 0
dt_fixed = opts.dt_fixed

# ----- DOMAIN ----- #

if useHex:
    nnx = 4 * refinement + 1
    nny=2*refinement+1
    hex=True
    domain = Domain.RectangularDomain(tank_dim)
elif structured:
    nnx = 4 * refinement
    nny = 2 * refinement
    domain = Domain.RectangularDomain(tank_dim)
    boundaryTags = domain.boundaryTags
else:
    domain = Domain.PlanarStraightLineGraphDomain()

# ----- TANK ----- #

tank = Tank2D(domain, tank_dim)
tank.facets = np.array([[[0, 1, 2, 3]]])
tank.facetFlags = np.array([1])

# ----- GAUGES ----- #
gauge_dx = tank_dim[0]/100.
probes=np.linspace(0., tank_dim[0], tank_dim[0]/gauge_dx+1)
PG=[]
PG2=[]
zProbes=water_depth*0.5
for i in probes:
    PG.append((i, zProbes, 0.),)
    PG2.append((i, water_depth, 0.),)

if opts.gauge_output:

    tank.attachPointGauges(
        'twp',
        gauges = ((('p',), PG),),
        activeTime=(0, opts.T),
        sampleRate=0,
        fileName='pointGauge_pressure.csv'
    )
    tank.attachPointGauges(
        'ls',
        gauges = ((('phi',), PG2),),
        activeTime=(0, opts.T),
        sampleRate=0,
        fileName='pointGauge_levelset.csv'
    )

# ----- EXTRA BOUNDARY CONDITIONS ----- #

if opts.opentop is True:
    tank.BC['y+'].setAtmosphere()
else:
    tank.BC['y+'].setFreeSlip()
tank.BC['y-'].setFreeSlip()
tank.BC['x+'].setFreeSlip()
tank.BC['x-'].setFreeSlip()

# ----- MESH CONSTRUCTION ----- #

he = opts.he
domain.MeshOptions.he = he
st.assembleDomain(domain)

# ----- STRONG DIRICHLET ----- #

ns_forceStrongDirichlet = False

# ----- NUMERICAL PARAMETERS ----- #
if useMetrics:
    ns_shockCapturingFactor  = 0.25
    ns_lag_shockCapturing = True
    ns_lag_subgridError = True
    ls_shockCapturingFactor  = 0.25
    ls_lag_shockCapturing = True
    ls_sc_uref  = 1.0
    ls_sc_beta  = 1.0
    vof_shockCapturingFactor = 0.25
    vof_lag_shockCapturing = True
    vof_sc_uref = 1.0
    vof_sc_beta = 1.0
    rd_shockCapturingFactor  = 0.25
    rd_lag_shockCapturing = False
    epsFact_density = epsFact_viscosity = epsFact_curvature \
        = epsFact_vof = ecH \
        = epsFact_consrv_dirac = epsFact_density \
        = 3.0
    epsFact_redistance = 0.33
    epsFact_consrv_diffusion = 0.1
    redist_Newton = False
    kappa_shockCapturingFactor = 0.25
    kappa_lag_shockCapturing = True#False
    kappa_sc_uref = 1.0
    kappa_sc_beta = 1.0
    dissipation_shockCapturingFactor = 0.25
    dissipation_lag_shockCapturing = True#False
    dissipation_sc_uref = 1.0
    dissipation_sc_beta = 1.0
else:
    ns_shockCapturingFactor  = 0.9
    ns_lag_shockCapturing = True
    ns_lag_subgridError = True
    ls_shockCapturingFactor  = 0.9
    ls_lag_shockCapturing = True
    ls_sc_uref  = 1.0
    ls_sc_beta  = 1.0
    vof_shockCapturingFactor = 0.9
    vof_lag_shockCapturing = True
    vof_sc_uref  = 1.0
    vof_sc_beta  = 1.0
    rd_shockCapturingFactor  = 0.9
    rd_lag_shockCapturing = False
    epsFact_density = epsFact_viscosity = epsFact_curvature \
        = epsFact_vof = ecH \
        = epsFact_consrv_dirac = epsFact_density \
        = 1.5
    epsFact_redistance = 0.33
    epsFact_consrv_diffusion = 10.0
    redist_Newton = False
    kappa_shockCapturingFactor = 0.9
    kappa_lag_shockCapturing = True#False
    kappa_sc_uref  = 1.0
    kappa_sc_beta  = 1.0
    dissipation_shockCapturingFactor = 0.9
    dissipation_lag_shockCapturing = True#False
    dissipation_sc_uref  = 1.0
    dissipation_sc_beta  = 1.0

# ----- NUMERICS: TOLERANCES ----- #
ns_nl_atol_res = max(1.0e-8,0.01*he**2)
vof_nl_atol_res = max(1.0e-8,0.01*he**2)
ls_nl_atol_res = max(1.0e-8,0.01*he**2)
rd_nl_atol_res = max(1.0e-8,0.01*he)
mcorr_nl_atol_res = max(1.0e-8,0.01*he**2)
kappa_nl_atol_res = max(1.0e-8,0.01*he**2)
dissipation_nl_atol_res = max(1.0e-8,0.01*he**2)

# ----- TURBULENCE MODELS ----- #
ns_closure=0 #1-classic smagorinsky, 2-dynamic smagorinsky, 3 -- k-epsilon, 4 -- k-omega
if useRANS == 1:
    ns_closure = 3
elif useRANS == 2:
    ns_closure = 4

##########################################
#                Solution                #
##########################################
# 3rd order solution

def eta0(x, t):
    eta_0 = np.sin(t)*np.cos(x)
    return eta_0

def phi0(x, y, t, h, w0):
    phi_0 = w0/(np.sinh(h))*np.cos(t)*np.cos(x)*np.cosh(y+h)
    return phi_0

def w0(h):
    w_0 = np.sqrt(np.tanh(h))
    return w_0

def eta1(x, t, w0):
    eta_1 = 1./8.*((w0**2-w0**(-2))+(w0**(-2)-3*w0**(-6))*np.cos(2.*t))*np.cos(2.*x)
    return eta_1

def phi1(x, y, t, h, w0):
    beta0 = 0  # arbitrary
    phi_1 = beta0+1./8.(w0-w0**(-3))*t-1./16.*(3*w0+w0**(-3))*np.sin(2*t)-3./(16.*np.cosh(2.*h))*(w0-w0**(-7))*np.sin(2*t)*np.cos(2*x)*np.cosh(2*(y+h))
    return phi_1

w1 = 0.

def eta2(x, t, w0):
    b11 = 1./32.*(3.*w0**(-8)+6.*w0**(-4)-5.+2.*w0**4)
    b13 = 3./128.*(9.*w0**(-8)+27.*w0**(-4)-15.+w0**4+2*w0**8)
    b31 = 1./128.*(3.*w0**(-8)+18.*w0**(-4)-5.)
    b33 = 3./128.*(-9.*w0**(-12)+3.*w0**(-8)-3.*w0**(-4)+1)
    eta_2 = b11*np.sin(t)*np.cos(x)+b13*np.sin(t)*np.cos(3*x)+b31*np.sin(3*t)*np.cos(x)+b33*np.sin(3*t)*np.cos(3*x)
    return eta_2

def phi2(x, y, t, h, w0):
    beta13 = 1./(128.*np.cosh(3.*h))*(1+3*w0**4)*(3*w0**(-9)-5.*w0**(-1)+2*w0**3)
    beta31 = 1./(128.*np.cosh(h))*(9.*w0**(-9)+62.*w0**(-5)-31.*w0**(-1))
    beta33 = 1./(128.*np.cosh(3.*h))*(1+3.*w0**4)*(-9.*w0**(-13)+22.*w0**(-9)-13.*w0**(-5))
    beta2 = 0.  # arbitrary
    phi_2 = beta2+beta13*np.cos(t)*np.cos(3*x)*np.cosh(3*(y+h))+beta31*np.cos(3*t)*np.cos(x)*np.cosh(y+h)+beta33*np.cos(3*t)*np.cos(3*x)*np.cosh(3*(y+h))
    return phi_2

def eps_eta(x, t, h, eps):
    w_0 = w0(h)
    eta_0 = eta0(x, t)
    eta_1 = eta1(x, t, w_0)
    eta_2 = eta2(x, t, w_0)
    epseta = eps*eta_0+eps**2*eta_1+0.5*eps**3*eta_2
    return epseta

def w2(w0):
    w_2 = 1./32.*(9.*w0**(-7)-12.*w0**(-3)-3*w0-2*w0**5)
    return w_2

def omega(h, eps):
    w_0 = w0(h)
    w_2 = w2(w_0)
    w = w_0+0.5*eps**2*w_2
    return w

def d_phi0_d_t(x, y, t, h, w0):
    return -w0/np.sinh(h)*np.sin(t)*np.cos(x)*np.cosh(y+h)

def d_phi0_d_x(x, y, t, h, w0):
    return -w0/np.sinh(h)*np.cos(t)*np.sin(x)*np.cosh(y+h)

def d_phi0_d_y(x, y, t, h, w0):
    return w0/np.sinh(h)*np.cos(t)*np.cos(x)*np.sinh(y+h)

def d_phi1_d_t(x, y, t, h, w0):
    return 1./8.*(w0-w0**(-3))-2./16.*(3.*w0+w0**(-3))*np.cos(2*t)-2.*3./(16.*np.cosh(2.*h))*(w0-w0**(-7))*np.cos(2.*t)*np.cos(2.*x)*np.cosh(2.*(y+h))

def d_phi1_d_x(x, y, t, h, w0):
    return 2.*3./(16.*np.cosh(2.*h))*(w0-w0**(-7))*np.sin(2.*t)*np.sin(2.*x)*np.cosh(2.*(y+h))

def d_phi1_d_y(x, y, t, h, w0):
    return 2.*3./(16.*np.cosh(2.*h))*(w0-w0**(-7))*np.sin(2.*t)*np.cos(2.*x)*np.sinh(2.*(y+h))

def d_phi2_d_t(x, y, t, h, w0):
    beta13 = 1./(128.*np.cosh(3.*h))*(1+3*w0**4)*(3*w0**(-9)-5.*w0**(-1)+2*w0**3)
    beta31 = 1./(128.*np.cosh(h))*(9.*w0**(-9)+62.*w0**(-5)-31.*w0**(-1))
    beta33 = 1./(128.*np.cosh(3.*h))*(1+3.*w0**4)*(-9.*w0**(-13)+22.*w0**(-9)-13.*w0**(-5))
    return -beta13*np.sin(t)*np.cos(3.*x)*np.cosh(3.*(y+h))-3.*beta31*np.sin(3.*t)*np.cos(x)*np.cosh(y+h)-3.*beta33*np.sin(3.*t)*np.cos(3.*x)*np.cosh(3.*(y+h))

def d_phi2_d_x(x, y, t, h, w0):
    beta13 = 1./(128.*np.cosh(3.*h))*(1+3*w0**4)*(3*w0**(-9)-5.*w0**(-1)+2*w0**3)
    beta31 = 1./(128.*np.cosh(h))*(9.*w0**(-9)+62.*w0**(-5)-31.*w0**(-1))
    beta33 = 1./(128.*np.cosh(3.*h))*(1+3.*w0**4)*(-9.*w0**(-13)+22.*w0**(-9)-13.*w0**(-5))
    return -3.*beta13*np.cos(t)*np.sin(3.*x)*np.cosh(3.*(y+h))-beta31*np.cos(3.*t)*np.sin(x)*np.cosh(y+h)-3.*beta33*np.cos(3.*t)*np.sin(3.*x)*np.cosh(3.*(y+h))

def d_phi2_d_y(x, y, t, h, w0):
    beta13 = 1./(128.*np.cosh(3.*h))*(1+3*w0**4)*(3*w0**(-9)-5.*w0**(-1)+2*w0**3)
    beta31 = 1./(128.*np.cosh(h))*(9.*w0**(-9)+62.*w0**(-5)-31.*w0**(-1))
    beta33 = 1./(128.*np.cosh(3.*h))*(1+3.*w0**4)*(-9.*w0**(-13)+22.*w0**(-9)-13.*w0**(-5))
    return 3.*beta13*np.cos(t)*np.cos(3.*x)*np.sinh(3.*(y+h))+beta31*np.cos(3.*t)*np.cos(x)*np.sinh(y+h)+3.*beta33*np.cos(3.*t)*np.cos(3.*x)*np.sinh(3.*(y+h))

def pressure(x, y, t, h, eps, rho, g, k, p0):
    x_ = x*k
    y_ = y*k
    h_ = h*k
    w_ = omega(h_, eps)
    t_ = (t+2*np.pi/(w_*np.sqrt(k*(-g[1])))*0.25)*(w_*np.sqrt(k*(-g[1])))
    w_0 = w0(h_)
    dt = eps*d_phi0_d_t(x_, y_, t_, h_, w_0)+eps**2*d_phi1_d_t(x_, y_, t_, h_, w_0)+0.5*eps**3*d_phi2_d_t(x_, y_, t_, h_, w_0)
    dx = eps*d_phi0_d_x(x_, y_, t_, h_, w_0)+eps**2*d_phi1_d_x(x_, y_, t_, h_, w_0)+0.5*eps**3*d_phi2_d_x(x_, y_, t_, h_, w_0)
    dy = eps*d_phi0_d_y(x_, y_, t_, h_, w_0)+eps**2*d_phi1_d_y(x_, y_, t_, h_, w_0)+0.5*eps**3*d_phi2_d_y(x_, y_, t_, h_, w_0)
    p = (-y_-dt*w_-0.5*dx**2-0.5*dy**2)*rho*abs(g[1])/k+p0
    return p

def eta(x, t):
    x_ = x*k
    h_ = h*k
    w_ = omega(h_, eps)
    t_ = (t+2*np.pi/(w_*np.sqrt(k*(-g[1])))*0.25)*(w_*np.sqrt(k*(-g[1])))
    eta = eps_eta(x_, t_, h_, eps)/k
    return eta

def signedDistance(x, t):
    #d=abs(x[1] - water_depth  - water_amplitude * cos(x[0]))
    d = x[1]-(water_depth+eta(x[0], t))
    return d

max_stars_repo_name: AntoninoScala/air-water-vv

max_issues_repo_name: AntoninoScala/air-water-vv

max_forks_repo_name: AntoninoScala/air-water-vv


text:
import numpy as np
from numbers import Number
from pytest import raises, warns
from hypothesis import given, strategies, unlimited
from hypothesis import settings as hyp_settings
from hypothesis import HealthCheck
from kernelmethods.numeric_kernels import DEFINED_KERNEL_FUNCS, PolyKernel, \
    GaussianKernel, LinearKernel, LaplacianKernel, SigmoidKernel, Chi2Kernel
from kernelmethods.utils import check_callable
from kernelmethods.base import KernelMatrix, KernelFromCallable, BaseKernelFunction
from kernelmethods.operations import is_positive_semidefinite

default_feature_dim = 10
range_feature_dim = [10, 50]
range_num_samples = [50, 100]

range_polynomial_degree = [2, 10] # degree=1 is tested in LinearKernel()

np.random.seed(42)

# choosing skip_input_checks=False will speed up test runs
# default values for parameters

num_tests_psd_kernel = 3

def gen_random_array(dim):
    """To better control precision and type of floats"""

    # TODO input sparse arrays for test
    return np.random.rand(dim)

def gen_random_sample(num_samples, sample_dim):
    """To better control precision and type of floats"""

    # TODO input sparse arrays for test
    return np.random.rand(num_samples, sample_dim)

def _test_for_all_kernels(kernel, sample_dim, check_PSDness=True):
    """Common tests that all kernels must pass."""

    x = gen_random_array(sample_dim)
    y = gen_random_array(sample_dim)

    try:
        result = kernel(x, y)
    except Exception:
        raise RuntimeError('{} unable to calculate!\n'
                           ' on x {}\n y{}'.format(kernel, x, y))

    if not isinstance(result, Number):
        raise ValueError('result {} of type {} is not a number!\n'
                         'x={}\ny={}\nkernel={}\n'
                         ''.format(result, type(result), x, y, kernel))

    if kernel(y, x) != result:
        raise ValueError('{} is not symmetric!'
                         'x={}\n y={}\n kernel={}\n'
                         ''.format(kernel.name, x, y, kernel))

    if check_PSDness:
        # ensuring it produces a PSD KM
        kernel.is_psd()

def test_kernel_design():
    """
    Every kernel must be
    1. must have a name defined
    2. must be callable with two samples
    3. returns a number

    """

    for kernel in DEFINED_KERNEL_FUNCS:

        # must be callable with 2 args
        check_callable(kernel, min_num_args=2)

        if not hasattr(kernel, 'name'):
            raise TypeError('{} does not have name attribute!'.format(kernel))

        # only numeric data is accepted and other dtypes must raise an error
        for non_num in ['string',
                        [object, object] ]:
            with raises(ValueError):
                _ = kernel(non_num, non_num)

def _test_func_is_valid_kernel(kernel, sample_dim, num_samples):
    """A func is a valid kernel if the kernel matrix generated by it is PSD.

    Not including this in tests for all kernels to allow for non-PSD kernels in the future

    """

    KM = KernelMatrix(kernel, name='TestKM')
    KM.attach_to(gen_random_sample(num_samples, sample_dim))
    is_psd = is_positive_semidefinite(KM.full, verbose=True)
    if not is_psd:
        raise ValueError('{} is not PSD'.format(str(KM)))

@hyp_settings(max_examples=num_tests_psd_kernel, deadline=None,
              timeout=unlimited, suppress_health_check=HealthCheck.all())
@given(strategies.integers(range_feature_dim[0], range_feature_dim[1]),
       strategies.integers(range_num_samples[0], range_num_samples[1]),
       strategies.integers(range_polynomial_degree[0], range_polynomial_degree[1]),
       strategies.floats(min_value=0, max_value=1e3,
                         allow_nan=False, allow_infinity=False))
def test_polynomial_kernel(sample_dim, num_samples,
                           poly_degree, poly_intercept):
    """Tests specific for Polynomial kernel."""

    poly = PolyKernel(degree=poly_degree, b=poly_intercept, skip_input_checks=False)
    _test_for_all_kernels(poly, sample_dim)
    _test_func_is_valid_kernel(poly, sample_dim, num_samples)

@hyp_settings(max_examples=num_tests_psd_kernel, deadline=None,
              timeout=unlimited, suppress_health_check=HealthCheck.all())
@given(strategies.integers(range_feature_dim[0], range_feature_dim[1]),
       strategies.integers(range_num_samples[0], range_num_samples[1]),
       strategies.floats(min_value=0, max_value=1e6,
                         allow_nan=False, allow_infinity=False))
def test_gaussian_kernel(sample_dim, num_samples, sigma):
    """Tests specific for Gaussian kernel."""

    gaussian = GaussianKernel(sigma=sigma, skip_input_checks=False)
    _test_for_all_kernels(gaussian, sample_dim)
    _test_func_is_valid_kernel(gaussian, sample_dim, num_samples)

@hyp_settings(max_examples=num_tests_psd_kernel, deadline=None,
              timeout=unlimited, suppress_health_check=HealthCheck.all())
@given(strategies.integers(range_feature_dim[0], range_feature_dim[1]),
       strategies.integers(range_num_samples[0], range_num_samples[1]))
def test_linear_kernel(sample_dim, num_samples):
    """Tests specific for Linear kernel."""

    linear = LinearKernel(skip_input_checks=False)
    _test_for_all_kernels(linear, sample_dim)
    _test_func_is_valid_kernel(linear, sample_dim, num_samples)

@hyp_settings(max_examples=num_tests_psd_kernel, deadline=None,
              timeout=unlimited, suppress_health_check=HealthCheck.all())
@given(strategies.integers(range_feature_dim[0], range_feature_dim[1]),
       strategies.integers(range_num_samples[0], range_num_samples[1]),
       strategies.floats(min_value=0, max_value=1e6,
                         allow_nan=False, allow_infinity=False))
def test_laplacian_kernel(sample_dim, num_samples, gamma):
    """Tests specific for Laplacian kernel."""

    laplacian = LaplacianKernel(gamma=gamma, skip_input_checks=False)
    _test_for_all_kernels(laplacian, sample_dim)
    _test_func_is_valid_kernel(laplacian, sample_dim, num_samples)

@hyp_settings(max_examples=num_tests_psd_kernel, deadline=None,
              timeout=unlimited, suppress_health_check=HealthCheck.all())
@given(strategies.integers(range_feature_dim[0], range_feature_dim[1]),
       strategies.integers(range_num_samples[0], range_num_samples[1]),
       strategies.floats(min_value=0, max_value=1e6,
                         allow_nan=False, allow_infinity=False),
       strategies.floats(min_value=0, max_value=1e6,
                         allow_nan=False, allow_infinity=False)
       )
def test_sigmoid_kernel(sample_dim, num_samples, gamma, offset):
    """Tests specific for sigmoid kernel."""

    sigmoid = SigmoidKernel(gamma=gamma, offset=offset, skip_input_checks=False)
    # sigmoid is not always PSD
    _test_for_all_kernels(sigmoid, sample_dim, check_PSDness=False)

@hyp_settings(max_examples=num_tests_psd_kernel, deadline=None,
              timeout=unlimited, suppress_health_check=HealthCheck.all())
@given(strategies.integers(range_feature_dim[0], range_feature_dim[1]),
       strategies.integers(range_num_samples[0], range_num_samples[1]),
       strategies.floats(min_value=0, max_value=1e6,
                         allow_nan=False, allow_infinity=False))
def test_chi2_kernel(sample_dim, num_samples, gamma):
    """Tests specific for Laplacian kernel."""

    chi2 = Chi2Kernel(gamma=gamma, skip_input_checks=False)
    _test_for_all_kernels(chi2, sample_dim)
    _test_func_is_valid_kernel(chi2, sample_dim, num_samples)

max_stars_repo_name: andreas-koukorinis/kernelmethods

max_issues_repo_name: andreas-koukorinis/kernelmethods

max_forks_repo_name: andreas-koukorinis/kernelmethods


text:
# Como executar:
# $ python maxmin.py <img_entrada> <img_saida_min> <imgsaida_max> <mask_size>
# <mask_size> é um número inteiro. Exemplo: Se mask_size=3 então a máscara possui tamanho 3x3.
# Gera duas imagens de saída.

import sys
import matplotlib.pyplot as plt
from scipy import misc
from skimage import img_as_float
from scipy.ndimage import filters

def loadImg(arg):
    return misc.imread(arg)

img_1 = loadImg(sys.argv[1])
saida_1 = sys.argv[2]+'.tif'
saida_2 = sys.argv[3]+'.tif'
mask_size = int(sys.argv[4])

# Converte os pixels em float, com valores entre 0 e 1
img_1 = img_as_float(img_1)

# Aplica o filtro de mínimo
img_saida_min = filters.minimum_filter(img_1, size = mask_size, mode = 'constant', cval = 0)

# Aplica o filtro de máximo
img_saida_max = filters.maximum_filter(img_1, size = mask_size, mode = 'constant')

# Faz o salvamento das imagens de saída após o processamento
misc.imsave(saida_1, img_saida_min)
misc.imsave(saida_2, img_saida_max)

# Organiza o plote das imagens
plt.figure()
plt.subplot(221);
plt.imshow(img_1, cmap='gray', interpolation='nearest');
plt.title('Imagem Entrada')
plt.subplot(222);
plt.imshow(img_saida_min, cmap='gray', interpolation='nearest')
plt.title('Imagem Saída Min.')
plt.subplot(223);
plt.imshow(img_saida_max, cmap='gray', interpolation='nearest')
plt.title('Imagem Saída Máx.')

# Plota as imagens de entrada e saída na tela
plt.tight_layout()
plt.show()

max_stars_repo_path: 3 - Spatial filtering for smoothing/4 - Max Min Filter/maxmin.py

reason: from scipy


text:
import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt
from scipy import io
from sklearn.model_selection import train_test_split
from tensorflow.contrib.distributions import Normal, Gamma
from stein.samplers import SteinSampler
from stein.optimizers import AdamGradientDescent

# Generate synthetic data.
X_train = np.random.uniform(size=(20, 1))
y_train = np.random.normal(np.cos(10 * X_train) * (5 * X_train), 0.1)

# Parameters for training such as the number of hidden neurons and the batch
# size to use during training, the total number of training iterations, and the
# number of particles to sample from the posterior.
n_hidden = 100
n_batch = 20
n_particles = 20
# Number of training data points.
n_train, n_feats = X_train.shape
# Extract learning rate.
learning_rate = 1e-1
# Precision prior parameters.
alpha, beta = 1., 0.01

# Define a Bayesian neural network model in TensorFlow.
with tf.variable_scope("model"):
    # Placeholder variables for data.
    model_X = tf.placeholder(tf.float32, shape=[None, n_feats])
    model_y = tf.placeholder(tf.float32, shape=[None, 1])
    # TensorFlow variables for the networks weights and biases, the precision
    # over the weights and biases, and the noise precision of the output.
    model_log_lambda = tf.Variable(tf.zeros([]))
    model_log_gamma = tf.Variable(tf.zeros([]))
    model_lambda = tf.exp(model_log_lambda)
    model_gamma = tf.exp(model_log_gamma)
    model_w_1 = tf.Variable(tf.zeros([n_feats, n_hidden]))
    model_b_1 = tf.Variable(tf.zeros([n_hidden]))
    model_w_2 = tf.Variable(tf.zeros([n_hidden, 1]))
    model_b_2 = tf.Variable(tf.zeros([]))

    # Compute the prediction from the network.
    with tf.variable_scope("prediction"):
        pred = tf.matmul(
            tf.nn.relu(tf.matmul(model_X, model_w_1) + model_b_1), model_w_2
        ) + model_b_2
    # Likelihood function.
    with tf.variable_scope("likelihood"):
        log_l_dist = Normal(pred, tf.reciprocal(tf.sqrt(model_gamma)))
        log_l = tf.reduce_sum(log_l_dist.log_prob(model_y))
    # Priors.
    with tf.variable_scope("priors"):
        prior_lambda = Gamma(alpha, beta)
        prior_gamma = Gamma(alpha, beta)
        prior_w_1 = Normal(
            tf.zeros([n_feats, n_hidden]),
            tf.reciprocal(tf.sqrt(model_lambda))
        )
        prior_b_1 = Normal(
            tf.zeros([n_hidden]),
            tf.reciprocal(tf.sqrt(model_lambda))
        )
        prior_w_2 = Normal(
            tf.zeros([n_hidden, 1]),
            tf.reciprocal(tf.sqrt(model_lambda))
        )
        prior_b_2 = Normal(
            tf.zeros([]),
            tf.reciprocal(tf.sqrt(model_lambda))
        )
    # Compute the log-posterior distribution.
    log_p = (
        # Rescaled log-likelihood (to account for batch updates).
        log_l * n_train / n_batch +
        # Variance priors.
        prior_lambda.log_prob(model_lambda) +
        prior_gamma.log_prob(model_gamma) +
        # Weight and bias priors.
        tf.reduce_sum(prior_w_1.log_prob(model_w_1)) +
        tf.reduce_sum(prior_w_2.log_prob(model_w_2)) +
        tf.reduce_sum(prior_b_1.log_prob(model_b_1)) +
        prior_b_2.log_prob(model_b_2)
    ) / n_train

# Gradient descent object.
gd = AdamGradientDescent(learning_rate=learning_rate, decay=0.999)
# Perform Stein variational gradient descent to sample from the posterior
# distribution of the Bayesian neural network.
sampler = SteinSampler(n_particles, log_p, gd)
# Number of iterations to perform before displaying diagnostics.
n_prog = 1000
# Perform learning iterations.
for i in range(10000):
    batch = np.random.choice(n_train, n_batch, replace=False)
    X, y = X_train[batch], y_train[batch]
    sampler.train_on_batch({model_X: X, model_y: y})
    if i % n_prog == 0:
        y_hat = sampler.function_posterior(pred, {model_X: X_train})
        mse = np.mean((y_train.ravel() - y_hat.mean(axis=0)) ** 2)
        print("Iteration: {}. Mean squared error: {:.4f}".format(i, mse))

if True:
    r = np.atleast_2d(np.linspace(0., 1.5, num=200)).T
    y_vis = sampler.function_posterior(pred, {model_X: r})
    plt.figure()
    plt.plot(X_train.ravel(), y_train.ravel(), "r.")
    for i in range(n_particles):
        plt.plot(r.ravel(), y_vis[i], "b-", alpha=0.3)
    plt.grid()
    plt.show()

max_stars_repo_name: JamesBrofos/Stein

max_issues_repo_name: JamesBrofos/Stein

max_forks_repo_name: JamesBrofos/Stein

reason: import numpy,from scipy


text:
import numpy as np
import time
import sys
from collections import deque
from agent.maddpg import MADDPG
from agent_examples import maddpg_agent_config_default, maddpg_agent_config_per
from utils import (
    log_training_info,
    save_checkpoint,
    get_training_env,
    load_agent_parameters,
)

def maddpg_runner(
    env,
    brain_name="Reacher",
    n_episodes=100,
    agent_config={},
    target_episodes=100,
    target_score=0.5,
    num_agents=20,
    print_log=False,
    run_mode="train",
    saved_checkpoint_path="model/checkpoint.pth",
):
    """
    Deep Deterministic Policy Gradients (DDPG).
    Parameters
    ----------
    env :
    brain_name :
    n_episodes :
    agent_config :
    target_episodes :
    target_score :
    num_agents :
    print_log :
    run_mode:
    saved_checkpoint_path:
    Returns
    -------

    """
    saved_good = False
    scores = []  # episodic scores
    training_info_log = []  # training time and meta data logger
    moving_avg = deque(maxlen=100)  # last 100 scores
    agent = MADDPG(**agent_config)

    if run_mode == "test":
        load_agent_parameters(agent, saved_checkpoint_path)
        for a in agent.agents:
            a.actor_local.eval()
            a.critic_local.eval()
    ## Perform n_episodes of training
    training_start_time = time.time()
    BETA_START = 0.5
    for i_episode in range(1, n_episodes + 1):
        beta = min(1.0, BETA_START + i_episode * (1.0 - BETA_START) / 1e6)
        states = env.reset(train_mode=run_mode == "train")[
            brain_name
        ].vector_observations
        scores_episode = np.zeros(num_agents)  # rewards per episode for each agent
        episode_start_time = time.time()
        agent.reset()
        while True:
            # Perform a step: S;A;R;S'
            actions = agent.act(states)  # select the next action for each agent
            env_info = env.step(actions)[
                brain_name
            ]  # send the actions to the environment
            rewards = env_info.rewards  # get the rewards
            next_states = env_info.vector_observations  # get the next states
            dones = env_info.local_done  # see if episode has finished
            # Send the results to the Agent
            if run_mode == "train":
                agent.step(states, actions, rewards, next_states, dones, beta)
            # Update the variables for the next iteration
            states = next_states
            scores_episode += rewards
            # break if any agents are done
            if np.any(dones):
                break

        if run_mode == "test":
            continue
        # Store the rewards and calculate the moving average
        score = np.max(scores_episode)
        scores.append(score)
        moving_avg.append(score)

        training_info_log.append(
            log_training_info(
                i_episode,
                training_start_time,
                episode_start_time,
                scores,
                moving_avg,
                print_log=print_log,
            )
        )
        ## Check if the environment has been solved
        # if np.mean(moving_avg) >= target_score and i_episode >= target_episodes:
        if i_episode % target_episodes == 0:
            save_checkpoint(
                i_episode,
                scores,
                moving_avg,
                training_info_log,
                agent,
                agent_config,
                training_start_time,
                target_episodes,
            )
        if np.mean(moving_avg) >= target_score and not saved_good:

            save_checkpoint(
                i_episode,
                scores,
                moving_avg,
                training_info_log,
                agent,
                agent_config,
                training_start_time,
                target_episodes,
            )
            saved_good = True
            print("\n checkpoint at episode {} saved.".format(i_episode))
            # break

    return scores, moving_avg, agent, training_info_log

if __name__ == "__main__":
    if not len(sys.argv) >= 5:
        print(
            "please path the path to environment file, n_episodes, max_t, print_log !"
        )
        exit(0)

    env_file_path = sys.argv[1]
    n_episodes = int(sys.argv[2])
    print_log = bool(int(sys.argv[3]))
    run_mode = sys.argv[4]
    saved_checkpoint_path = ""
    if run_mode == "test":
        saved_checkpoint_path = sys.argv[5]
    env, num_agents, brain_name, state_size, action_size = get_training_env(
        env_file_path, test_mode=run_mode == "test"
    )
    maddpg_agent_config_default["action_size"] = action_size
    maddpg_agent_config_default["state_size"] = state_size
    scores, moving_avg, agent, training_info_log = maddpg_runner(
        env=env,
        brain_name=brain_name,
        n_episodes=n_episodes,
        agent_config=maddpg_agent_config_default,
        num_agents=num_agents,
        print_log=print_log,
        run_mode=run_mode,
        saved_checkpoint_path=saved_checkpoint_path,
    )


text:
import numpy as np
import pyfits
from matplotlib import rc
from matplotlib import cm
import matplotlib as mpl
rc('font', family='serif')
rc('text', usetex=True)
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

print("reading in apogee labels")
# read in apogee labels

tr_IDs = np.load("../tr_id.npz")['arr_0']
labels_apogee = np.load("../tr_label.npz")['arr_0']
apogee_teff = labels_apogee[:,0]
apogee_logg = labels_apogee[:,1]
apogee_feh = labels_apogee[:,2]

# read in lamost labels

print("reading in lamost labels")
a = pyfits.open("../../make_lamost_catalog/lamost_catalog_training.fits")
b = a[1].data
a.close()
IDs_lamost = b['lamost_id']
IDs_lamost = np.array([val.strip() for val in IDs_lamost])
teff_all_lamost = b['teff_1']
logg_all_lamost = b['logg_1']
feh_all_lamost = b['feh']
inds = np.array([np.where(IDs_lamost==a)[0][0] for a in tr_IDs])
lamost_teff = teff_all_lamost[inds]
lamost_logg = logg_all_lamost[inds]
lamost_feh = feh_all_lamost[inds]

# histogram of LAMOST values and APOGEE values for Teff, logg, Fe/H

fig, (ax0,ax1,ax2) = plt.subplots(ncols=3, figsize=(10,4),
                              sharex=False, sharey=False)
plt.subplots_adjust(wspace=0.3)
bins = 50
col_apogee = 'darkorchid'
col_lamost = 'darkorange'
ax0.hist(
    lamost_teff, color=col_lamost, histtype='step',
    bins=bins, range=(3500,6000), normed=True, label="LAMOST")
ax0.hist(
    apogee_teff, color=col_apogee, histtype='step',
    bins=bins, range=(3500,6000), normed=True, label="APOGEE")
ax0.set_xlabel("${\mbox{T}_{\mbox{eff}}}$ [K]", fontsize=14)
ax0.legend(loc="upper left")
ax0.set_ylabel("Num. Training Objects, Normalized", fontsize=14)
ax0.set_ylim(0,0.003)
ax0.locator_params(nbins=5)
ax0.tick_params(axis='both', labelsize=14)
ax1.hist(
    lamost_logg, color=col_lamost, histtype='step',
    bins=bins, range=(0.2,4.8), normed=True, label="LAMOST")
ax1.hist(
    apogee_logg, color=col_apogee, histtype='step',
    bins=bins, range=(0.2,4.8), normed=True, label="APOGEE")
ax1.legend(loc="upper left")
ax1.set_ylim(0,2)
ax1.locator_params(nbins=5)
ax1.set_xlabel("log g [dex]", fontsize=14)
ax1.tick_params(axis='both', labelsize=14)
ax2.hist(
    lamost_feh, color=col_lamost, histtype='step',
    bins=bins, normed=True, range=(-2,0.5), label="LAMOST")
ax2.hist(
    apogee_feh, color=col_apogee, histtype='step',
    bins=bins, normed=True, range=(-2, 0.5), label="APOGEE")
ax2.legend(loc="upper left")
ax2.set_xlabel("[Fe/H] [dex]")
ax2.locator_params(nbins=5)
ax2.set_ylim(0,2)
ax2.tick_params(axis='both', labelsize=14)

#plt.show()
plt.tight_layout()
plt.savefig("hist_label_dist.png")
plt.close()

max_stars_repo_name: aceilers/TheCannon

max_issues_repo_name: aceilers/TheCannon

max_forks_repo_name: aceilers/TheCannon


text:
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.autograd import Variable
from torch import optim

from tempfile import gettempdir
from scipy.stats import pearsonr
from tqdm import tqdm

from torch_utils import *
from reader_utils import *
from skempi_lib import *
from grid_utils import *

VERBOSE = True
USE_CUDA = True
BATCH_SIZE = 4
LR = 0.01

def prepare_siamese_batch(batch):
    v1, v2, oh, ddg = zip(*batch)
    if USE_CUDA:
        dat_var1 = Variable(torch.FloatTensor(v1)).cuda()
        dat_var2 = Variable(torch.FloatTensor(v2)).cuda()
        ddg_var = Variable(torch.FloatTensor(ddg)).cuda()
        oh_var = Variable(torch.FloatTensor(oh)).cuda()
    else:
        dat_var1 = Variable(torch.FloatTensor(v1))
        dat_var2 = Variable(torch.FloatTensor(v2))
        ddg_var = Variable(torch.FloatTensor(ddg))
        oh_var = Variable(torch.FloatTensor(oh))
    return dat_var1, dat_var2, oh_var, ddg_var

def prepare_batch(batch):
    vox, oh, ddg = zip(*batch)
    if USE_CUDA:
        dat_var = Variable(torch.FloatTensor(vox)).cuda()
        ddg_var = Variable(torch.FloatTensor(ddg)).cuda()
        oh_var = Variable(torch.FloatTensor(oh)).cuda()
    else:
        dat_var = Variable(torch.FloatTensor(vox))
        ddg_var = Variable(torch.FloatTensor(ddg))
        oh_var = Variable(torch.FloatTensor(oh))
    return dat_var, oh_var, ddg_var

def batch_generator(loader, prepare_batch, batch_size=BATCH_SIZE):

    stop = False
    while not stop:
        batch = []
        while len(batch) < batch_size:
            try:
                batch.append(next(loader))
            except StopIteration:
                stop = True
                break
        if len(batch) == 0:
            break
        yield prepare_batch(batch)

class CNN3d(nn.Module):

    def __init__(self, dropout=0.2):
        super(CNN3d, self).__init__()

        self.features = nn.Sequential(
            nn.Conv3d(8, 64, kernel_size=(5, 5, 5)),
            nn.ReLU(inplace=True),
            nn.MaxPool3d((2, 2, 2)),
            nn.Conv3d(64, 128, kernel_size=(3, 3, 3)),
            nn.ReLU(inplace=True),
            nn.MaxPool3d((2, 2, 2)),
            nn.Conv3d(128, 256, kernel_size=(3, 3, 3)),
            nn.ReLU(inplace=True),
        )
        self.compression = nn.Sequential(
            nn.Linear(256, 100),
            nn.ReLU(inplace=True),
            nn.Dropout(dropout),
        )

        self.fc = nn.Sequential(
            nn.Linear(120, 1),
            nn.Dropout(dropout),
        )

    def compress(self, x):
        x = self.features(x)
        x = x.view(len(x), -1)
        x = self.compression(x)
        return x

    def forward(self, x, oh):
        x = self.compress(x)
        x = torch.cat([x, oh], 1)
        x = self.fc(x)
        return x.view(-1)

def eval_model(model, variables):
    vox, oh, ddg = variables
    ddg = ddg.view(-1)
    ddg_hat = model(vox, oh)
    loss = nn.MSELoss().cuda()(ddg_hat, ddg) if USE_CUDA else nn.MSELoss()(ddg_hat, ddg)
    return loss, ddg, ddg_hat

def evaluate(model, batch_generator, length_xy, batch_size=BATCH_SIZE):
    model.eval()
    pbar = tqdm(total=length_xy, desc="calculation...")
    err, i = 0, 0
    preds = np.zeros((length_xy,))
    truth = np.zeros((length_xy,))
    for i, vars in enumerate(batch_generator):
        loss, ddg, ddg_hat = eval_model(model, vars)
        err += loss.item()
        start = i * batch_size
        end = i * batch_size + len(ddg)
        preds[start:end] = ddg_hat.data.cpu().numpy()
        truth[start:end] = ddg.data.cpu().numpy()
        pbar.update(len(ddg))
        e = err / (i + 1)
        pbar.set_description("Validation Loss:%.6f" % (e,))
    pbar.close()
    return e, preds, truth

def train(model, opt, adalr, batch_generator, length_xy):
    model.train()
    pbar = tqdm(total=length_xy, desc="calculating...")
    err = 0
    for i, vars in enumerate(batch_generator):
        opt.zero_grad()
        loss, ddg, ddg_hat = eval_model(model, vars)
        adalr.update(loss.item())
        err += loss.item()
        loss.backward()
        opt.step()
        lr, e = adalr.lr, err/(i + 1)
        pbar.set_description("Training Loss:%.6f, LR: %.6f" % (e, lr))
        pbar.update(len(ddg))
    pbar.close()

def add_arguments(parser):
    parser.add_argument('-r', '--resume', default='', type=str, metavar='PATH',
                        help='path to latest checkpoint (default: none)')
    parser.add_argument("-e", "--eval_every", type=int, default=1,
                        help="How often to evaluate on the validation set.")
    parser.add_argument('-n', "--num_epochs", type=int, default=20,
                        help="How many epochs to train the model?")
    parser.add_argument("-o", "--out_dir", type=str, required=False,
                        default=gettempdir(), help="Specify the output directory.")
    parser.add_argument("-s", "--seed", type=int, default=9898,
                        help="Sets the seed for generating random number.")
    parser.add_argument("-d", "--device", type=str, choices=["0", "1", "2", "3"],
                        default="3", help="Choose a device to run on.")
    parser.add_argument("-g", '--debug', action='store_true', default=True,
                        help="Run in debug mode.")

if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser()
    add_arguments(parser)
    args = parser.parse_args()

    os.environ["CUDA_DEVICE_ORDER"] = "PCI_BUS_ID"
    os.environ["CUDA_VISIBLE_DEVICES"] = args.device

    net = CNN3d()

    # opt = optim.SGD(net.parameters(), lr=LR, momentum=0.9, nesterov=True)
    opt = optim.Adam(net.parameters(), lr=LR)

    if os.path.isfile(args.resume):
        print("=> loading checkpoint '%s'" % args.resume)
        checkpoint = torch.load(args.resume, map_location=lambda storage, loc: storage)
        init_epoch = checkpoint['epoch']
        loaded_dict = checkpoint['net']
        net.load_state_dict(loaded_dict, strict=False)
        opt.load_state_dict(checkpoint['opt'])
    else:
        print("=> no checkpoint found at '%s'" % args.resume)

    ckptpath = args.out_dir
    model_summary(net)

    non_blocking_producer = non_blocking_interface_producer

    indx1 = skempi_df.Protein.isin(G1 + G2 + G3 + G4 + G5)
    indx2 = (skempi_df.num_muts == 1)
    training_set = skempi_df.Protein[indx1 & indx2]
    validation_set = skempi_df.Protein[~indx1 & indx2]
    records = load_skempi(skempi_df[indx1].reset_index(), PDB_PATH, True, False)
    ddg_mapper = DDGMapper(records)
    rotations = get_xyz_rotations(0.5)
    # rotations = [None]
    reader_trn = SkempiReader(non_blocking_producer, training_set, ddg_mapper, rotations, num_voxels=20)
    loader_trn = PdbLoader(reader_trn, 2000, 19 * len(rotations))
    reader_val = SkempiReader(non_blocking_producer, validation_set, ddg_mapper, [None], num_voxels=20)
    loader_val = PdbLoader(reader_val, len(validation_set), 19)

    init_epoch = 0
    num_epochs = args.num_epochs

    adalr = AdaptiveLR(opt, LR, num_iterations=2000)

    # Move models to GPU
    if USE_CUDA:
        net = net.cuda()
    if USE_CUDA and args.resume:
        optimizer_cuda(opt)

    for epoch in range(init_epoch, num_epochs):

        train(net, opt, adalr, batch_generator(loader_trn, prepare_batch), len(loader_trn))

        if epoch < num_epochs - 1 and epoch % args.eval_every != 0:
            continue

        loss, preds, truth = evaluate(net, batch_generator(loader_val, prepare_batch), len(loader_val))
        cor, _ = pearsonr(preds, truth)

        if VERBOSE:
            print("[Epoch %d/%d] (Validation Loss: %.6f, PCC: %.6f" %
                  (epoch + 1, num_epochs, loss, cor))

        save_checkpoint({
            'lr': adalr.lr,
            'epoch': epoch,
            'net': net.state_dict(),
            'opt': opt.state_dict()
        }, loss, "skempi", ckptpath)

        loader_val.reset()
        loader_trn.reset()

max_stars_repo_name: yotamfr/skempi

max_issues_repo_name: yotamfr/skempi

max_forks_repo_name: yotamfr/skempi


text:
import numpy as np
import pandas as pd
import tempfile
from cycif_db.markers import Markers, Marker_Comparator
from cycif_db.model import Marker

cyc_markers = Markers()

def test_load_known_markers():
    markers, other_features = cyc_markers.markers, cyc_markers.other_features

    assert len(markers) > 80, len(markers)
    assert len(other_features) == 12, len(other_features)

    assert 'asma' in markers
    assert markers['asma'] == 268, markers['asma']

    assert 'area' in other_features
    assert other_features['area'] == 'area'

def test_get_marker_db_entry():
    assert cyc_markers.get_marker_db_entry('aSMA') == \
        ('aSMA', '', '', ''), \
        cyc_markers.get_marker_db_entry('aSMA')
    assert cyc_markers.get_marker_db_entry('alpha-SMA') == \
        ('aSMA', '', '', ''), \
        cyc_markers.get_marker_db_entry('alpha-SMA')

def test_get_other_feature_db_name():
    assert cyc_markers.get_other_feature_db_name('Area') == 'area'
    assert cyc_markers.get_other_feature_db_name('cellID') == 'sample_cell_id'
    assert cyc_markers.get_other_feature_db_name('X_centroid') == \
        'x_centroid'

def test_update_stock_markers():
    cyc_markers = Markers()
    o_markers = cyc_markers.markers
    o_other_features = cyc_markers.other_features

    new_markers = [
        ('ABC1', None, None, 1, 'ABC_1'),
        ('ABC1', 'eF660', None, None, 'ABC1_660'),
        ('ABC_2', None, None, None, 'ABC-2, ABC2'),
        ('CK14', None, None, None, 'CK14a'),
        ('CD3', 'eF450', 'goat', None, 'CD3a, CD3b')
    ]

    with tempfile.NamedTemporaryFile() as tmp:
        to_path = tmp.name
        cyc_markers.update_stock_markers(
            new_markers, toplace=to_path, reload=True)
        df = pd.read_csv(to_path, sep='\t').fillna('')

    other_features = cyc_markers.other_features

    assert cyc_markers._path == to_path
    assert len(other_features) == len(o_other_features), other_features

    assert len(cyc_markers.markers) - len(o_markers) == 7, \
        set(cyc_markers.markers.keys()) - set(o_markers.keys())
    entry = tuple(df.iloc[-2])
    assert (entry) == ('ABC_2', '', '', '', 'ABC-2, ABC2'), entry
    entry = tuple(df.iloc[-1])
    assert (entry) == ('CD3', 'eF450', 'goat', '', 'CD3a, CD3b'), entry

    masks = [df[x] == (y or '') for x, y in
             zip(cyc_markers.unique_keys, new_markers[-2])]
    masks = np.logical_and.reduce(masks)
    entry = list(df.loc[masks, 'aliases'])[-1]
    assert (entry) == 'CK14, CK_14, CK-14, CK14a', entry

def test_marker_comparator():
    m1 = Marker(name='CD4')
    m2 = Marker(name='CD4', fluor='ef570')
    m3 = Marker(name='CD4', fluor='ef570', anti='goat')
    m4 = Marker(name='CD4', anti='goat')
    m5 = Marker(name='CD4', duplicate='1')
    m6 = Marker(name='CD4', duplicate='2')
    m7 = Marker(name='CD45')

    assert Marker_Comparator(m1) != Marker_Comparator(m7)
    assert Marker_Comparator(m1) != Marker_Comparator(m2)
    assert Marker_Comparator(m1, fluor_sensitive=False) == \
        Marker_Comparator(m2, fluor_sensitive=False)
    assert Marker_Comparator(m1) == Marker_Comparator(m4)
    assert Marker_Comparator(m1, anti_sensitive=True) != \
        Marker_Comparator(m4, anti_sensitive=True)
    assert Marker_Comparator(m2) == Marker_Comparator(m3)
    assert Marker_Comparator(m2, anti_sensitive=True) != \
        Marker_Comparator(m3, anti_sensitive=True)

    assert Marker_Comparator(m1) == Marker_Comparator(m5)
    assert Marker_Comparator(m1) == Marker_Comparator(m6)

    assert repr(Marker_Comparator(m1)) == 'CD4', repr(Marker_Comparator(m1))
    assert repr(Marker_Comparator(m2)) == 'CD4_ef570', \
        repr(Marker_Comparator(m2))
    assert repr(Marker_Comparator(m2, fluor_sensitive=False)) == 'CD4', \
        repr(Marker_Comparator(m2, fluor_sensitive=False))
    assert repr(Marker_Comparator(m4)) == 'CD4', \
        repr(Marker_Comparator(m4))
    assert repr(Marker_Comparator(m4, anti_sensitive=True)) == 'CD4_goat', \
        repr(Marker_Comparator(m4, anti_sensitive=True))
    assert repr(Marker_Comparator(m5)) == 'CD4_1', repr(Marker_Comparator(m5))
    assert repr(Marker_Comparator(m6)) == 'CD4_2', repr(Marker_Comparator(m6))

    marker_set = set([Marker_Comparator(m1), Marker_Comparator(m5),
                      Marker_Comparator(m6)])
    assert len(marker_set) == 1

    assert Marker_Comparator(m1) in marker_set
    assert Marker_Comparator(m5) in marker_set
    assert Marker_Comparator(m6) in marker_set

max_stars_repo_name: ohsu-comp-bio/cycIF-DB

max_issues_repo_name: ohsu-comp-bio/cycIF-DB

max_forks_repo_name: ohsu-comp-bio/cycIF-DB