text:
# -*- coding: utf-8 -*-
# --------------------------------------
# @Time    : 2021/11/1 下午2:03
# @Author  : Jingyi Wang
# @Email   : Wangjingyi70@gmail.com
# @File    : detect.py
# @Software: PyCharm
# Description : None
# --------------------------------------
import json
import os
import cv2
import numpy as np

import os
import sys
import time
import yaml                                         # yaml配置模块
import torch
import logging                                      # 日志模块
import argparse
import os.path as osp
from tqdm import tqdm
from pathlib import Path                            # 路径操作模块
from torchvision import transforms
from torch.utils.tensorboard import SummaryWriter
from valid import validate
from dataset.dataloader import create_train_dataloader, create_val_dataloader
from utils.misc import AverageMeter, ProgressMeter
from loss.loss import KeyPointsMSELoss
from loss.optimizer import create_optimizer, create_lr_scheduler
from models.kpnet import build_model
from models.basenet import parameters_num
from utils.general import to_python_float, get_latest_run, check_file, colorstr, increment_dir
from utils.evaluate import calc_acc
from utils.torch_utils import select_device, save_checkpoint, get_saved_state
from utils.postprocess import *
from utils.visualize import *

mean_rgb = np.array([0.485, 0.456, 0.406], np.float32).reshape(1, 1, 3)
std_rgb = np.array([0.229, 0.224, 0.225], np.float32).reshape(1, 1, 3)
def normalize_image(image):
    return (image / 255. - mean_rgb) / std_rgb

def predict_keypoints_from_roi(model, img_roi):
    model.eval()
    with torch.no_grad():
        heatmaps = model(torch.FloatTensor(img_roi))
        preds, _ = get_max_preds(heatmaps.detach().cpu().numpy())
        preds = preds[0,: 4,:] *4 + 2
        return preds

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='Train keypoints network')
    parser.add_argument('--config', type=str, default='./config/hyp.yaml',
                        help='the path of the pretrained weights')
    parser.add_argument('--device', type=str, default='cpu', help='number of select device')
    opt = parser.parse_args()

    hyp = check_file(opt.config)
    assert len(hyp), '--hyp file must be specified'

    # 载入初始超参
    with open(hyp, encoding='UTF-8') as f:
        hyp = yaml.load(f, Loader=yaml.SafeLoader)

    device = select_device(opt.device, batch_size=hyp['batch_size'])
    model = build_model(hyp['pretrained'], hyp['num_keypoints'], is_train=False)
    model.load_state_dict(torch.load("./saved_weights/Epoch20_RGB.pth"))

    test_image = hyp["test_image"]
    test_label = hyp['test_label']
    image_size = np.array([1664,512])
    input_size = np.array(hyp["input_size"])
    heatmap_size = np.array(hyp["heatmap_size"])

    frame_cnt = 0
    for filename in os.listdir(test_label):
        if not filename.endswith(".json"):
            continue

        frame_cnt += 1

        image_path = os.path.join(test_image, '{}.png'.format(os.path.splitext(filename)[0]))
        image_bgr = cv2.imread(image_path)
        image = cv2.cvtColor(image_bgr, cv2.COLOR_BGR2RGB)

        with open(os.path.join(test_label, filename), 'r') as json_file:
            data = json.load(json_file)

            # Find all labels with shape "LShape" in one image
            for obj in data['objects']:
                if obj['shape'] != "LShape":
                    continue

                # Sequence: su -> sl -> ul -> lr
                keypoints = np.zeros((4, 2), dtype=float)

                keypoints[0] = obj['su']
                keypoints[1] = obj['sl']
                keypoints[2] = obj['ul']
                keypoints[3] = obj['lr']

                # Coordinate of labels less than zero
                if (keypoints < 0).sum() > 0:
                    continue
                # Coordinate of labels out of the image
                if (keypoints[:, 0] >= image_size[0]).sum() > 0:
                    continue
                if (keypoints[:, 1] >= image_size[1]).sum() > 0:
                    continue

                # Height of the car less than 50, too small
                if keypoints[3][1] - keypoints[2][1] < 50:
                    continue

                # Find the ROI region and add 16 pixels for each edge
                left_most = int(max(keypoints[:, 0].min() - 16, 0))
                right_most = int(min(keypoints[:, 0].max() + 16, image_size[0] - 1))
                width_roi = right_most - left_most

                top_most = int(max(keypoints[:, 1].min() - 16, 0))
                bottom_most = int(min(keypoints[:, 1].max() + 16, image_size[1] - 1))
                height_roi = bottom_most - top_most

                # Crop and resize the input image
                img_roi = image[top_most:bottom_most, left_most:right_most]
                img_roi = cv2.resize(img_roi, tuple(input_size))

                # Apply norm for the cropped image
                img_roi = normalize_image(img_roi)
                img_roi = img_roi.transpose(2, 0, 1).reshape([1, 3, input_size[1], input_size[0]])

                # Predict!
                keypoints = predict_keypoints_from_roi(model, img_roi)
                keypoints[:, 0] = keypoints[:, 0] / input_size[0] * width_roi + left_most
                keypoints[:, 1] = keypoints[:, 1] / input_size[1] * height_roi + top_most

                image_bgr = draw_3D_on_cv_image(image_bgr, keypoints)

        cv2.imshow("Current Frame", image_bgr)
        cv2.waitKey(50)
        # cv2.imwrite(os.path.join("Test", os.path.splitext(filename)[0] + ".png"), image_bgr)


text:
# coding=utf-8
# Copyright 2021 The jax_verify Authors.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# Lint as: python3
# pylint: disable=invalid-name
"""Crown bound propagation used in SDP verification."""

import functools
import jax.numpy as jnp
import jax_verify
from jax_verify.src.nonconvex import duals
from jax_verify.src.nonconvex import nonconvex
from jax_verify.src.nonconvex import optimizers
from jax_verify.src.nonconvex.optimizers import LinesearchFistaOptimizer as FistaOptimizer
from jax_verify.src.sdp_verify import utils
IntBound = utils.IntBound

def boundprop(params, x, epsilon, input_bounds, boundprop_type,
              **extra_boundprop_kwargs):
  """Computes interval bounds for NN intermediate activations.

  Args:
    params: Parameters for the NN.
    x: Batch of inputs to NN (dimension 2 for MLP or 4 for CNN)
    epsilon: l-inf perturbation to the input.
    input_bounds: Valid lower and upper for the NN as a tuple -- e.g. (0., 1.)
    boundprop_type: string, indicating method used for bound propagation, e.g.
      'crown_ibp' or 'nonconvex'
    **extra_boundprop_kwargs: any additional kwargs, passed directly to
      underlying boundprop method

  Returns:
    layer_bounds: upper and lower bounds across the layers of the NN as a list
    of IntBound-s.
  """
  boundprop_type_to_method = {
      'crown_ibp': _crown_ibp_boundprop,
      'nonconvex': _nonconvex_boundprop,
  }
  assert boundprop_type in boundprop_type_to_method, 'invalid boundprop_type'
  boundprop_method = boundprop_type_to_method[boundprop_type]
  return boundprop_method(params, x, epsilon, input_bounds,
                          **extra_boundprop_kwargs)

def _crown_ibp_boundprop(params, x, epsilon, input_bounds):
  """Runs CROWN-IBP for each layer separately."""
  def get_layer_act(layer_idx, inputs):
    act = utils.predict_cnn(params[:layer_idx], inputs)
    return act

  initial_bound = jax_verify.IntervalBound(
      jnp.maximum(x - epsilon, input_bounds[0]),
      jnp.minimum(x + epsilon, input_bounds[1]))

  out_bounds = [IntBound(
      lb_pre=None, ub_pre=None, lb=initial_bound.lower, ub=initial_bound.upper)]
  for i in range(1, len(params) + 1):
    fwd = functools.partial(get_layer_act, i)
    bound = jax_verify.crownibp_bound_propagation(fwd, initial_bound)
    out_bounds.append(
        IntBound(lb_pre=bound.lower, ub_pre=bound.upper,
                 lb=jnp.maximum(0, bound.lower),
                 ub=jnp.maximum(0, bound.upper)))
  return out_bounds

def _nonconvex_boundprop(params, x, epsilon, input_bounds,
                         nonconvex_boundprop_steps=100,
                         nonconvex_boundprop_nodes=128):
  """Wrapper for nonconvex bound propagation."""
  # Get initial bounds for boundprop
  init_bounds = utils.init_bound(x, epsilon, input_bounds=input_bounds,
                                 add_batch_dim=False)

  # Build fn to boundprop through
  all_act_fun = functools.partial(utils.predict_cnn, params,
                                  include_preactivations=True)

  # Collect the intermediate bounds.
  input_bound = jax_verify.IntervalBound(init_bounds.lb, init_bounds.ub)

  def _create_optimizer():
    return optimizers.OptimizingConcretizer(
        FistaOptimizer(num_steps=nonconvex_boundprop_steps),
        max_parallel_nodes=nonconvex_boundprop_nodes)
  nonconvex_bound_propagation = functools.partial(
      nonconvex.build_nonconvex_formulation,
      duals.WolfeNonConvexBound, _create_optimizer)
  _, intermediate_nonconvex_bounds = nonconvex_bound_propagation(
      all_act_fun, input_bound)

  bounds = [init_bounds]
  for nncvx_bound in intermediate_nonconvex_bounds:
    bounds.append(utils.IntBound(lb_pre=nncvx_bound.lower,
                                 ub_pre=nncvx_bound.upper,
                                 lb=jnp.maximum(nncvx_bound.lower, 0),
                                 ub=jnp.maximum(nncvx_bound.upper, 0)))
  return bounds

ext: py

reason: import jax,from jax


text:
import numpy as np
import cv2
import pickle
import requests
import serial
import math
import time
import json

face_cascade = cv2.CascadeClassifier('cascades/data/haarcascade_frontalface_alt2.xml')
recognizer = cv2.face.LBPHFaceRecognizer_create()
recognizer.read("trainer.yml")
labels = {"person_name": 1}
with open("labels.pickle", "rb") as f:
    og_labels = pickle.load(f)
    labels = {v:k for k, v in og_labels.items()}

arduino = serial.Serial('COM19', 115200)

horizontalOffset = 0
verticalOffset = 0

def move3DPoint():
    jsonInput = open('data.json', 'r')
    decoded = json.load(jsonInput)
    if decoded['pinch'] > 0.7:
        updatePosition(-decoded['position'][0], decoded['position'][1], decoded['position'][2])

def setX(angle):
        angle = truncate(angle, 5, 175)
        command = "<x," + str(angle) + ">"
        #print(command)
        arduino.write(command.encode())

def setY(angle):
    angle = truncate(angle, 5, 175)
    command = "<y," + str(angle) + ">"
    #print(command)
    arduino.write(command.encode())

def truncate(value, lowerBound, upperBound):
    value = max(lowerBound, value)
    value = min(value, upperBound)

    return value

def updatePosition(horizontal, vertical, depth):
    depth = min(depth, 1)

    xAngle = math.atan2(depth, horizontal+horizontalOffset) * 180 / math.pi
    yAngle = math.atan2(depth, (vertical+verticalOffset)) * 180 / math.pi

    setX(xAngle)
    setY(yAngle)

def face():
    updatePosition(0, -1, 1)
    setLed(20)

def setLed(pwm):
    command = "<z," + str(pwm) + ">"
    #print(command)
    arduino.write(command.encode())

def resetPosition():
    updatePosition(0, 0, 1)
    setLed(255)

def testLed():
    for pwm in range(64):
        setLed(4 * pwm)

def testServo(pause):
    values = [-1, 0 ,1]

cap = cv2.VideoCapture(1)
resetPosition()
#url = "http://100.64.72.17:8080/shot.jpg"#changes almost everytime you use it

while (True):
    #capturing data frame by frame
    ret, frame = cap.read()
    #frames_resp = requests.get(url)
    #frame_arr = np.array(bytearray(frames_resp.content), dtype=np.uint8)
    #frame = cv2.imdecode(frame_arr, -1)
    gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
    faces = face_cascade.detectMultiScale(gray, scaleFactor = 1.5, minNeighbors = 5)
    for (x, y, w, h) in faces:
        roi_gray = gray[y:y+h, x:x+w]
        roi_color = frame[y:y+h, x:x+w] #(ycord_start-height, ycord_end-height)

        id, conf = recognizer.predict(roi_gray)
        #print("confidence is", conf)

        if labels[id] == "parth-sareen" and conf > 05:
            #print(id)
            #print(labels[id])
            font = cv2.FONT_HERSHEY_SIMPLEX
            name = "parth is saved"
            color = (255,255,255)
            box_color = (0, 255, 0)
            stroke = 2

            cv2.putText(frame, name, (x,y), font, 1, color, stroke, cv2.LINE_AA)
            end_cord_x = x + w
            end_cord_y = y + h
            cv2.rectangle(frame, (x, y), (end_cord_x, end_cord_y), box_color, stroke)
            face()

        elif labels[id] == "parth-sareen":
            font = cv2.FONT_HERSHEY_SIMPLEX
            name = "Unsure"
            color = (255,255,255)
            box_color = (0, 0, 255)
            stroke = 2

            cv2.putText(frame, name, (x,y), font, 1, color, stroke, cv2.LINE_AA)
            end_cord_x = x + w
            end_cord_y = y + h
            cv2.rectangle(frame, (x, y), (end_cord_x, end_cord_y), box_color, stroke)
            setLed(255)
            move3DPoint()

        img_item = "7.png" #creates image for testing

        cv2.imwrite(img_item, roi_color)
        #print(x,y,w,h)

        #color = (0, 0, 255) #BGR
        #stroke = 2
        #end_cord_x = x + w
        #end_cord_y = y + h
        #cv2.rectangle(frame, (x, y),(end_cord_x, end_cord_y), color, stroke)

    #display frame
    cv2.imshow('frame',frame)
    if cv2.waitKey(20) & 0xFF == ord('q'):
        break

#when proccesses are done executing, release VideoCapture
cap.release()
cv2.destroyAllWindows()

max_stars_repo_name: ParthSareen/SmartLamp

max_issues_repo_name: ParthSareen/SmartLamp

max_forks_repo_name: ParthSareen/SmartLamp

reason: import numpy


text:
import os
import sys
import time
import torch
import torch.nn.functional as F
import torch.nn as nn
from torch.optim import lr_scheduler
from torch.autograd import Variable
import numpy as np
import collections
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import random

class Optimizer():
    def __init__(self, ndownsamples=4, currn=5, avgn=20, lambda1=0.4, lr=0.2, eps=0.01, min_iter=20, max_iter=1000):
        self.ndownsamples = ndownsamples
        self.currn = currn
        self.avgn = avgn
        self.lambda1 = lambda1
        self.lr = lr
        self.eps = eps
        self.identities = {}
        self.min_iter = min_iter
        self.max_iter = max_iter

    @staticmethod
    def center(var, dims, d):
        if not isinstance(d, collections.Sequence):
            d = [d for i in range(len(dims))]
        for idx, dim in enumerate(dims):
            if d[idx] == 0:
                continue
            var = var.narrow(dim, int(d[idx]/2), int(var.size()[dim] - d[idx]))
        return var

    def get_identity_grid(self, dim, cache=True):
        if dim not in self.identities:
            gx, gy = np.linspace(-1, 1, dim), np.linspace(-1, 1, dim)
            I = np.stack(np.meshgrid(gx, gy))
            I = np.expand_dims(I, 0)
            I = torch.FloatTensor(I)
            I = torch.autograd.Variable(I, requires_grad=False)
            I = I.permute(0,2,3,1).cuda()
            self.identities[dim] = I
        if cache:
            return self.identities[dim]
        else:
            return self.identities[dim].clone()

    def jacob(self, fields):
        def dx(f):
            p = Variable(torch.zeros((1,1,f.size(1),2))).cuda()
            return torch.cat((p, f[:,2:,:,:] - f[:,:-2,:,:], p), 1)
        def dy(f):
            p = Variable(torch.zeros((1,f.size(1),1,2))).cuda()
            return torch.cat((p, f[:,:,2:,:] - f[:,:,:-2,:], p), 2)
        fields = sum(map(lambda f: [dx(f), dy(f)], fields), [])
        field = torch.sum(torch.cat(fields, -1) ** 2, -1)
        return field

    def penalty(self, fields, mask=1):
        jacob = self.jacob(fields)
        jacob = torch.mul(jacob, mask)
        return torch.sum(jacob)

    def render(self, src, field):
        src, field = torch.FloatTensor(src).cuda(), torch.FloatTensor(field).cuda()
        src, field = Variable(src).unsqueeze(0).unsqueeze(0), Variable(field).unsqueeze(0)
        #print(src, field)
        y =  F.grid_sample(src, field + self.get_identity_grid(field.size(2)), mode='bilinear')
        return  y.data.cpu().numpy()

    def process(self, s, t, crop=0, mask=1):
        print(s.shape, t.shape)
        downsample = lambda x: nn.AvgPool2d(2**x,2**x, count_include_pad=False) if x > 0 else (lambda y: y)
        upsample = nn.Upsample(scale_factor=2, mode='bilinear')
        s, t = torch.FloatTensor(s), torch.FloatTensor(t)
        src = Variable((s - torch.mean(s)) / torch.std(s)).cuda().unsqueeze(0).unsqueeze(0)
        target = Variable((t - torch.mean(t)) / torch.std(t)).cuda().unsqueeze(0).unsqueeze(0)
        mask = Variable(torch.FloatTensor(mask)).cuda().unsqueeze(0)
        dim = int(src.size()[-1] / (2 ** (self.ndownsamples - 1)))
        field = Variable(torch.zeros((1,dim,dim,2))).cuda().detach()
        field.requires_grad = True
        updates = 0
        masking = not list(mask.shape)[-1] == 1
        for downsamples in reversed(range(self.ndownsamples)):
            src_, target_ = downsample(downsamples)(src).detach(), downsample(downsamples)(target).detach()
            mask_ = downsample(downsamples)(mask).detach() if masking else mask.detach()
            mask_.requires_grad = False
            src_.requires_grad = False
            target_.requires_grad = False
            field = field.detach()
            field.requires_grad = True
            opt = torch.optim.SGD([field], lr=self.lr/(downsamples+1))
            #sched = lr_scheduler.StepLR(opt, step_size=1, gamma=0.995)
            costs = []
            start_updates = updates
            print(downsamples)
            while True:
                updates += 1
                pred = F.grid_sample(src_, field + self.get_identity_grid(field.size(2)), mode='bilinear')
                if masking:
                    penalty1 = self.penalty([self.center(field, (1,2), 128 / (2**downsamples))], self.center(mask_, (1,2), 128 / (2**downsamples)))
                else:
                    penalty1 = self.penalty([self.center(field, (1,2), 128 / (2**downsamples))])
                diff = torch.mean(self.center((pred - target_)**2, (-1,-2), 128 / (2**downsamples)))
                cost = diff + penalty1 * self.lambda1/(downsamples+1)
                print(cost.data.cpu().numpy())
                costs.append(cost)
                cost.backward()
                opt.step()
                #sched.step()uniform
                opt.zero_grad()
                if len(costs) > self.avgn + self.currn and len(costs)>self.min_iter:
                    hist = sum(costs[-(self.avgn+self.currn):-self.currn]).data[0] / self.avgn
                    curr = sum(costs[-self.currn:]).data[0] / self.currn
                    if abs((hist-curr)/hist) < self.eps/(2**downsamples) or len(costs)>self.max_iter:
                        break
            #print downsamples, updates - start_updates
            if downsamples > 0:
                field = upsample(field.permute(0,3,1,2)).permute(0,2,3,1)
        #print(cost.data[0], diff.data[0], penalty1.data[0])
        print('done:', updates)
        print(field.shape)
        return self.center(field, (1,2), crop*2).data.cpu().numpy()[0]

if __name__ == '__main__':
    o = Optimizer()
    print('Testing...')
    s = np.random.uniform(0, 1, (256,256)).astype(np.float32)
    t = np.random.uniform(0, 1, (256,256)).astype(np.float32)

    flow = o.process(s, t)
    print(flow.shape)
    assert flow.shape == (1,256,256,2)

    flow = o.process(s, t, crop=10)
    assert flow.shape == (1,236,236,2)

    print ('All tests passed.')

max_stars_repo_name: emitch/SEAMLeSS

max_issues_repo_name: emitch/SEAMLeSS

max_forks_repo_name: emitch/SEAMLeSS


text:
import pytest
import jax
from jax import scipy as jscipy
from pyscfad.lib import numpy as np
from pyscfad.lib import linalg_helper as linalg

def test_eigh():
    a = np.ones((2,2))
    b = np.eye(2)

    w0, v0 = jscipy.linalg.eigh(a)
    w, v = linalg.eigh(a)

    assert abs(w-w0).max() < 1e-10
    assert abs(v-v0).max() < 1e-10

    jac0 = jax.jacfwd(jscipy.linalg.eigh)(a)
    jac = jax.jacfwd(linalg.eigh)(a)

    assert abs(jac[0] - jac0[0]).max() < 1e-10
    assert abs(jac[1] - jac0[1]).max() < 1e-10

    #finite difference
    m = -0.0005
    p = 0.0005
    b_m = np.array([[1., 0.],[0., 1.+m]])
    w_m, v_m = linalg.eigh(a, b_m)
    b_p = np.array([[1, 0.],[0., 1.+p]])
    w_p, v_p = linalg.eigh(a, b_p)

    g0 = (-v_p - v_m) / 0.001 # -v_p due to specific gauge

    jac = jax.jacfwd(linalg.eigh, argnums=1)(a, b)
    g = jac[1][:,:,1,1]
    assert abs(g-g0).max() < 1e-7

max_stars_repo_name: yangdatou/pyscfad

max_issues_repo_name: yangdatou/pyscfad

max_forks_repo_name: yangdatou/pyscfad

reason: import numpy,import scipy,import jax,from jax


text:
"""
spincnstr2.py - constrained robust optimization of the spin system
"""

from argparse import ArgumentParser
from copy import copy
import os

from autograd.extend import Box
import autograd.numpy as anp
import jax.numpy as jnp
from filelock import FileLock, Timeout
import h5py
import numpy as np
from qoc import (
    grape_schroedinger_discrete,
    evolve_schroedinger_discrete,
    lqr,
)
from qoc.core.common import(
    clip_control_norms,
)
from qoc.standard import (
    TargetStateInfidelity,
    conjugate_transpose,
    gram_schmidt,
    matmuls,
    project,
    get_annihilation_operator,
    get_creation_operator,
    SIGMA_Z, SIGMA_X,
    generate_save_file_path,
    SGD, Adam,
)
from qoc.models import Cost
from qoc.core.mathmethods import interpolate_linear_set, integrate_rkdp5
from qutip import essolve, mesolve, Qobj

# Define paths.
EXPERIMENT_META = "spin"
EXPERIMENT_NAME = "spincnstr2"
WDIR_PATH = os.environ["ROBUST_QOC_PATH"]
SAVE_PATH = os.path.join(WDIR_PATH, "out", EXPERIMENT_META, EXPERIMENT_NAME)
SAVE_FILE = EXPERIMENT_NAME

# Define compute resources.
CORE_COUNT = 8
os.environ["MKL_NUM_THREADS"] = "{}".format(CORE_COUNT)
os.environ["OPENBLAS_NUM_THREADS"] = "{}".format(CORE_COUNT)

# Define experimental constants.
MAX_AMP_0 = 2 * np.pi * 3e-1
OMEGA = 2 * np.pi * 1e-2

# Define the system.
HILBERT_SIZE = 2
H_SYSTEM_0 = SIGMA_Z / 2
H_CONTROL_0 = SIGMA_X / 2
hamiltonian = lambda controls, hargs, time: (
    hargs[0] * H_SYSTEM_0
    + controls[0] * H_CONTROL_0
)
HAMILTONIAN_ARGS = np.array([OMEGA])
MAX_CONTROL_NORMS = np.array([MAX_AMP_0])

# Define the optimization.
OPTIMIZER = Adam()
COMPLEX_CONTROLS = False
CONTROL_COUNT = 1
EVOLUTION_TIME = 120
CONTROL_EVAL_COUNT = 2 * int(EVOLUTION_TIME) + 1
CONTROL_EVAL_TIMES = np.linspace(0, EVOLUTION_TIME, CONTROL_EVAL_COUNT)
ITERATION_COUNT = 100

# Define the problem.
INITIAL_STATE_0 = np.array([[1], [0]], dtype=np.float64)
TARGET_STATE_0 = np.array([[0], [1]], dtype=np.float64)
INITIAL_STATES = np.stack((INITIAL_STATE_0,),)
TARGET_STATES = np.stack((TARGET_STATE_0,),)
TARGET_STATES_DAGGER = conjugate_transpose(TARGET_STATES)
FIDELITY_CONSTRAINT = 1e-3
FIDELITY_MULTIPLIER = 5
FIDELITY_MULTIPLIER_SCALE = 5
INITIAL_FIDELITY_ROBUSTNESS = 0

INITIAL_ASTATE = np.hstack([
    INITIAL_STATES.ravel(),
],)
states_count = INITIAL_STATES.shape[0]
states_offset = 0
states_shape = INITIAL_STATES.shape
states_size = np.prod(states_shape)
get_states = lambda astate: anp.reshape(astate[states_offset:states_offset + states_size], states_shape)

# Define the problem.
class Fidelity(Cost):
    """
    A cost function to promote fidelity.
    """
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

    def cost(self, controls, final_astate):
        final_states = get_states(final_astate)
        inner_products = anp.matmul(TARGET_STATES_DAGGER, final_states)[:, 0, 0]
        fidelities = anp.real(inner_products * anp.conjugate(inner_products))
        fidelity_normalized = anp.sum(fidelities) / states_count
        cost_ = 1 - fidelity_normalized

        augmented_cost = self.augment_cost(cost_)
        return augmented_cost

class FidelityRobustness(Cost):
    """
    A cost function to promote robustness of fidelity.
    """
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

    def cost(self, controls, final_astate):
        cost_ = anp.abs(final_astate[3][0][1][0])
        augmented_cost = self.augment_cost(cost_)
        return augmented_cost

COSTS = [
    Fidelity(),
    # FidelityRobustness(),
]

def rhs(astate, controls, time):
    delta_astate = anp.zeros_like(astate)
    # controls_ = interpolate_linear_set(time, CONTROL_EVAL_TIMES, controls)
    controls_ = controls[0]
    hamiltonian_ = hamiltonian(controls_, HAMILTONIAN_ARGS, time)
    states = get_states(astate)
    delta_states = -1j * anp.matmul(hamiltonian_, states)
    # # dpsi_dw
    # delta_dpsi_dw = -1j * anp.matmul(H_SYSTEM_0, astate[0]) -1j * anp.matmul(hamiltonian_, astate[1])
    # # d2psi_dw2
    # delta_d2psi_d21 = -2j * anp.matmul(H_SYSTEM_0, astate[1]) -1j * anp.matmul(hamiltonian_, astate[2])
    # # d2f_dw2
    # delta_d2f_dw2 = anp.sum(
    #     -1j * anp.matmul(TARGET_STATES_DAGGER,
    #                      (2 * anp.matmul(H_SYSTEM_0, astate[1])
    #                       + anp.matmul(hamiltonian_, astate[2]))
    #     )[:, 0, 0])

    delta_astate = anp.hstack([
        delta_states.ravel(),
    ])
    return delta_astate

# Define controls update procedure.
# Generate orthonormal basis for zero net flux constraint manifold.
ZF_BASIS = np.zeros((CONTROL_EVAL_COUNT - 1, CONTROL_EVAL_COUNT))
for i in range(CONTROL_EVAL_COUNT - 1):
    ZF_BASIS[i][0] = 1
    ZF_BASIS[i][i + 1] = -1
ZF_BASIS = gram_schmidt(ZF_BASIS)

def impose_control_conditions(controls):
    # Impose zero at boundaries.
    # controls[0, :] = 0
    # controls[-1, :] = 0

    # Project onto zero net flux constraint manifold.
    # controls_ = anp.zeros_like(controls)
    # for control_index in range(controls.shape[1]):
    #     controls_[:, control_index] = project(controls[:, control_index], ZF_BASIS)
    # controls = controls_

    # Impose control norm constraints.
    # clip_control_norms(controls, MAX_CONTROL_NORMS)

    return controls

# Choose initial controls.
GRAB_CONTROLS = False
GEN_CONTROLS = True
if GRAB_CONTROLS:
    controls_file_path = os.path.join(SAVE_PATH, "00005_{}.h5".format(EXPERIMENT_NAME))
    controls_lock_file_path = "{}.lock".format(controls_file_path)
    try:
        with FileLock(controls_lock_file_path):
            with h5py.File(controls_file_path) as f:
                index = np.argmin(f["error"][()])
                controls_ = f["controls"][index][()]
    except Timeout:
        print("Timeout on {}."
              "".format(controls_lock_file_path))
elif GEN_CONTROLS:
    controls_ = np.ones((CONTROL_EVAL_COUNT, CONTROL_COUNT))
    # mid_index = int(np.floor(CONTROL_EVAL_COUNT / 2))
    # controls_[:mid_index] = -controls_[:mid_index]
    # controls_ = controls_ * (3 * MAX_CONTROL_NORMS / 4)
else:
    controls_ = None
INITIAL_CONTROLS = controls_

# Define output.
LOG_ITERATION_STEP = 1
SAVE_ITERATION_STEP = 1

LQR_CONFIG = {
    "complex_controls": COMPLEX_CONTROLS,
    "costs": COSTS,
    "evolution_time": EVOLUTION_TIME,
    "impose_control_conditions": impose_control_conditions,
    "initial_astate": INITIAL_ASTATE,
    "initial_controls": INITIAL_CONTROLS,
    "iteration_count": ITERATION_COUNT,
    "rhs": rhs,
}

def do_lqr():
    config = copy(LQR_CONFIG)
    result = lqr(**config)
    final_states = get_states(result.final_astate)
    while isinstance(final_states, Box):
        final_states = final_states._value
    c = np.matmul(TARGET_STATES_DAGGER, final_states)[0, 0, 0]
    print("final_states:\n{}\nc:\n{}"
          "".format(final_states, c))

def do_qutip():
    h0_qutip = Qobj(OMEGA * H_SYSTEM_0)
    c0_qutip = np.ones_like(INITIAL_CONTROLS[:, 0])
    h1_qutip = Qobj(H_CONTROL_0)
    c1_qutip = INITIAL_CONTROLS[:, 0]
    h_list = [
        [h0_qutip, c0_qutip],
        [h1_qutip, c1_qutip],
    ]
    initial_state_qutip = Qobj(INITIAL_STATE_0)
    tlist = np.linspace(0, EVOLUTION_TIME, CONTROL_EVAL_COUNT)
    c_op_list = []
    e_ops = []
    result_qutip = mesolve(h_list, initial_state_qutip, tlist, c_op_list, e_ops)
    final_state_qutip = result_qutip.states[-1].full()
    ip = np.matmul(conjugate_transpose(final_state_qutip), final_state_qutip)
    print("final_state_qutip:\n{}\nip:\n{}"
          "".format(final_state_qutip, ip))

def do_evolve():
    controls = INITIAL_CONTROLS
    pc = np.polyfit(CONTROL_EVAL_TIMES, controls, 10)
    initial_state = INITIAL_STATE_0
    hargs = HAMILTONIAN_ARGS
    initial_time = 0
    final_times = np.array([EVOLUTION_TIME])

def rhs_(time, state):
    # controls_ = interpolate_linear_set(time, CONTROL_EVAL_TIMES, controls)
    hamiltonian_ = hamiltonian(controls_, hargs, time)
    return -1j * np.matmul(hamiltonian_, state)
    # return -1j * np.matmul(OMEGA * H_SYSTEM_0, state)
    # result = integrate_rkdp5(rhs_, final_times, initial_time, initial_state)
    # final_state = result
    # ip = np.matmul(conjugate_transpose(final_state), final_state)[0, 0]
    # print("final_state:\n{}\nip:\n{}"
    #       "".format(final_state, ip))

def main():
    parser = ArgumentParser()
    parser.add_argument("--evolve", dest="evolve", action="store_true", default=False)
    parser.add_argument("--lqr", dest="lqr", action="store_true", default=False)
    parser.add_argument("--qutip", dest="qutip", action="store_true", default=False)
    args = vars(parser.parse_args())
    run_evolve = args["evolve"]
    run_lqr = args["lqr"]
    run_qutip = args["qutip"]

    if run_evolve:
        do_evolve()
    if run_lqr:
        do_lqr()
    if run_qutip:
        do_qutip()

if __name__ == "__main__":
    main()

max_stars_repo_name: SchusterLab/rbqoc

max_issues_repo_name: SchusterLab/rbqoc

max_forks_repo_name: SchusterLab/rbqoc

reason: import numpy,import jax


text:
#  Created byMartin.cz
#  Copyright (c) Martin Strohalm. All rights reserved.

import numpy
import perrot

# prepare data
count = 25
data = []

for x in numpy.linspace(-numpy.pi, numpy.pi, count):
    for y in numpy.linspace(-numpy.pi, numpy.pi, count):
        z = (numpy.sin(x) + numpy.sin(y))/2
        data.append((x/numpy.pi, y/numpy.pi, z))

data = numpy.array(data)
z_min = numpy.min(data[:, 2])
z_max = numpy.max(data[:, 2])

# init plot
plot = perrot.Plot(
    x_axis_title = "pi",
    x_axis_ticker_major_step = 0.5,
    y_axis_title = "pi",
    y_axis_ticker_major_step = 0.5,)

# init color bar
color_bar = perrot.ColorBar(
    position = perrot.RIGHT,
    gradient = perrot.colors.YlOrBr,
    margin = (0, 0, 0, 15))

# init z-axis
z_axis = perrot.LinAxis(
    title = "sin(x) + sin(y)",
    position = perrot.RIGHT,
    scale_in_range = (z_min, z_max),
    level = 3,
    margin = 0,
    show_line = False)

# add bar and axis
plot.add(color_bar)
plot.add(z_axis)

# map bar scale to z-axis
plot.map(color_bar, 'scale', z_axis, 'scale')

# add series
series = plot.rects(
    data = data,
    x = lambda d: d[0],
    y = lambda d: d[1],
    width = 2/(count-1),
    height = 2/(count-1),
    margin = 0,
    line_width = 0,
    fill_color = lambda d: color_bar.get_color(d[2]))

# show plot
plot.zoom()
plot.view("Heatmap Plot", 660, 500)

max_stars_repo_name: xxao/perrot

max_issues_repo_name: xxao/perrot

max_forks_repo_name: xxao/perrot


text:
import numpy as np
import glob
import cv2
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import pickle
import os

# Visualization function
def line_visualization(binary_warped, nonzerox, nonzeroy, left_lane_inds, right_lane_inds, left_fitx, right_fitx, ploty, margin):
    # Create an image to draw on and an image to show the selection window
    out_img = np.dstack((binary_warped, binary_warped, binary_warped))*255
    window_img = np.zeros_like(out_img)
    # Color in left and right line pixels
    out_img[nonzeroy[left_lane_inds], nonzerox[left_lane_inds]] = [255, 0, 0]
    out_img[nonzeroy[right_lane_inds], nonzerox[right_lane_inds]] = [0, 0, 255]

    # Generate a polygon to illustrate the search window area (margins)
    # And recast the x and y points into usable format for cv2.fillPoly()
    left_line_window1 = np.array([np.transpose(np.vstack([left_fitx-margin, ploty]))])
    left_line_window2 = np.array([np.flipud(np.transpose(np.vstack([left_fitx+margin,
                              ploty])))])
    left_line_pts = np.hstack((left_line_window1, left_line_window2))
    right_line_window1 = np.array([np.transpose(np.vstack([right_fitx-margin, ploty]))])
    right_line_window2 = np.array([np.flipud(np.transpose(np.vstack([right_fitx+margin,
                              ploty])))])
    right_line_pts = np.hstack((right_line_window1, right_line_window2))

    # Draw the lane onto the warped blank image
    cv2.fillPoly(window_img, np.int_([left_line_pts]), (0,255, 0))
    cv2.fillPoly(window_img, np.int_([right_line_pts]), (0,255, 0))
    result = cv2.addWeighted(out_img, 1, window_img, 0.3, 0)

    # Plot the polynomial lines onto the image
    plt.plot(left_fitx, ploty, color='yellow')
    plt.plot(right_fitx, ploty, color='yellow')

    # View your output
    plt.imshow(result)

# Curvature computation
def compute_curvature(img, allx, ally):
    ym_per_pix = 20 / 720
    xm_per_pix = 3.7 / 800

    fit_cr = np.polyfit(ally * ym_per_pix, allx * xm_per_pix, 2)
    curvature = ((1 + (2*fit_cr[0]*img.shape[0]*ym_per_pix + fit_cr[1])**2)**1.5) / np.absolute(2*fit_cr[0])

    return curvature

# Function to fit polynomial, used by subsequent functions
def fit_poly(img_shape, allx, ally):
    ### Fit a second order polynomial to each with np.polyfit() ###
    fit = np.polyfit(ally, allx, 2)
    # Generate x and y values for plotting
    ploty = np.linspace(0, img_shape[0]-1, img_shape[0])
    ### Calc polynomial using ploty and fit ###
    fitx = fit[0]*ploty**2 + fit[1]*ploty + fit[2]

    return fit, fitx, ploty

# Function to draw the lane
def drawLane(undist, binary_warped, ploty, left_fitx, right_fitx, curvature, distance, Minv):
    # Create an image to draw the lines on
    warp_zero = np.zeros_like(binary_warped).astype(np.uint8)
    color_warp = np.dstack((warp_zero, warp_zero, warp_zero))

    # Recast the x and y points into usable format for cv2.fillPoly()
    pts_left = np.array([np.transpose(np.vstack([left_fitx, ploty]))])
    pts_right = np.array([np.flipud(np.transpose(np.vstack([right_fitx, ploty])))])
    pts = np.hstack((pts_left, pts_right))

    # Draw the lane onto the warped blank image
    cv2.fillPoly(color_warp, np.int_([pts]), (0,255, 0))

    # Warp the blank back to original image space using inverse perspective matrix (Minv)
    newwarp = cv2.warpPerspective(color_warp, Minv, (color_warp.shape[1], color_warp.shape[0]))
    # Combine the result with the original image
    result = cv2.addWeighted(undist, 1, newwarp, 0.3, 0)

    result = cv2.putText(result, "Radius of the curvature: %.0fm" % curvature, (50, 50), cv2.FONT_HERSHEY_SIMPLEX,
                   1, (255, 0, 0), 2, cv2.LINE_AA)
    if(distance < 0):
        result = cv2.putText(result, "Distance: %.2fm on the left of center" % np.absolute(distance), (50, 100),
                             cv2.FONT_HERSHEY_SIMPLEX, 1, (255, 0, 0), 2, cv2.LINE_AA)
    elif(distance > 0):
        result = cv2.putText(result, "Distance: %.2fm on the right of center" % np.absolute(distance), (50, 100),
                             cv2.FONT_HERSHEY_SIMPLEX, 1, (255, 0, 0), 2, cv2.LINE_AA)
    else:
        result = cv2.putText(result, "Distance: %.2fm" % np.absolute(distance), (50, 100),
                             cv2.FONT_HERSHEY_SIMPLEX, 1, (255, 0, 0), 2, cv2.LINE_AA)

    return result


text:
from __future__ import absolute_import
from __future__ import print_function
from __future__ import division

import sys, os, glob
import json
import argparse
import numpy as np
import open3d as o3d
import pickle

import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.nn import init
import torch.optim as optim
from torch.optim import lr_scheduler

import smplx
from human_body_prior.tools.model_loader import load_vposer
import chamfer_pytorch.dist_chamfer as ext

from cvae import BodyParamParser, HumanCVAES1, GeometryTransformer

class TestOP:
    def __init__(self, testconfig):
        for key, val in testconfig.items():
            setattr(self, key, val)

        if not os.path.exists(self.ckpt_dir):
            print('--[ERROR] checkpoints do not exist')
            sys.exit()

        #define model

        if self.use_cont_rot:
            n_dim_body=72+3
        else:
            n_dim_body=72

        self.model_h_latentD = 256
        self.model_h = HumanCVAES1(latentD=self.model_h_latentD,
                                  n_dim_body=n_dim_body)

        ### body mesh model
        self.vposer, _ = load_vposer(self.vposer_ckpt_path, vp_model='snapshot')
        self.body_mesh_model = smplx.create(self.human_model_path,
                                            model_type='smplx',
                                            gender='neutral', ext='npz',
                                            num_pca_comps=12,
                                            create_global_orient=True,
                                            create_body_pose=True,
                                            create_betas=True,
                                            create_left_hand_pose=True,
                                            create_right_hand_pose=True,
                                            create_expression=True,
                                            create_jaw_pose=True,
                                            create_leye_pose=True,
                                            create_reye_pose=True,
                                            create_transl=True,
                                            batch_size=self.n_samples
                                            )

    def data_preprocessing(self, img, modality, target_domain_size=[128, 128],
                                filename=None):

        """
        input:
            - img (depthmap or semantic map): [height, width].
            - modality: 'depth' or 'seg'
        output:
            canvas: with shape of target_domain_size, where the input is in the
                    center tightly, with shape target_domain_size
            factor: the resizing factor
        """

        # prepare the canvas
        img_shape_o = img.shape
        canvas = torch.zeros([1,1]+target_domain_size, dtype=torch.float32,
                             device=self.device)

        # filter out unavailable values
        if modality == 'depth':
            img[img>6.0]=6.0

        if modality == 'seg':
            img[img>41] = 41

        max_val = torch.max(img)
        _img = 2* img / max_val - 1.0

        # put _img to the canvas
        if img_shape_o[0]>= img_shape_o[1]:
            factor = float(target_domain_size[0]) / img_shape_o[0]
            target_height = target_domain_size[0]
            target_width = int(img_shape_o[1] * factor) //2 *2

            # for depth map we use bilinear interpolation in resizing
            # for segmentation map we use bilinear interpolation as well.
            # note that float semantic label is not real in practice, but
            # helpful in our work
            target_size = [target_height, target_width]

            _img = _img.view(1,1,img_shape_o[0],img_shape_o[1])
            img_resize = F.interpolate(_img, size=target_size, mode='bilinear',
                                        align_corners=False)

            na = target_width
            nb = target_domain_size[1]
            lower = (nb //2) - (na //2)
            upper = (nb //2) + (na //2)

            canvas[:,:,:, lower:upper] = img_resize

        else:
            factor = float(target_domain_size[1]) / img_shape_o[1]

            target_height = int(factor*img_shape_o[0]) //2 *2
            target_width = target_domain_size[1]

            target_size = [target_height, target_width]
            _img = _img.view(1,1,img_shape_o[0],img_shape_o[1])
            img_resize = F.interpolate(_img, size=target_size, mode='bilinear',
                                        align_corners=False)

            na = target_height
            nb = target_domain_size[0]
            lower = (nb //2) - (na //2)
            upper = (nb //2) + (na //2)

            canvas[:,:,lower:upper, :] = img_resize

        return canvas, factor, max_val

    def test(self):

        self.model_h.eval()
        self.model_h.to(self.device)
        self.vposer.to(self.device)
        self.body_mesh_model.to(self.device)

        ## load checkpoints
        ckp_list = sorted(glob.glob(os.path.join(self.ckpt_dir,'epoch-*.ckp')),
                            key=os.path.getmtime)
        ckp_path = ckp_list[-1]
        checkpoint = torch.load(ckp_path)
        print('[INFO] load checkpoints: ' + ckp_path)
        self.model_h.load_state_dict(checkpoint['model_h_state_dict'])

        ## read data from here!
        cam_file_list = sorted(glob.glob(self.test_data_path+'/cam_*'))

        for ii, cam_file in enumerate(cam_file_list):

            cam_params = np.load(cam_file, allow_pickle=True, encoding='latin1').item()

            depth0 = torch.tensor(np.load(cam_file.replace('cam','depth')), dtype=torch.float32, device=self.device)
            seg0 = torch.tensor(np.load(cam_file.replace('cam','seg')),dtype=torch.float32, device=self.device)

            cam_ext = torch.tensor(cam_params['cam_ext'], dtype=torch.float32, device=self.device).unsqueeze(0) #[1,4,4]
            cam_int = torch.tensor(cam_params['cam_int'], dtype=torch.float32, device=self.device).unsqueeze(0) # [1,3,3]

            depth, _, max_d = self.data_preprocessing(depth0, 'depth', target_domain_size=[128, 128]) #[1,1,128,128]
            max_d = max_d.view(1)
            seg, _, _ = self.data_preprocessing(seg0, 'depth', target_domain_size=[128, 128])

            xs = torch.cat([depth, seg], dim=1)
            xs_batch = xs.repeat(self.n_samples, 1,1,1)
            max_d_batch = max_d.repeat(self.n_samples)
            cam_int_batch = cam_int.repeat(self.n_samples, 1,1)
            cam_ext_batch = cam_ext.repeat(self.n_samples, 1,1)

            xhnr_gen = self.model_h.sample(xs_batch)
            xhn_gen = GeometryTransformer.convert_to_3D_rot(xhnr_gen)
            xh_gen = GeometryTransformer.recover_global_T(xhn_gen, cam_int_batch, max_d_batch)

            body_param_list = BodyParamParser.body_params_encapsulate(xh_gen)

            if not os.path.exists(self.outdir):
                os.makedirs(self.outdir)

            print('[INFO] save results to: '+self.outdir)
            for jj, body_param in enumerate(body_param_list):
                body_param['cam_ext'] = cam_ext_batch.detach().cpu().numpy()
                body_param['cam_int'] = cam_int_batch.detach().cpu().numpy()
                outfilename = os.path.join(self.outdir, 'body_gen_{:06d}.pkl'.format(self.n_samples*ii+jj))
                outfile = open(outfilename, 'wb')
                pickle.dump(body_param, outfile)
                outfile.close()

if __name__ == '__main__':

    scene_list = ['17DRP5sb8fy-bedroom', '17DRP5sb8fy-familyroomlounge',
                  '17DRP5sb8fy-livingroom', 'sKLMLpTHeUy-familyname_0_1',
                  'X7HyMhZNoso-livingroom_0_16', 'zsNo4HB9uLZ-bedroom0_0',
                  'zsNo4HB9uLZ-livingroom0_13']

    mp3dr_path = '/is/cluster/yzhang/mp3d-rooms'

    for fileid in range(len(scene_list)):
        print('[INFO] processing: '+scene_list[fileid])
        testconfig={
            'outdir': 'results_habitat_stage1_sceneloss/virtualcams/'+scene_list[fileid],
            'ckpt_dir': 'checkpoints_prox_modelS1_batch32_epoch30_LRS0.0003_LRH0.0003_LossVposer0.001_LossKL0.1_LossContact0.001_LossCollision0.01',
            'human_model_path': '/is/ps2/yzhang/body_models/VPoser',
            'vposer_ckpt_path': '/is/ps2/yzhang/body_models/VPoser/vposer_v1_0',
            'device': torch.device("cuda" if torch.cuda.is_available() else "cpu"),
            'test_data_path': os.path.join(mp3dr_path ,scene_list[fileid]+'-sensor'),
            'n_samples': 200,
            'use_cont_rot':True
        }

        test_op = TestOP(testconfig)

        test_op.test()


text:
# A model for visualizing the light output from a star + planet or an eclipsing binary system of two stars.

from argparse import ZERO_OR_MORE
from re import S
import sys
import math
from tkinter.constants import ANCHOR
import numpy as np
import tkinter
from tkinter import CENTER, NW, StringVar, messagebox
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use("TkAgg")
from matplotlib.figure import Figure
from matplotlib.backends.backend_tkagg import (FigureCanvasTkAgg)
from matplotlib.ticker import (MultipleLocator)
from systemModel import SystemModel
from os.path import exists

# planet radius RJupiter to RSun
RJupiterToRSun = 0.10047

# create the system model
pm = SystemModel(
    1.0, # star radius in units of Solar radius
    2.0, # planet orbital radius in units of R_Star
    1.0 * RJupiterToRSun, # planet radius in units of Jupiter radius
    0.8, # planet luminosity(0 = planet, > 1 eclipsing binary star)
    15.0, # inclination view angle (degrees)
    0.5 # star limb-darkening factor (0 = none)
    )

class App_Window(tkinter.Tk):
    def __init__(self, parent):
        tkinter.Tk.__init__(self,parent)
        self.initialize()
    def appQuit(self):
        sys.exit()
    def loadInputs(self):
        # Load the model inputs from a local text file.
        # First ensure that the modelInputs file exists
        if exists("modelInputs.txt"):
            inFile = open("modelInputs.txt","r")
            for line in inFile:
                # ignore comment lines
                if line.index('#') == 0:
                    continue
                split1 = line.split('#')
                split2 = split1[0].split(':')
                valFloat = float(split2[1])
                if line.startswith('diskOuterRadius'):
                    pm.diskOuterRadius = valFloat
                elif line.startswith('diskInnerRadius'):
                    pm.diskInnerRadius = valFloat
                elif line.startswith('diskThick'):
                    pm.diskThick = valFloat
                elif line.startswith('diskBright'):
                    pm.diskBright = valFloat
                elif line.startswith('diskInclination'):
                    pm.diskInclination = valFloat
                elif line.startswith('diskOrbitAngle'):
                    pm.diskOrbitAngle = valFloat
                elif line.startswith('planRadius'):
                    pr = valFloat
                    pm.planRadius = pr*RJupiterToRSun
                    self.planRadEntry.config(textvariable=tkinter.StringVar(value=pr))
                elif line.startswith('planOrbitRadius'):
                    pm.planOrbitRadius = valFloat
                    self.planOrbitRadiusEntry.config(textvariable = tkinter.StringVar(value=pm.planOrbitRadius))
                elif line.startswith('planBright'):
                    pm.planBright = valFloat
                    self.planBrightEntry.config(textvariable = tkinter.StringVar(value=pm.planBright))
                elif line.startswith('starRadius'):
                    pm.starRadius = valFloat
                    self.starRadEntry.config(textvariable=tkinter.StringVar(value=pm.starRadius))
                elif line.startswith('LDC'):
                    pm.LDC = valFloat
                    self.LDCEntry.config(textvariable=tkinter.StringVar(value=pm.LDC))
                else:
                    print("unrecognized line",line)
            inFile.close()
            print("Loaded system model inputs")
            self.plotLDC()
            self.resetCanvas2()
    def initialize(self):
        self.title("ExoTransit_Visualization V2.1")
        self.geometry("1100x800")
        # scaling from brightness to pixel value = uint16/2
        self.pixelScale = 65535
        # define arrays for the orbital angle and the total luminosity
        pm.systemLuminosity = pm.starLuminosity + pm.planetLuminosity
        self.luminosities = np.zeros(pm.nAngleBins, float)
        # the starting angle is chosen so that the planet is at the minimum X position, at -orbitalRadius
        step = 360./pm.nAngleBins
        self.angles = np.arange(0.,360.,step,float)
       # place the entry fields
        ypo = 30
        xpo1 = 0
        xpo2 = 220
        xpo3 = 320
        self.lbl0 = tkinter.Label(self, text = "Star radius (R_star)")
        self.lbl01 = tkinter.Label(self, text = "in units of R_Sun")
        self.starRadEntry = tkinter.Entry(bd=3, textvariable = tkinter.StringVar(value=pm.starRadius), width = 8)
        self.starRadEntry.bind("<Return>",self.validateStarRadius)
        self.lbl0.place(x=xpo1,y=ypo)
        self.starRadEntry.place(x=xpo2, y=ypo)
        self.lbl01.place(x=xpo3,y=ypo)
        ypo += 25
        self.lbl1 = tkinter.Label(self, text = "Limb-Darkening Coefficient (LDC)")
        self.lbl11 = tkinter.Label(self, text = "0 (none) < LDC < 1")
        self.LDCEntry = tkinter.Entry(bd=3, textvariable = tkinter.StringVar(value=pm.LDC), width = 8)
        self.LDCEntry.bind("<Return>",self.validateLDC)
        self.lbl1.place(x=xpo1, y=ypo)
        self.LDCEntry.place(x=xpo2, y=ypo)
        self.lbl11.place(x=xpo3, y=ypo)
        #
        ypo += 25
        self.lbl3 = tkinter.Label(self, text = "Planet radius (R_planet)")
        self.lbl31 = tkinter.Label(self, text = "R_Jupiter. (R_Jupiter= 0.1 R_Sun)")
        self.planRadEntry = tkinter.Entry(bd=3, textvariable = tkinter.StringVar(value=pm.planRadius/RJupiterToRSun), width = 8)
        self.planRadEntry.bind("<Return>",self.validatePlanRadius)
        self.lbl3.place(x=xpo1, y=ypo)
        self.planRadEntry.place(x=xpo2, y=ypo)
        self.lbl31.place(x=xpo3, y=ypo)
        #
        ypo += 25
        self.lbl2 = tkinter.Label(self, text = "Planet orbital radius")
        self.lbl21 = tkinter.Label(self, text = "R_Star. Note: R_Mercury = 83 R_Sun")
        self.planOrbitRadiusEntry = tkinter.Entry(bd=3, textvariable = tkinter.StringVar(value=pm.planOrbitRadius), width = 8)
        self.planOrbitRadiusEntry.bind("<Return>",self.validatePlanOrbitRadius)
        self.lbl2.place(x=xpo1, y=ypo)
        self.planOrbitRadiusEntry.place(x=xpo2, y=ypo)
        self.lbl21.place(x=xpo3, y=ypo)
        #
        ypo += 25
        self.lbl4 = tkinter.Label(self, text = "Planet central brightness")
        self.lbl41 = tkinter.Label(self, text = "0 for planet or >0 for a binary companion star\u2020")
        self.planBrightEntry = tkinter.Entry(bd=3, textvariable = tkinter.StringVar(value=pm.planBright), width = 8)
        self.planBrightEntry.bind("<Return>",self.validatePlanBright)
        self.lbl4.place(x=xpo1, y=ypo)
        self.planBrightEntry.place(x=xpo2, y=ypo)
        self.lbl41.place(x=xpo3, y=ypo)
        #
        ypo += 25
        self.lbl5 = tkinter.Label(self, text = "Planet orbit inclination angle")
        self.lbl51 = tkinter.Label(self, text = "degrees")
        self.planInclEntry = tkinter.Entry(bd=3, textvariable = tkinter.StringVar(value=pm.inclination), width = 8)
        self.planInclEntry.bind("<Return>",self.validateInclination)
        self.lbl5.place(x=xpo1, y=ypo)
        self.planInclEntry.place(x=xpo2, y=ypo)
        self.lbl51.place(x=xpo3, y=ypo)
        # noise
        ypo += 25
        self.lbl6 = tkinter.Label(self, text = "Noise rms")
        self.lbl61 = tkinter.Label(self, text = "(about 0.0005 for real NGTS data)")
        self.noiseEntry = tkinter.Entry(bd=3, textvariable = tkinter.StringVar(value=pm.noise), width = 8)
        self.noiseEntry.bind("<Return>",self.validateNoise)
        self.lbl6.place(x=xpo1, y=ypo)
        self.noiseEntry.place(x=xpo2, y=ypo)
        self.lbl61.place(x=xpo3, y=ypo)
        # run button
        ypo += 30
        self.runBut=tkinter.Button(self, text='Run', command=self.run)
        self.runBut.place(x=xpo1, y=ypo)
        # transit zoom
        self.zoomText = tkinter.StringVar()
        self.zoomText.set("Zoom")
        self.zoombutton=tkinter.Button(self, textvariable=self.zoomText, command=self.zoomTransit)
        self.zoombutton.place(x=xpo1+50, y=ypo)
        # add features like dust disks and stellar luminosity asymmetry
        self.featureText = tkinter.StringVar()
        self.featureText.set("Add feature")
        self.featurebutton=tkinter.Button(self, textvariable=self.featureText, command=self.addFeature)
        self.featurebutton.place(x=xpo1+200, y=ypo)
        # add the export button
        self.exportText = tkinter.StringVar()
        self.exportText.set("Export data")
        self.exportText=tkinter.Button(self, textvariable=self.exportText, command=self.export)
        self.exportText.place(x=xpo1+300, y=ypo)
        # quit button at the top
        self.b3=tkinter.Button(self, text="Quit", command=self.appQuit)
        self.b3.pack(anchor="nw")
        self.b4=tkinter.Button(self, text="Load inputs", command=self.loadInputs)
        self.b4.place(x=190,y=0)
        # read in system inputs
        # add some informational text fields
        ypo += 30
        daggerTxt = tkinter.Label(self, text = "\u2020 Note: Star central brightness, B(r=0) = 0, by definition.  (Luminosity = \u222B B(r) LDC(r) dA)")
        daggerTxt.place(x=xpo1, y=ypo)
        ypo += 20
        self.starLumTxt = tkinter.Label(self, text = "Star luminosity = UNDEFINED")
        self.starLumTxt.place(x=xpo1, y=ypo)
        self.planLumEntryTxt = tkinter.Label(self, text = "Planet luminosity = UNDEFINED")
        self.planLumEntryTxt.place(x=xpo1+200, y=ypo)
        ypo += 20
        self.systLumEntryTxt = tkinter.Label(self, text = "System luminosity(0 deg) = UNDEFINED")
        self.systLumEntryTxt.place(x=xpo1, y=ypo)
        ypo += 20
        self.ang90LumTxt = tkinter.Label(self, text = "System luminosity(90 deg) = UNDEFINED")
        self.ang90LumTxt.place(x=xpo1, y=ypo)
        self.dipTxt = tkinter.Label(self, text = "Normalized luminosity dip = UNDEFINED")
        self.dipTxt.place(x=xpo1+250, y=ypo)
        # add the Brightness vs Radius figure
        Fig = Figure(figsize=(5.5,3.5),dpi=100)
        FigSubPlot = Fig.add_subplot(111)
        FigSubPlot.set_xlabel("Radius Fraction")
        FigSubPlot.set_ylabel("LDC Brightness")
        FigSubPlot.set_xlim(0, 1.1)
        FigSubPlot.set_ylim(0, 1.1)
        Fig.tight_layout()
        ra = np.linspace(0, 1, num=50, endpoint=True)
        ya = np.zeros(50)
        indx = 0
        for r in ra:
            r2 = r * r
            ya[indx] = self.limbDarkFactor(r2)
            indx += 1
        self.line1, = FigSubPlot.plot(ra,ya,"r-")
        self.canvas = matplotlib.backends.backend_tkagg.FigureCanvasTkAgg(Fig, master=self)
        self.canvas.get_tk_widget().pack(side=tkinter.TOP,anchor="ne", expand=True)
        # create the system cartoon and transit figure in another canvas
        Fig2 = Figure(figsize=(10,5))
        self.canvas2 = matplotlib.backends.backend_tkagg.FigureCanvasTkAgg(Fig2, master=self)
        self.F2Cartoon = Fig2.add_subplot(211)
        self.updateStarImage(pm)
        self.updatePlanetImage(pm)
        self.updateSystemImage(pm,0)
        self.ims = self.F2Cartoon.imshow(pm.systemImage, cmap="gray", vmin=0, vmax=65535)
        self.F2Cartoon.axis("off")
        self.F2Transit = Fig2.add_subplot(212)
        self.F2Transit.autoscale(True,axis='both',tight=True)
        angla = np.linspace(0, 360, num=60, endpoint=True)
        intena = np.zeros(60)
        self.line2, = self.F2Transit.plot(angla,intena,"o")
        self.canvas2.get_tk_widget().pack(side=tkinter.BOTTOM, expand=True)
        self.canvas2.draw()
        self.UpdateRunInfo()
        # finish up
        self.resizable(True, True)
#        self.update()
    def UpdateRunInfo(self):
        # Updates the detailed results below the Run button
        self.starLumTxt.configure(text="Star luminosity = {:.4f}".format(pm.starLuminosity))
        self.planLumEntryTxt.configure(text="Planet luminosity = {:.4f}".format(pm.planetLuminosity))
        lumMax = np.double(pm.starLuminosity + pm.planetLuminosity)
        self.systLumEntryTxt.configure(text="System luminosity(0 deg) = {:.4f}".format(lumMax))
        self.updateSystemImage(pm,90.)
        pm.luminosity90 = pm.systemLuminosity
        self.ang90LumTxt.configure(text="System luminosity(90 deg) = {:.4f}".format(pm.luminosity90))
        dip = np.double((lumMax - pm.luminosity90) / lumMax)
        self.dipTxt.configure(text="Normalized luminosity dip = {:.4f}".format(dip))
        # move the planet back to 0
        self.updateSystemImage(pm,0.)
    def updateStarImage(self, pm):
        # there is no need to initialize the star image to 0's because the radius is always 1
        centerRow = int(pm.starImage.shape[0]/2)
        centerCol = int(pm.starImage.shape[1]/2)
        maxBright = 1.
        if pm.planBright > maxBright: maxBright = pm.planBright
        if pm.diskBright > maxBright: maxBright = pm.diskBright
        self.pixelScale = 65535./maxBright
        # determine the step size per pixel for a circle of radius = 1
        step = pm.stepSize
        xs = np.double(-1)
        while xs < 1:
            ymax = np.sqrt(1 - xs*xs)
            ys = -ymax
            xPixel = int(centerRow + xs/step)
            while ys < ymax:
                # distance^2 from star center
                rs2 = xs*xs + ys*ys
                yPixel = centerCol + int(ys/step)
                pixelValue = int(self.pixelScale*self.limbDarkFactor(rs2))
                if pixelValue <= 0: pixelValue = 1
                elif pixelValue > 65535: pixelValue = 65535
                pm.starImage[xPixel][yPixel] = np.uint16(pixelValue)
                ys += step
            xs += step
        # define the luminosity normalization if it hasn't been done
        if pm.luminosityNorm < 0:
            pm.luminosityNorm = np.double(1/(self.pixelScale*np.double(np.count_nonzero(pm.starImage))))
        pm.starLuminosity = np.double(np.sum(pm.starImage))*pm.luminosityNorm
    def updatePlanetImage(self,pm):
        # reset to 0
        pm.planetImage.fill(0)
        maxBright = 1.
        if pm.planBright > maxBright: maxBright = pm.planBright
        if pm.diskBright > maxBright: maxBright = pm.diskBright
        self.pixelScale = 65535./maxBright
        centerRow = pm.planetImage.shape[0]/2
        centerCol = pm.planetImage.shape[1]/2
        step = pm.stepSize
        # scale the planet radius to the star radius
        planRad = pm.planRadius / pm.starRadius
        pr2 = planRad * planRad
        xp = -planRad
        isStar = (pm.planBright > 0)
        while xp < planRad:
            xPixel = int(centerRow + xp/step)
            ymax = np.sqrt(pr2 - xp*xp)
            yp = -ymax
            while yp < ymax:
                yPixel = int(centerCol + yp/step)
                # make the planet not-totally-dark to distinguish it from a deep-space background
                pixelValue = 1
                if isStar:
                    # apply limb-darkening to the 2ndry star
                    rp2 = xp*xp + yp*yp
                    pixelValue = int(self.pixelScale*self.limbDarkFactor(rp2/pr2)*pm.planBright)
                    if pixelValue <= 0: pixelValue = 1
                    elif pixelValue > 65535: pixelValue = 65535
                pm.planetImage[xPixel][yPixel] = np.uint16(pixelValue)
                yp += step
            xp += step
        # Add a dustdisk?
        if pm.featureState == 1 and pm.diskOuterRadius > 0:
            # construct a 3D rotation matrix using Euler angles.
            # the angle wrt the planet orbital plane, alpha
            alpha = pm.diskOrbitAngle * np.pi/180.
            ca = np.cos(alpha)
            sa = np.sin(alpha)
            # the inclination angle of the planet orbital plane, beta
            beta = (90 - pm.diskInclination) * np.pi/180.
            cb = np.cos(beta)
            sb = np.sin(beta)
            # the third rotation angle, gamma, is around the disk axis (Z) which isn't relevant
            cg = 1.0
            sg = 0.0
            # Ref https://en.wikipedia.org/wiki/Rotation_matrix
            rotMat = np.array([[cb*cg, sa*sb*cg-ca*sg, ca*sb*cg+sa*sg],
                               [cb*sg, sa*sb*sg+ca*cg, ca*sb*sg-sa*cg],
                               [-sb, sa*cb, ca*cb]])
            # find the image pixel loop bounds
            diskRad = pm.diskOuterRadius * pm.starRadius
            dr2 = diskRad*diskRad
            diskIn = pm.diskInnerRadius * pm.starRadius
            di2 = diskIn*diskIn
            xd = -diskRad
            # disk thickness
            zmax = 0.5 * pm.diskThick * diskRad
            if (zmax < step): zmax = step
            # this assumes that the disk is uniformly bright
            pixelValue = int(self.pixelScale*pm.diskBright)
            if pixelValue <= 0: pixelValue = 1
            elif pixelValue > 65535: pixelValue = 65535
            pixelVal16 = np.uint16(pixelValue)
            # iterate over all the pixels in the 3D ring model
            while xd < diskRad:
                ymax = np.sqrt(dr2 - xd*xd)
                yd = -ymax
                while yd < ymax:
                    rd2 = yd*yd + xd*xd
                    if rd2 < di2:
                        # this point is less than the ring inner radius so ignore it
                        yd += step
                        continue
                    # iterate over the disk thickness
                    zd = -zmax
                    ringInFront = bool(xd<0)
                    while zd < zmax:
                        # Rotate the position (xd, yd, zd) in 3D
                        pixelPos = np.array([xd, yd, zd])
                        # and project into 2D
                        rotPixelPos = np.dot(rotMat, pixelPos)
                        # Convert to pixel location offset from the image center
                        xPixel = int(centerRow + rotPixelPos[0]/step)
                        yPixel = int(centerCol + rotPixelPos[1]/step)
                        if pm.planetImage[xPixel][yPixel] == 0:
                            pm.planetImage[xPixel][yPixel] = pixelVal16
                        elif ringInFront:
                            pm.planetImage[xPixel][yPixel] = pixelVal16
                        zd += step
                    yd += step
                xd += step
        # dust disk done
        sum = np.sum(pm.planetImage)
        pixCnt = np.count_nonzero(pm.planetImage)
        # subtract the pixel value = 1 dum = pixCnt and normalize to the number of pixels in the star
        pm.planetLuminosity = np.double(sum - pixCnt)*pm.luminosityNorm
        # normalize to the star radius
        pm.planetLuminosity /= pm.starRadius
    def setAngleRange(self):
        # Finds one or two transits without updating the image arrays
        # Default for zoomState == 0
        if pm.zoomState == 0 or (pm.zoomState == 2 and pm.planetLuminosity == 0):
            self.angles = np.arange(0., 360., 360./pm.nAngleBins)
        else:
            start = -1
            nt = 0
            minAngle = 0.
            maxAngle = 360.
            if pm.planetLuminosity == 0:
                maxAngle = 180
            # adjust the step size so that there will be about 20 occluded points in the transit curve
            step = pm.starRadius/pm.planOrbitRadius
            step /= 20.
            step *= 180. / np.pi
            for angleDeg in np.arange(0.,maxAngle,step):
                angle = angleDeg * np.pi/180.
                xcp = -pm.planOrbitRadius * np.cos(angle)
                ycp = pm.planOrbitRadius * np.sin(angle) * np.tan(pm.inclination*np.pi/180)
                rp2 = xcp*xcp + ycp*ycp
                radCut2 = (1.0+pm.planRadius)*(1.0+pm.planRadius)
                if rp2 < radCut2 and start < 0:
                    # the planet center position is within the star disk
                    start = angleDeg
                elif rp2 > radCut2 and start > 0:
                    nt += 1
                    if nt == pm.zoomState:
                        # pad the angle range
                        pad = 0.5 * (angleDeg - start)
                        start = round(start-pad)
                        end = round(angleDeg+pad)
                        if start == angleDeg:
                            print("Failed to find an expected dip")
                            return
                        step = (end - start)/pm.nAngleBins
                        for bin in np.arange(0,pm.nAngleBins,1):
                            self.angles[bin] = start + bin * step
                        break
                    start = -1
        # update the transit curve
#        print("angle range",self.angles[0],self.angles[pm.nAngleBins-1])
        self.line2.set_xdata(self.angles)
        self.F2Transit.set_xlim([self.angles[0],self.angles[pm.nAngleBins-1]])
        self.resetTransitCurve()
        self.canvas2.draw()
        self.canvas2.flush_events()
    def updateSystemImage(self, pm, angleDeg):
        # updates the image of the star/planet system and sets the system luminosity
        # The first step is to see if the planet image is within the cartoon canvas.
        step = pm.stepSize
        planetInFront = bool(angleDeg < 180)
        angle = angleDeg * np.pi/180.
        # x, y position of the planet center viewed at the inclination angle in
        # coordinate system where RStar = 1.
        xcp = -pm.planOrbitRadius * np.cos(angle)
        # see if the planet center is outside the cartoon window
        planetEdge = abs(xcp - pm.planRadius)
        if planetEdge > 2.5:
            pm.systemImage = pm.starImage.copy()
            pm.systemLuminosity = pm.starLuminosity + pm.planetLuminosity
            return
        ycp = pm.planOrbitRadius * np.sin(angle) * np.tan(pm.inclination*np.pi/180)
        # find the number of pixels to shift the planet position into the system image
        iyShift = int(xcp/step) + int(pm.systemImage.shape[1]/2 - pm.planetImage.shape[1]/2)
        ixShift = int(ycp/step)
        # get a tuple array of non-zero planet pixel indices
        nzPixels = np.nonzero(pm.planetImage)
        maxX = pm.systemImage.shape[0]
        maxY = pm.systemImage.shape[1]
        if planetInFront:
            # copy the star image into the system image
            pm.systemImage = pm.starImage.copy()
            # shift and paste non-zero planet pixels
            for ii in range(nzPixels[0].shape[0]):
                ix = nzPixels[0][ii]
                iy = nzPixels[1][ii]
                ixs = ix+ixShift
                iys = iy+iyShift
                if ixs < 0 or ixs >= maxX: continue
                if iys < 0 or iys >= maxY: continue
                pm.systemImage[ixs][iys] = pm.planetImage[ix][iy]
        else:
            # The planet is behind the star. Clear the system image
            pm.systemImage.fill(0)
            # paste the shifted planet image pixels
            for ii in range(nzPixels[0].shape[0]):
                ix = nzPixels[0][ii]
                iy = nzPixels[1][ii]
                ixs = ix+ixShift
                iys = iy+iyShift
                if ixs < 0 or ixs >= maxX: continue
                if iys < 0 or iys >= maxY: continue
                pm.systemImage[ixs][iys]= pm.planetImage[ix][iy]
            # Get a list of the non-zero star pixels
            nzPixels = np.nonzero(pm.starImage)
            # paste them into the system image
            for ii in range(nzPixels[0].shape[0]):
                ix = nzPixels[0][ii]
                iy = nzPixels[1][ii]
                pm.systemImage[ix][iy]= pm.starImage[ix][iy]
        pm.systemLuminosity = np.double(np.sum(pm.systemImage))*pm.luminosityNorm
    def run(self):
        # Run the model and update the screen with the results.
        if pm.zoomState > 1 and pm.planetLuminosity == 0:
            return
        self.runBut.configure(text= '...')
        self.resetTransitCurve()
        sysCnt = np.count_nonzero(pm.starImage)+np.count_nonzero(pm.planetImage)
        isPlanet = bool(pm.planetLuminosity == 0)
        self.luminosities.fill(pm.systemLuminosity)
        # guess at the transit curve Y limits
        ymax = pm.systemLuminosity + 0.001
        self.updateSystemImage(pm,90.)
        ymin = pm.systemLuminosity
        if pm.planetLuminosity > 0:
            # find low value for a bright companion
            self.updateSystemImage(pm,270.)
            if pm.systemLuminosity < ymin:
                ymin = pm.systemLuminosity
        ymin -= 0.001
        self.F2Transit.set_ylim(ymin, ymax)
        self.canvas2.draw()
        self.canvas2.flush_events()
        cnt = 0
        for angle in self.angles:
            if isPlanet and angle > 180:
                break
            self.updateSystemImage(pm,angle)
            self.ims.set_data(pm.systemImage)
            self.luminosities[cnt] = pm.systemLuminosity
            self.line2.set_ydata(self.luminosities)
            self.canvas2.draw()
            self.canvas2.flush_events()
            cnt += 1
        self.UpdateRunInfo()
        self.runBut.configure(text='Run')
        self.updateTransitCurve()
    def resetCanvas2(self):
        self.updateStarImage(pm)
        self.updatePlanetImage(pm)
        self.updateSystemImage(pm, 0.)
        self.ims.set_data(pm.systemImage)
        self.UpdateRunInfo()
        self.resetTransitCurve()
        self.setAngleRange()
    def resetTransitCurve(self):
        self.luminosities.fill(0.)
        self.updateTransitCurve()
    def updateTransitCurve(self):
        # update the transit curve
        self.line2.set_xdata(self.angles)
        transitY = self.luminosities
        if pm.noise > 0:
            transitY = self.luminosities + np.random.normal(0, pm.noise, pm.nAngleBins)
        self.line2.set_ydata(transitY)
        maxLum = np.amax(transitY)
        minLum = np.amin(transitY)
        if (maxLum == minLum):
            maxLum = minLum + 0.01
        yRange = maxLum - minLum
        loY = minLum - 0.1*yRange
        hiY = maxLum + 0.1*yRange
        self.F2Transit.set_ylim(loY, hiY)
        self.canvas2.draw()
        self.canvas2.flush_events()
    def plotLDC(self):
        # Get the entries in the text fields and validate
        ra = np.linspace(0, 1, 50)
        ya = np.zeros(50)
        indx = 0
        for r in ra:
            r2 = r * r
            ya[indx] = self.limbDarkFactor(r2)
            indx += 1
        self.line1.set_ydata(ya)
        self.canvas.draw()
        self.canvas.flush_events()
    def addFeature(self):
        pm.featureState = (pm.featureState+1) % 3
        if pm.featureState == 0:
            self.featureText.set("No feature")
        elif pm.featureState == 1:
            if pm.diskOuterRadius == 0.:
                print("Dust disk model inputs haven't been loaded.")
            else:
                self.featureText.set("Dust disk")
        else:
            self.featureText.set("Future feature")
        self.resetCanvas2()
    def export(self):
        outFile = open("transit.csv","w")
        cnt = 0
        for angle in self.angles:
            outFile.write("{:.2f}".format(angle))
            outFile.write(",")
            outFile.write(" {:.5f}".format(self.luminosities[cnt]))
            outFile.write("\n")
            cnt += 1
        outFile.close()
        print("Exported transit curve to transit.csv")
    def zoomTransit(self):
        pm.zoomState = (pm.zoomState+1) % 3
        if pm.zoomState == 0:
            self.zoomText.set("Zoom")
        elif pm.zoomState == 1:
            self.zoomText.set("Zoom Primary transit")
        elif pm.zoomState == 2:
            if pm.planetLuminosity > 0:
                self.zoomText.set("Zoom 2ndry transit")
            else:
                print("No secondary transit is expected")
                pm.zoomState = 0
                self.zoomText.set("Zoom")
        # Set the angle range for the current zoom state
        self.setAngleRange()
    def validateInclination(self, event):
        # inclination angle
        try:
            inclin = float(self.planInclEntry.get())
        except:
            self.planInclEntry.config(bg="red")
            return False
        if (inclin < 0 or inclin > 90):
            print("The inclination angle must be between 0 and 90 degrees")
            self.planInclEntry.config(bg="red")
            return False
        maxInclination = 180*math.atan((pm.starRadius+pm.planRadius)/pm.planOrbitRadius)/np.pi
        if (inclin > maxInclination):
            mess = "The inclination angle should be less than {:.0f} degrees to see a transit".format(maxInclination)
            print(mess)
            self.planInclEntry.config(bg="red")
            return False
        self.planInclEntry.config(bg="white")
        pm.inclination = inclin
        self.planInclEntry.config(textvariable=tkinter.StringVar(value=pm.inclination))
        self.resetTransitCurve()
        return True
    def validateNoise(self, event):
        # noise
        try:
            noise = float(self.noiseEntry.get())
        except:
            self.noiseEntry.config(bg="red")
            return False
        if (noise < 0 or noise > 0.005):
            print("The noise rms must be between 0 and 0.005")
            self.noiseEntry.config(bg="red")
            return False
        self.noiseEntry.config(bg="white")
        pm.noise = noise
        self.noiseEntry.config(textvariable=tkinter.StringVar(value=pm.noise))
        self.updateTransitCurve()
        return True
    def validatePlanBright(self, event):
        # planet brightness
        try:
            planBright = float(self.planBrightEntry.get())
        except:
            self.planBrightEntry.config(bg="red")
            return False
        if (planBright < 0):
            print("The planet brightness can't be < 0")
            self.planBrightEntry.config(bg="red")
            return False
        pm.planBright = planBright
        self.planBrightEntry.config(bg="white")
        self.planBrightEntry.config(textvariable=tkinter.StringVar(value=pm.planBright))
        self.resetCanvas2()
        return True
    def validatePlanOrbitRadius(self, event):
        # orbital radius
        try:
            orbRad = float(self.planOrbitRadiusEntry.get())
        except:
            self.planOrbitRadiusEntry.config(bg="red")
            return False
        # ensure that the planet isn't inside the star
        if orbRad < pm.starRadius + pm.planRadius:
            print("The planet orbital radius must be > R_Star + R_Planet")
            self.planOrbitRadiusEntry.config(bg="red")
            return False
        pm.planOrbitRadius = orbRad
        self.planOrbitRadiusEntry.config(bg="white")
        self.planOrbitRadiusEntry.config(textvariable=tkinter.StringVar(value=orbRad))
        self.resetCanvas2()
        return True
    def validatePlanRadius(self, event):
        # planet radius
        try:
            planRad = float(self.planRadEntry.get())
        except:
            self.planRadEntry.config(bg="red")
            return False
        planRad *= RJupiterToRSun
        if planRad < 0 or planRad > pm.starRadius:
            print("The planet radius must in the range 0 to R_star")
            self.planRadEntry.config(bg="red")
            return False
        self.planRadEntry.config(bg="white")
        pm.planRadius = planRad
        self.planRadEntry.config(textvariable=tkinter.StringVar(value=pm.planRadius/RJupiterToRSun))
        self.resetCanvas2()
        return True
    def validateStarRadius(self, event):
        try:
            srad = float(self.starRadEntry.get())
        except:
            self.starRadEntry.config(bg="red")
            return False
        # next ensure that the radius value is valid
        if (srad < 0.1):
            print("The star radius should be greater than 0.1 R_Sun")
            self.starRadEntry.config(bg="red")
            return False
        else:
            self.starRadEntry.config(bg="white")
        pm.starRadius = srad
        self.starRadEntry.config(textvariable=tkinter.StringVar(value=srad))
        self.resetCanvas2()
        return True
    def validateLDC(self, event):
        try:
            LDC = float(self.LDCEntry.get())
        except:
            self.LDCEntry.config(bg="red")
            return False
        if (LDC < 0 or LDC > 1):
            print("The limb-darkening coefficient must be between 0 and 1")
            self.LDCEntry.config(bg="red")
            return False
        else:
            self.LDCEntry.config(bg="white")
        print("LDC old",pm.LDC,"new",LDC)
        pm.LDC = LDC
        self.LDCEntry.config(textvariable=tkinter.StringVar(value=LDC))
        self.resetCanvas2()
        self.plotLDC()
        return True
    def limbDarkFactor(self, radialDistance2):
        # calculate the limb-darkening factor at given radial distance^2 with a given limbDark factor.
        if (radialDistance2 < 0 or radialDistance2 > 1): return 0
        arg = 1 - radialDistance2
        arg = 1 - pm.LDC * (1 - np.sqrt(arg))
        return arg

if __name__ == "__main__":
    MainWindow = App_Window(None)
    '''
    str = "Hit the return key after changing a variable to update the model.\n\n"
    str += "Hit the Run button after you are done to start the animation.\n\n"
    str += "Hitting the return key in the 'noise rms' field randomizes the noise."
    tkinter.messagebox.showinfo(message=str)
    '''
    MainWindow.mainloop()


text:
from __global__ import *
import matplotlib
try:
    matplotlib.use('pdf')
except Exception:
    pass
import matplotlib.pyplot as plt
import numpy as np
import scipy.constants
import time
import logging
import pickle

log = logging.getLogger(LOG_NAME+'.'+__name__)
log.setLevel(logging.INFO)

def safeGet(key,dataCont):

    if key in dataCont.keys():
       return np.array(dataCont[key])
    else:

        log.warning('Key %s not found for worker %s' %(key,dataCont['workerId']))
        return np.array([])

def processData(conf,data,startTime = 0,log_folder = None):
    saveMsg = []

    if log_folder:
        # dump all the data in one file
        filetimeStamp = time.strftime("%Y_%m_%d_%H_%M_", time.gmtime(startTime))
        dataName = '%s/%sstat_data%s' %(log_folder,
                                            str(filetimeStamp),
                                            '.dat')
        with open(dataName,'wb') as f:
            pickle.dump(data,f,protocol=pickle.HIGHEST_PROTOCOL)
            f.close()

    for d in data:
        saveMsg.append(processWorkerData(conf,d,startTime,log_folder))
    return saveMsg

def processWorkerData(conf,dataCont,startTime = 0, log_folder = None):

    workerId = safeGet('workerId',dataCont)
    doppler = safeGet('doppler',dataCont)
    doppler_std = safeGet('doppler_std',dataCont)
    timestamp = safeGet('timestamp',dataCont)
    spSymEst = safeGet('spSymEst',dataCont)
    SNR = safeGet('SNR',dataCont)
    numSyncSig = safeGet('numSyncSig',dataCont)
    packetSuc = safeGet('packetSuc',dataCont)
    packetFail = safeGet('packetFail',dataCont)
    numSlaves = safeGet('numSlaves',dataCont).astype(np.int)

    baudRate = safeGet('baudRate',dataCont)
    if len(baudRate) > 0:
        baudRate = baudRate[0]

    # We need to find the proper frequency for the radio.  If not found, use a default value (results in wrong rangerate estimates)
    keys = [k for k in conf['Radios']['Rx'].keys()]
    keyNames = [k.split('-')[0] for k in keys]
    Fc = None
    for i,k in enumerate(keyNames):
        if k in str(workerId):
            print('true')
            radioConf = conf['Radios']['Rx'][keys[i]]
            Fc = radioConf['frequency_Hz']
            samplesPerSlice = 2**conf['GPU'][radioConf['CUDA_settings']]['blockSize']-2**conf['GPU'][radioConf['CUDA_settings']]['overlap']
            log.info('Rx frequency for {} found as {}'.format(workerId,Fc))
            break
    if Fc == None:
        Fc = 1e8
        log.warning('No Rx frequency found for {}. {}'.format(workerId,Fc))

    try:
        rangerateMax = conf['Radios']['rangeRateMax']
    except Exception as e:
        log.warning('Could not find rangeRateMax in log, reason:\n{}'.format(e))
        rangerateMax = 12000

    log.debug('rangerateMax read to {}'.format(rangerateMax))

    spSym = np.mean(spSymEst)
    pkts = packetSuc + packetFail

    plt.ioff()                  # turn of interactive mode such that we can save the figure without showing
    Nplot = len(doppler)
    timeX = timestamp-startTime # makes all workers use the same reference time
    f = plt.figure()
    rangeRate = doppler*scipy.constants.speed_of_light/Fc
    ax1 = f.add_subplot(311)
    ax1.plot(timeX,rangeRate[:Nplot],label='Doppler')
    # ax1.plot(timeX,txFreqOffset[:Nplot],'--',label='Tx Freq offset') # estimated frequency offset
    tmp = np.where(pkts[:Nplot]>0)[0]
    tmpSuc = np.where(packetSuc[:Nplot]>0)[0]
    tmpErr = np.where(packetFail[:Nplot]>0)[0]

    packetNumVotes = numSlaves[tmp]
    colors = ['#800000','#ff0000','#ff6600','#ff9933','#ffcc66','#ffff66']
    # plot number of voters for each packet in a different colour
    if len(packetNumVotes) > 0:
        for i in range(int(np.max(packetNumVotes)+1)):
            tmpV = tmp[np.where(packetNumVotes==i)[0]]
            if len(tmpV) > 0:
                ax1.plot(timeX[tmpV],rangeRate[tmpV],'o',color=colors[i],label='packet {} votes'.format(i+1),markersize=2)

    ax1.plot(timeX[tmpSuc],rangeRate[tmpSuc],'kx',label='decoded',markersize=1)
    ax1.set_title('%s -- Doppler' %(workerId), fontsize='x-small')
    ax1.grid()
    ax1.legend(loc='lower center', mode='expand', fontsize='xx-small',ncol=4)
    ax1.set_ylabel('Range rate [m/s]', fontsize='x-small')
    # ax1.set_xlabel('time [s]')
    ax1.tick_params(axis = 'both', which = 'major',labelsize = 7)
    ax1.tick_params(axis = 'both', which = 'minor',labelsize = 5)

    rangerateStart = -1.35 * rangerateMax
    rangerateStop = 1.05 * rangerateMax
    ax1.set_ylim(rangerateStart,rangerateStop) # fix the limits of the plot

    ## This plot shows the number of sync sigs and voters
    ax3 = f.add_subplot(312,sharex=ax1)
    fixSyncsigAxis = False # hardcodes the axes in the doppler plot
    if 'Stats' in conf.keys():
        if 'FixSyncsigAxis' in conf['Stats'].keys():
            fixSyncsigAxis = conf['Stats']['FixSyncsigAxis']

    ax3.bar(timeX,numSyncSig[:Nplot])
    if fixSyncsigAxis:
        ax3.set_ylim([0,baudRate/8*1.1])
    ax3b = ax3.twinx()
    ax3b.plot(timeX,numSlaves+1,'.',color='red',markersize=1) # #003366
    ax3b.set_ylabel('number of votes', fontsize='x-small')
    ax3b.set_ylim([0,8])
    ax3b.tick_params(axis = 'both', which = 'major',labelsize = 7)
    ax3b.tick_params(axis = 'both', which = 'minor',labelsize = 5)
    # ax3.bar(tmp,numSyncH[tmp],'rx')
    ax3.grid()
    ax3.set_ylabel('Number of sync signals/second', fontsize='x-small')
    ax3.set_xlabel('time [s]', fontsize='x-small')
    ax3.tick_params(axis = 'both', which = 'major',labelsize = 7)
    ax3.tick_params(axis = 'both', which = 'minor',labelsize = 5)
    # ax1.set_xticklabels( () )

    ## Plot SNR
    ax4 = f.add_subplot(313,sharex=ax1)
    ax4.plot(timeX,SNR)
    ax4.set_ylim([-5,50])
    ax4.set_ylabel('SNR [dB]', fontsize='x-small')
    ax4.grid()
    ax4.tick_params(axis = 'both', which = 'major',labelsize = 7)
    ax4.tick_params(axis = 'both', which = 'minor',labelsize = 5)

    ax1.set_xlim(0,timeX[-1])

    if log_folder:

        figPrefix = str(time.strftime("%Y_%m_%d_%H_%M_", time.gmtime(startTime)))

        figName = '%s/%sdoppler_%s.pdf' %(log_folder,
                                          figPrefix,
                                          workerId)
        dataName = '%s/%sdoppler_%s' %(log_folder,
                                       figPrefix,
                                       workerId)

        f.set_size_inches(16.53,11.69) # Save as A3
        # f.tight_layout()
        f.savefig(figName,bbox_inches='tight',format='pdf')
        plt.close(f)

        # save numpy data
        np.savez(dataName,
                 startTime = startTime,
                 timeX = timeX+startTime,
                 doppler = doppler[:Nplot],
                 rawDoppler = doppler[:Nplot],
                 doppler_std = doppler_std[:Nplot],
                 rangeRate =rangeRate[:Nplot],
                 packets = tmp,
                 packetsError = tmpErr,
                 packetsSuccess = tmpSuc,
                 numSyncSignals = numSyncSig[:Nplot],
        )

    log.info('Figure for worker %s saved as %s' %(workerId, figName ))
    return 'Figure for worker %s saved as %s' %(workerId, figName )

max_stars_repo_name: mugpahug/pycu-sdr

max_issues_repo_name: mugpahug/pycu-sdr

max_forks_repo_name: mugpahug/pycu-sdr

reason: import numpy,import scipy


text:
import csv
import datetime
import os
import os.path
import pickle
import re

import numpy as np
import pandas as pd

def storeSportDataParticipant1(filename, data):
  curYear = ""
  curMonth = ""
  curDay = ""
  with open(filename, newline="") as csvfile:
    spamreader = csv.reader(csvfile)
    count = 0
    for row in spamreader:
        count = count + 1
        if count > 1 and len(row):
            if len(row[0]):
                curYear = row[0]
            if len(row[1]):
                curMonth = row[1]
                if len(curMonth) == 1:
                    curMonth = "0" + curMonth
            if len(row[2]):
                curDay = row[2]
                if len(curDay) == 1:
                    curDay = "0" + curDay
            date = curYear + "-" + curMonth + "-" + curDay
            dict = {
                "Walk": "walk",
                "Road Bike": "roadBike",
                "Mt Bike": "mountainBike",
                "Swimming": "swimming",
                "Surfing": "surfing",
                "Climbing": "climbing",
                "Via Ferrata": "viaFerrata",
                "Alpi Ski": "alpiSki",
                "Down Ski": "downSki",
            }
            if row[3] in dict:
                data.loc[date, dict[row[3]]] = 1

            if row[3] == 'Climbing':
                try:
                    data.loc[date, 'climbingDenivelation'] = float(row[7])
                except:
                    data.loc[date, 'climbingDenivelation'] = 0
                try:
                    data.loc[date, 'climbingMaxEffortIntensity'] = float(row[13])
                except:
                    data.loc[date, 'climbingMaxEffortIntensity'] = 0
                try:
                    data.loc[date, 'climbingMeanEffortIntensity'] = float(row[15])
                except:
                    data.loc[date, 'climbingMeanEffortIntensity'] = 0

            if row[3] == 'Swimming':
                try:
                    data.loc[date, 'swimmingKm'] = float(row[6])
                except:
                    data.loc[date, 'swimmingKm'] = 0

  return data

max_stars_repo_name: oliviermirat/Scientizen

max_issues_repo_name: oliviermirat/Scientizen

max_forks_repo_name: oliviermirat/Scientizen


text:
import numpy as np
from PyQt5.QtWidgets import QComboBox, QListWidget

from piafedit.gui.common.template_widget import TemplateWidget
from piafedit.model.libs.filters import erode, dilate, edge_detection
from piafedit.ui_utils import gui_app

operators = {
    'erode': erode,
    'dilate': dilate,
    'edge': edge_detection,
}

class PipelinePanel(TemplateWidget):
    def __init__(self, parent=None):
        super().__init__('pipeline', parent)
        self.operator_combo.addItems(operators.keys())
        self.addBeforeButton.clicked.connect(self.add_before)
        self.addAfterButton.clicked.connect(self.add_after)
        self.removeButton.clicked.connect(self.remove_selected)
        self.clearButton.clicked.connect(self.clear_pipeline)
        self.show()

    def pipeline(self):
        keys = [self.pipeline_list.item(i).text() for i in range(self.pipeline_list.count())]
        ops = [
            operators[key]
            for key in keys
        ]

        def pipe(buffer: np.ndarray):
            for op in ops:
                buffer = op(buffer)
            return buffer

        return pipe

    def clear_pipeline(self):
        while self.pipeline_list.count() > 0:
            self.pipeline_list.takeItem(0)
        self.request_update()

    def add_before(self):
        text = self.operator_combo.currentText()
        index = 0
        for item in self.pipeline_list.selectedItems():
            index = self.pipeline_list.row(item)
        self.pipeline_list.insertItem(index, text)
        self.request_update()

    def add_after(self):
        text = self.operator_combo.currentText()
        index = self.pipeline_list.count()
        for item in self.pipeline_list.selectedItems():
            index = self.pipeline_list.row(item)
        index = min(index+1, self.pipeline_list.count())
        self.pipeline_list.insertItem(index, text)
        self.request_update()

    def remove_selected(self):
        for item in self.pipeline_list.selectedItems():
            self.pipeline_list.takeItem(self.pipeline_list.row(item))
        self.request_update()

    @property
    def pipeline_list(self) -> QListWidget:
        return self.pipelineList

    @property
    def operator_combo(self) -> QComboBox:
        return self.operatorCombo

if __name__ == '__main__':
    with gui_app():
        pipe = PipelinePanel()
        pipe.show()

max_stars_repo_name: flegac/piaf-edit

max_issues_repo_name: flegac/piaf-edit

max_forks_repo_name: flegac/piaf-edit


text:
import numpy as np
import math

def rotation_matrix(axis, theta):
    """
    Return the rotation matrix associated with counterclockwise rotation about
    the given axis by theta radians.

    example
    v = [3, 5, 0]
    axis = [4, 4, 1]
    theta = 1.2

    np.dot(rotation_matrix(axis,theta), v)
    # [ 2.74911638  4.77180932  1.91629719]

    Ry = rotation_matrix([0,1,0], np.deg2rad(45))
    np.dot(Ry, [1,456,1])
    # [  1.41421356e+00   4.56000000e+02  -1.11022302e-16]

    """
    axis = np.asarray(axis)
    axis = axis/math.sqrt(np.dot(axis, axis))
    a = math.cos(theta/2.0)
    b, c, d = -axis*math.sin(theta/2.0)
    aa, bb, cc, dd = a*a, b*b, c*c, d*d
    bc, ad, ac, ab, bd, cd = b*c, a*d, a*c, a*b, b*d, c*d
    return np.array([[aa+bb-cc-dd, 2*(bc+ad), 2*(bd-ac)],
                     [2*(bc-ad), aa+cc-bb-dd, 2*(cd+ab)],
                     [2*(bd+ac), 2*(cd-ab), aa+dd-bb-cc]])

max_stars_repo_name: ggruszczynski/pyVLM

max_issues_repo_name: ggruszczynski/pyVLM

max_forks_repo_name: ggruszczynski/pyVLM