text:
from sklearn.linear_model import LinearRegression, Ridge, RidgeCV
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from scipy.optimize import curve_fit
from sklearn.base import RegressorMixin
from sklearn.preprocessing import FunctionTransformer
from sklearn.model_selection import train_test_split
from sklearn.metrics import explained_variance_score
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LinearRegression
from mpl_toolkits.mplot3d import Axes3D
import Data_Clean_Functions
from scipy import signal
import numpy as np
import matplotlib.pyplot as plt
import pickle
import pandas as pd  # pandas for data handling
import os
import warnings
from sklearn.ensemble import RandomForestRegressor
%matplotlib qt
warnings.simplefilter(action='ignore', category=(FutureWarning, DeprecationWarning))

funcs = Data_Clean_Functions.Funcs()

def low_pass(sr, dt=10, Hz=100):
    'Apply lowpass filter to a data series'
    win = signal.hann(int(Hz*dt)+1)
    pad = np.zeros(len(win)//2) * np.NaN
    sig = signal.convolve(sr, win, mode='valid') / sum(win)
    return pd.Series(np.r_[pad, sig, pad])

# def read_power(folder):
#     df = pd.read_csv(folder+'/fort.9', delim_whitespace=True,
#                      names=["time", "CFL", "power", "phi"])
#     df.drop(df.index[:3], inplace=True)
#     df['power_smooth'] = low_pass(df.power, dt=0.03, Hz=len(df.power))
#     df.Theta = 0.14
#     return df

#
# def average(df, var):
#     return np.trapz(df[var], df.time)/(df.time.iloc[-1]-df.time.iloc[0])
#
#
# plt.figure(figsize=(6.3, 3), dpi=100)
# plt.xlabel(r'$z/L$', fontsize=12)
# plt.ylabel(r'$C_P$', rotation=0, labelpad=15, fontsize=12)
#
# for i in range(4):
#     folder = 'phase/case{:02d}'.format(i)
#     df = read_power(folder)
#     plt.plot(df.time, df.power_smooth/df.Theta**3, label=str(i)+'/8')
#
# bar, amp = average(df, 'power')/df.Theta**3, 1.5748*df.power.mad()/df.Theta**3
# plt.axhline(y=bar, color='grey', ls='--')
# plt.axhline(y=bar+amp, color='grey', ls='-.')
# plt.axhline(y=bar-amp, color='grey', ls='-.')
#
# plt.legend(title=r'$t/T$')
# plt.tight_layout()
# plt.savefig('phase_Cp.png')
# plt.show()
# # %% codecell
# df = pd.read_csv('grid_out.csv').query('r<1')
# df = pd.concat([df, pd.read_csv('grid_out_2.csv')])
#
# fig, axa = plt.subplots(1, 2, figsize=(7, 3), dpi=100, sharey=True)
# for ax in axa:
#     ax.set_yscale('log')
#     ax.set_ylim([1e-2, 1])
# axa[0].set_ylabel(r'$\overline{C_P}$', rotation=0, labelpad=15, fontsize=12)
# axa[0].scatter(df.r, df.bar_p/df.Theta**3, c=180/np.pi*df.Theta)
# axa[0].set_xlabel(r'$r/R$', fontsize=12)
# im = axa[1].scatter(180/np.pi*df.alpha, df.bar_p/df.Theta**3, c=180/np.pi*df.Theta)
# axa[1].set_xlabel(r'$\alpha$ (deg)', fontsize=12)
# fig.colorbar(im, label=r'$\Theta$ (deg)')
# plt.tight_layout()
# plt.savefig('grid_bar_Cp.png')
#
# fig, axa = plt.subplots(1, 2, figsize=(7, 3), dpi=100, sharey=True)
# for ax in axa:
#     ax.set_yscale('log')
#     ax.set_ylim([1e-3, 1])
# axa[0].set_ylabel(r'$\Theta |C_P|$', rotation=0, labelpad=15, fontsize=12)
# axa[0].scatter(df.r, df.amp_p/df.Theta**2, c=180/np.pi*df.Theta)
# axa[0].set_xlabel(r'$r/R$', fontsize=12)
# im = axa[1].scatter(180/np.pi*df.alpha, df.amp_p/df.Theta**2, c=180/np.pi*df.Theta)
# axa[1].set_xlabel(r'$\alpha$ (deg)', fontsize=12)
# fig.colorbar(im, label=r'$\Theta$ (deg)')
# plt.tight_layout()
# plt.savefig('grid_amp_Cp.png')
#
# fig = plt.figure(figsize=(8, 3), dpi=100)
# ax = fig.add_subplot(121, projection='3d')
# ax.scatter(df.r, 180/np.pi*df.alpha, df.bar_p/df.Theta**3, c=180/np.pi*df.Theta)
# ax.set_xlabel(r'$r/R$', fontsize=12)
# ax.set_ylabel(r'$\alpha$ (deg)', fontsize=12)
# ax.set_zlabel(r'$\overline{C_P}$', fontsize=12)
#
# ax = fig.add_subplot(122, projection='3d')
# ax.scatter(df.r, 180/np.pi*df.alpha, df.amp_p/df.Theta**2, c=180/np.pi*df.Theta)
# ax.set_xlabel(r'$r/R$', fontsize=12)
# ax.set_ylabel(r'$\alpha$ (deg)', fontsize=12)
# ax.set_zlabel(r'$\Theta |C_P|$', fontsize=12)
# plt.tight_layout()
# plt.savefig('grid_3D_Cp.png')
# plt.show()

strCWD = os.getcwd()

# ------------------------------------------------------------------------------
# Read in data

df = pd.read_csv(os.path.join(strCWD, 'Data', 'Jul28_Labeled.csv'))
df2 = pd.read_csv(os.path.join(strCWD, 'Data', 'Jul29_Labeled.csv'))
df3 = pd.read_csv(os.path.join(strCWD, 'Data', 'Jul30_Labeled.csv'))

df = df.append(df2)
df = df.append(df3)

df['Utc'] = df['Utc'].astype(str)
df['Utc'] = pd.to_datetime(df['Utc'])
df['Day'] = df.Utc.dt.day

df = df.reset_index().drop('index', axis=1)
index = list(df[(df['Segment'] == 9) & (df['Day'] == 30)].index)
df.loc[index, 'label'] = 'NR'
index = list(df[(df['Segment'] == 41) & (df['Day'] == 30)].index)
df.loc[index, 'label'] = 'NR'
index = list(df[(df['Segment'] == 5) & (df['Day'] == 28)].index)
df.loc[index, 'label'] = 'NR'

df['abs_Twa'] = abs(df.Twa)
df = df.reset_index().drop('index', axis=1)
df = funcs.manoeuver_label(df, man_len=50)
df.label.unique()

df = df[(df['label'] == 'UW') | (df['label'] == 'DW')]

# need to eliminate angles  and wind speeds that the VPP do not capture
df = df[(df['Twa'] > -150) & (df['Twa'] < 150)]
df = df[(df['Tws'] < 20)]
df = df[df['Bsp'] > 5]
df = df[df['abs_Twa'] > 30]

len(df)
df.head()

# -------------------------------------------------------------------------------
# Pipeline
poly_model = make_pipeline(StandardScaler(),
                           PolynomialFeatures(15),
                           LinearRegression())

rf_model = make_pipeline(StandardScaler(),
                         RandomForestRegressor(n_estimators=100))
# Model based on all the data
X = df[['Tws', 'abs_Twa']]
y = df.Bsp
poly_model.fit(X, y)
rf_model.fit(X, y)
score_poly = explained_variance_score(y, poly_model.predict(X))
score_rf = explained_variance_score(y, rf_model.predict(X))
print('4647 point polynomial fit score:{}'.format(score_poly))
print('4647 point polynomial fit score:{}'.format(score_rf))

# Data slice
# Theta_surf = df.Theta.unique()[5]
# X_surf = df[df.Theta == Theta_surf][['Theta', 'alpha', 'r']]
# y_dat = df[df.Theta == Theta_surf].bar_p/Theta_surf**3
# y_fit = poly_model.predict(X_surf)
# df_slice = df[(df['Tws'] < 10.2) & (df['Tws'] > 9.8)]
df_slice = df
X_surf = df_slice[['Tws', 'abs_Twa']]
y_dat = df_slice.Bsp
y_fit_poly = poly_model.predict(X_surf)
y_fit_rf = rf_model.predict(X_surf)

# Phi: Simplified models based on physical arguments
# X[0]=Theta, X[1]=alpha, X[2]=r

#
# def trig(X, a0=1, a1=1, a2=np.pi/2):
#     return a0*np.sin(a1*X[:, 1])*np.cos(a2*X[:, 2])
#
#
# def hyper(X, a0=1, a1=2, a2=4):
#     return a0*(1+np.tanh(a1*X[:, 1]))*np.tanh(a2*(1-X[:, 2]))

pickleIn = open('Models/ORC_Swan45_Rbf.pickle', 'rb')
rbf_swan = pickle.load(pickleIn)
pickleIn.close()

def ORC(X):
    return rbf_swan(X[:, 0], X[:, 1])

IM = ORC

y_IM = IM(X_surf.values)
score_IM = explained_variance_score(y, IM(X.values))
print('Phi model score:{}'.format(score_IM))

# plot
fig = plt.figure(figsize=(8, 4), dpi=100)
ax = fig.add_subplot(121, projection='3d')
# ax.plot_trisurf(df.Tws, df.abs_Twa, y_fit_poly, alpha=0.8)
# ax.scatter(df.Tws, df.abs_Twa, y_dat, c=y_dat-y_fit_poly, alpha=1)
ax.plot_trisurf(df.Tws, df.abs_Twa, y_fit_rf, alpha=0.8)
ax.scatter(df.Tws, df.abs_Twa, y_dat, c=y_dat-y_fit_rf, alpha=1)
ax.set_xlabel('Tws (Knots)', fontsize=12)
ax.set_ylabel('abs_Twa (deg)', fontsize=12)
ax.set_zlabel('Bsp', fontsize=12)

ax = fig.add_subplot(122, projection='3d')
ax.plot_trisurf(df.Tws, df.abs_Twa, y_IM, color='C1', alpha=0.8)
ax.scatter(df.Tws, df.abs_Twa, y_dat, c=y_dat-y_IM, alpha=1)
ax.set_xlabel('Tws (Knots)', fontsize=12)
ax.set_ylabel('abs_Twa (deg)', fontsize=12)
ax.set_zlabel('Bsp', fontsize=12)
plt.show()
# plt.savefig('surf_3D_Cp.png')
# But the polynomial fit assumes we know all $N=1000$ data points. What if we only have $N=100$... or $N=10$? Lets sub-sample the data and to see how this approach generalizes.
warnings.simplefilter(action='ignore', category=(FutureWarning, DeprecationWarning))

# Make the model "pipeline"
poly_model = make_pipeline(StandardScaler(),
                           PolynomialFeatures(4),
                           #                            LinearRegression()) # overfit until N~1000
                           Ridge())  # best for N>100
#                            RidgeCV()) # over-regularized unless N<50
rf_model = make_pipeline(StandardScaler(),
                         #                            LinearRegression()) # overfit until N~1000
                         RandomForestRegressor(n_estimators=100))
# Evaulate model performance

def subfit(n, model, state=1):  # split -> fit -> test
    X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=n, random_state=state)
    model.fit(X_train, y_train)
    y_fit = model.predict(X_test)
    return {'score': explained_variance_score(y_test, y_fit), 'n': n, 'state': state}

def evaluate(model):  # loop through the data set sizes
    return pd.DataFrame([subfit(n, model, state)
                         for n in np.logspace(1, 3.6, 20).astype(int)
                         for state in range(20)])

poly_results = evaluate(poly_model)
rf_results = evaluate(rf_model)
rf_results
poly_results

fig = plt.figure(figsize=(6, 3), dpi=100)
ax = fig.add_subplot(111)
plt.axhline(y=score_poly, ls='--', c='k')
poly_results.plot.scatter(x='n', y='score', ax=ax, alpha=0.5)
poly_results.groupby('n').median().reset_index().plot(x='n', y='score', ax=ax, label='Polynomial Ridge')
rf_results.plot.scatter(x='n', y='score', ax=ax, alpha=0.5)
rf_results.groupby('n').median().reset_index().plot(x='n', y='score', ax=ax, label='RandomForest')

plt.ylim(0, 1)
plt.xscale('log')
plt.xlim(10, 4000)
plt.legend()
plt.title('Explained Bsp variance')
plt.xlabel('Training size')
plt.tight_layout()
plt.show()

# plt.savefig('poly-ridge.png')
# %% markdown
# If we only use `LinearRegression` the results are terrible until $N\approx 1000$. `Ridge` regression looks pretty good for $N\ge100$. `RidgeCV` adds cross-validation, but doesn't consistently improve results - it just avoid disasters for $N<40$.
#
# Lets fit the simple physics-based $\Phi$ model as a semi-emperical model instead.

# Use the IM as a SemiEmperical model

#
# class SemiEmperical(RegressorMixin):
#     def fit(self, X=None, y=None):
#         self.alpha, _ = curve_fit(IM, X.values, y.values, p0=[1, 1, np.pi/2])
#
#     def predict(self, X=None):
#         # Give back the mean of y, in the same
#         # length as the number of X observations
#         return IM(X.values, *self.alpha)
#
#
# semi_results = evaluate(SemiEmperical())
#
# fig = plt.figure(figsize=(6, 3), dpi=100)
# ax = fig.add_subplot(111)
# plt.axhline(y=score_1K, ls='--', c='k')
# poly_results.plot.scatter(x='n', y='score', ax=ax, alpha=0.5)
# poly_results.groupby('n').median().reset_index().plot(x='n', y='score', ax=ax, label='Polynomial Ridge')
# semi_results.plot.scatter(x='n', y='score', ax=ax, c='C1', alpha=0.5)
# semi_results.groupby('n').median().reset_index().plot(x='n', y='score', ax=ax, c='C1', label='Semi-emperical')
#
# plt.ylim(0, 1)
# plt.xscale('log')
# plt.xlim(10, 1e3)
# plt.legend()
# plt.title('Explained Bsp variance')
# plt.xlabel('Training size')
# plt.tight_layout()
# plt.show()
# plt.savefig('semi-emp.png')
# %% markdown
# So this is much better for low data sizes (at least on average), but it has no real learning capacity.
#
# Let's use $\Phi$ to create a basis instead - simple concatenation is enough in this case.
# %% codecell
# Concatenate input

def Phi_generator(X):
    return np.column_stack([X, IM(X)])
#     return np.column_stack([X,hyper(X)])

# PBLM pipeline
PBLM_Ridge_model = make_pipeline(FunctionTransformer(Phi_generator, validate=True),
                                 PolynomialFeatures(2),
                                 RidgeCV(alphas=[1e-3, 1e-2, 1e-1, 1, 1e2, 1e3]))

PBLM_rf_model = make_pipeline(FunctionTransformer(Phi_generator, validate=True),
                              RandomForestRegressor(n_estimators=50))

PBLM_Ridge_results = evaluate(PBLM_Ridge_model)
PBLM_rf_results = evaluate(PBLM_rf_model)

fig = plt.figure(figsize=(6, 3), dpi=100)
ax = fig.add_subplot(111)
plt.axhline(y=score_1K, ls='--', c='k')
poly_results.plot.scatter(x='n', y='score', ax=ax, alpha=0.5)
poly_results.groupby('n').median().reset_index().plot(x='n', y='score', ax=ax, label='Polynomial Ridge')
rf_results.plot.scatter(x='n', y='score', ax=ax, alpha=0.5)
rf_results.groupby('n').median().reset_index().plot(x='n', y='score', ax=ax, label='Random Forest')
# semi_results.plot.scatter(x='n', y='score', ax=ax, c='C1', alpha=0.5)
# semi_results.groupby('n').median().reset_index().plot(x='n', y='score', ax=ax, c='C1', label='Semi-emperical')
PBLM_Ridge_results.plot.scatter(x='n', y='score', ax=ax, c='C4', alpha=0.5)
PBLM_Ridge_results.groupby('n').median().reset_index().plot(x='n', y='score', ax=ax, c='C4', label='$\Phi$-basis ridge')
PBLM_rf_results.plot.scatter(x='n', y='score', ax=ax, c='C5', alpha=0.5)
PBLM_rf_results.groupby('n').median().reset_index().plot(x='n', y='score', ax=ax, c='C5', label='$\Phi$-basis rf')

plt.ylim(0, 1)
plt.xscale('log')
plt.xlim(10, 2e3)
plt.legend()
plt.title('Explained Bsp variance')
plt.xlabel('Training size')
plt.tight_layout()
# plt.savefig('phi-basis.png')
plt.show()

# So the results are **much** more robust than the black-box (polynomial ridge regression) to sparse data - both the median and variance are vastly improved. But there is also enough learning capacity to achieve a perfect fit when $N>100$. Note that *adaptive* regularization is key to the $\Phi$-basis working well across the full range.

reason: import numpy,from scipy


text:
import os
import pandas as pd
import pickle
import numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.base import TransformerMixin #gives fit_transform method for free

DATASET_PATH = "/home/avin/QM/material_data/CPM/data"

def load_material_data(dataset_path=DATASET_PATH):

	csv_path = os.path.join(dataset_path, "4000_dielectric_materials.csv")
	return pd.read_csv(csv_path)

material_data = load_material_data()

material_data = material_data.apply(pd.to_numeric, errors='coerce')
material_data = material_data.dropna(subset=['Sintering', 'er', 'Qf', 'Tcf'])
material_data = material_data[['Sintering', 'er', 'Qf', 'Tcf']]
# print material_data.info()

train_set, test_set = train_test_split(material_data, test_size=0.2, random_state=42)

train_data = train_set.drop('Sintering', axis=1)
train_labels = train_set['Sintering'].copy()
test_data = test_set.drop('Sintering', axis=1)
test_labels = test_set['Sintering'].copy()

#To handle fit_transform() problem
class MyLabelBinarizer(TransformerMixin):
    def __init__(self, *args, **kwargs):
        self.encoder = LabelBinarizer(*args, **kwargs)
    def fit(self, x, y=0):
        self.encoder.fit(x)
        return self
    def transform(self, x, y=0):
        return self.encoder.transform(x)

class DataFrameSelector(BaseEstimator, TransformerMixin):
	def __init__(self, attribute_names):
		self.attribute_names = attribute_names
	def fit(self, X, y=None):
		return self
	def transform(self, X):
		return X[self.attribute_names].values

attribs = ['er', 'Qf', 'Tcf']
pipeline = Pipeline([
						# ('selector', DataFrameSelector(attribs)),
						('std_scaler', StandardScaler())
						])

pipeline.fit_transform(train_data)

#trained model
pkl_filename = "saved_model_sintering/pickle_model.pkl"

# Load from file
with open(pkl_filename, 'r') as file:
    sintering_model = pickle.load(file)

def predict_sintering(input_df):

	sintering_prediction = sintering_model.predict(pipeline.transform(input_df))
	return sintering_prediction[0]

# predict_sintering(np.array([65, 20000, 35]).reshape(1,-1))
# predict_sintering(np.array([45, 20000, 30]).reshape(1,-1))

ext: py

max_stars_repo_name: ihalage/AIMat

max_issues_repo_name: ihalage/AIMat

max_forks_repo_name: ihalage/AIMat

reason: import numpy


text:
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Fri Oct 19 14:51:38 2018

@author: bernier2, sfli
"""
import numpy as np
import tensorflow as tf

epsf = np.finfo(float).eps
sqrt_epsf = np.sqrt(epsf)

def unit_vector(invec):
    invec = np.atleast_2d(invec)
    npts, dim = invec.shape

    invec[np.isnan(invec)] = 0.  # supress RuntimeWarning for invaild values

    nrm = np.sqrt(np.sum(invec*invec, axis=1))
    zidx = nrm < epsf
    nrm[zidx] = 1.

    outvec = invec*np.tile(1./nrm, (dim, 1)).T
    outvec[zidx, :] = 0.
    return outvec

def ray_plane(rays, rmat, tvec):
    """
    tvec: Is the normal defining the plane.
    rvec
    """
    nhat = rmat[:, 2].reshape(3, 1)

    rays = np.atleast_2d(rays)  # shape (npts, 3)
    npts = len(rays)

    numerator = np.tile(np.dot(tvec, nhat), npts)
    denominator = np.dot(rays, nhat).flatten()
    can_intersect = denominator < 0

    u = np.tile(numerator[can_intersect] / denominator[can_intersect], (3, 1))

    output = np.nan*np.ones_like(rays)
    output[can_intersect, :] = rays[can_intersect] * u.T

    return np.dot(output - tvec, rmat)[:, :2]

def tf_ray_plane(rays, rmat, tvec):
    """
    In general, the plane equation is defined by
    $$\hat{n} \cdot (\vec{x} - \vec{x}_0)  = 0$$

    where $\vec{x}_0$ defines the origin of rotation for the plane, $\hat{n}$
    defines the plane normal.

    For a given ray, $\vec{r}(t) = \vec{p}_0 + \hat{u} t$

    the ray-plane intersection is defined by

    $$t = \frac{\hat{n} \cdot (\vec{x}_0 - \vec{p}_0)} {\hat{n} \cdot \hat{u}} $$

    To parameterize the detector plane, we will give the origin point a
    translation vector.

    This gives us

    $$ \mathbf{R} \hat{n} \cdot (\vec{x} - (\vec{x}_0 + \vec{d}))  = 0$$

    where $\vec{d}$ is the translation vector and $\mathbf{R}$ is the rotation
    of the detector about its origin.

    For this function, we'll take \vec{x_0} to be at the origin. We will define
    the translation and rotation as `tvec` and `rmat`

    :param tvec: Translation vector for the origin of the plane.
    :param rmat: Rotation of the detector plane about its origin.
    :param rays: Ray vector with $\vec{p}_0$ taken to be coordinate origin.

    :return: The list of x, y intersection location of the plane.

    As a convention, for tvec = [0, 0, 0], and R = np.eye(3), the detector
    plane is located at the origin, with its normal pointing in the [0, 0, 1]
    direction.

    Similarly, by convention, the origin of the ray, $\vec{p}_0$ is taken to
    be the origin as well.

    """

    # TODO(Frankie): Static type check to ensure that rmat, tvec, and ray
    # are already type, `tf.tensor`.

    # Detector normal, rotated about its origin.
    nhat = rmat[:, 2]

    numerator = tf.broadcast_to(tf.tensordot(tvec, nhat, axes=1), [tf.shape(rays)[0]])
    denominator = tf.tensordot(rays, nhat, axes=1)
    tmp = numerator / denominator

    cannot_intersect_idx = tf.where(denominator > 0)
    # Equivilent of...
    # output[cannot_intersect, :] = np.nan
    tmp = tf.tensor_scatter_nd_update(tmp,
                                      cannot_intersect_idx,
                                      tf.ones(cannot_intersect_idx.shape[0]) * np.nan)
    u = tf.broadcast_to(tmp, [3, tf.shape(rays)[0]])
    output = rays * tf.transpose(u)

    # tf.tensordot with axis [[1], [0]] is just matrix multiply.
    # We only take the x, y, the detector plane coordinates.
    return tf.tensordot(output - tvec, rmat, axes=[[1],[0]])[:, :2]

def test_ray_plane():
    EPS = 0.001

    import numpy as np
    num_samples = 100
    random_angles = np.random.uniform(-np.pi, np.pi,
                                      (num_samples, 3)).astype(np.float32)

    rmats = rotation_matrix_3d.from_euler(random_angles)

    for i in np.arange(num_samples):
        classic_result = ray_plane(rays.numpy(), rmats[i, :, :].numpy(), tvec.numpy())
        tf_result = tf_ray_plane(rays, rmats[i, :, :], tvec).numpy()

        # Assert that the intersected are the same
        classic_intersected = np.logical_not(np.isnan(classic_result).any(axis=1))
        tf_intersected = np.logical_not(np.isnan(tf_result).any(axis=1))

        assert not np.any(np.logical_xor(classic_intersected, tf_intersected)), 'Failed, inteserction check not identical'
        classic_result = classic_result[np.ix_(classic_intersected, [0, 1])]
        tf_result = tf_result[np.ix_(tf_intersected, [0, 1])]

        res = classic_result - tf_result
        # Assert that the different in magnitude is small
        l = np.linalg.norm(res[:])
        print(f'residual {l}')
        assert l < EPS, 'Failed, coorindates check not identical'

    print('Pass')

def gradient_ray_plane(rays, rmat, tvec):
  # Call `GradientTape` with `persistent=True` if we want to
  # keep the gradient after calling it. If that's the case,
  # an explicit destruction by calling `del` on `tape` is needed.
  with tf.GradientTape() as tape:
    tape.watch(rays)
    intersections = tf_ray_plane(rays, rmat, tvec)
  """ Take the derivative of the 'tf_ray_plane' with respect
  to the input variables, 'rays' evaluated at the
  'rays.'

  See: https://www.tensorflow.org/api_docs/python/tf/GradientTape
  """
  return tape.gradient(intersections, rays)

def ray_cylinder(rays, radius, rmat, tvec):
    """
    return points P in ref frame (origin [0, 0, 0]) associated with rays n

    parameterize ray as p = u*d for scalar u and unit vector d

    in ref frame, vector to point on cylinder is l = p - t
    """
    rays = np.atleast_2d(rays)  # shape (npts, 3)
    npts = len(rays)

    output = np.nan * np.ones_like(rays)

    r_0 = rmat[:, 0]
    r_2 = rmat[:, 2]

    a = np.dot(rays, r_0)
    b = np.tile(np.dot(tvec, r_0), npts)
    c = np.dot(rays, r_2)
    d = np.tile(np.dot(tvec, r_2), npts)

    aq = a**2 + c**2
    bq = -2.*(a*b + c*d)
    cq = b**2 + d**2 - radius**2

    discriminant = bq**2 - 4.*aq*cq
    can_intersect = discriminant >= sqrt_epsf
    discriminant[~can_intersect] = np.nan
    aq[aq == 0] = 1.    # nans in discriminant will persist

    # the two roots for the ray vector scaling
    root1 = 0.5*(-bq + np.sqrt(discriminant))/aq
    root2 = -0.5*(bq + np.sqrt(discriminant))/aq
    root1[np.isnan(root1)] = 0.
    root2[np.isnan(root2)] = 0.
    root1[root1 <= 0] = np.nan
    root2[root2 <= 0] = np.nan

    p_both = np.vstack([np.tile(root1, (3, 1)).T*rays,
                        np.tile(root2, (3, 1)).T*rays])
    l_both = p_both - tvec

    rdp = np.sum(unit_vector(p_both)*unit_vector(l_both), axis=1)
    rdp[np.isnan(rdp)] = 0.

    idx_arr = np.tile(range(npts), (2, 1))
    hits_inner_surf = (rdp > 0).reshape(2, npts)
    if np.any(np.sum(hits_inner_surf, axis=0) > 1):
        raise RuntimeError("something is wrong...")
    output[idx_arr[hits_inner_surf], :] = np.dot(
            l_both[hits_inner_surf.flatten(), :],
            rmat)
    return output

max_stars_repo_name: joelvbernier/tf-transforms

max_issues_repo_name: joelvbernier/tf-transforms

max_forks_repo_name: joelvbernier/tf-transforms


text:
import numpy as np
from numba import njit, float64, int64
from numba.types import Tuple
from ..structures import Trimesh

@njit(Tuple((float64[:,::1], int64[:,::1]))(float64[:,::1],int64[:,::1], int64[:,::1]), cache=True)
def _mid_point_subdivision_trimesh(vertices, faces, edges):

    vtx_dictionary = dict()
    support_set = set()
    vtx_dictionary[(-1.,-1.,-1.)] = -1
    support_set.add((-1.,-1.,-1.))
    n_vertices = vertices.shape[0]

    new_vertices = np.empty((edges.shape[0], 3), dtype=np.float64)
    new_faces    = np.empty((faces.shape[0]*4, 3), dtype=np.int64)

    j = 0
    for i in range(faces.shape[0]):

        new_vtx1 = (vertices[edges[i*3][0]] + vertices[edges[i*3][1]]) / 2
        new_vtx2 = (vertices[edges[i*3+1][0]] + vertices[edges[i*3+1][1]]) / 2
        new_vtx3 = (vertices[edges[i*3+2][0]] + vertices[edges[i*3+2][1]]) / 2

        index0 = 0
        index1 = 0
        index2 = 0

        v1 = (new_vtx1[0], new_vtx1[1], new_vtx1[2])
        v2 = (new_vtx2[0], new_vtx2[1], new_vtx2[2])
        v3 = (new_vtx3[0], new_vtx3[1], new_vtx3[2])

        if v1 not in support_set:

            new_vertices[j] = new_vtx1
            vtx_dictionary[v1] = j+n_vertices
            index0 = j+n_vertices
            j+=1

        else:
            index0 = vtx_dictionary[v1]

        if v2 not in support_set:

            new_vertices[j] = new_vtx2
            vtx_dictionary[v2] = j+n_vertices
            index1 = j+n_vertices
            j+=1

        else:
            index1 = vtx_dictionary[v2]

        if v3 not in support_set:

            new_vertices[j] = new_vtx3
            vtx_dictionary[v3] = j+n_vertices
            index2 = j+n_vertices
            j+=1

        else:
            index2 = vtx_dictionary[v3]

        new_faces[i*4] = np.array([faces[i][0], index0, index2])
        new_faces[i*4+1] = np.array([index0, index1, index2])
        new_faces[i*4+2] = np.array([index0, faces[i][1], index1])
        new_faces[i*4+3] = np.array([index2, index1, faces[i][2]])

    return np.concatenate((vertices, new_vertices[:j]), axis=0), new_faces

def mid_point_subdivision(mesh):

    v, f = _mid_point_subdivision_trimesh(mesh.vertices, mesh.faces, mesh.edges)

    return Trimesh(vertices=v, faces=f)

reason: import numpy,from numba


text:
import unittest
import networkx
from python_paris.paris import *

class TestParis(unittest.TestCase):

    def setUp(self):
        self.weighted_graph = nx.Graph()
        self.weighted_graph.add_nodes_from([0, 1, 2, 3, 4, 5])
        self.weighted_graph.add_weighted_edges_from([(0, 1, 1), (0, 2, 1), (1, 2, 1), (2, 3, 1), (3, 4, 1), (3, 5, 1),
                                                (4, 5, 1)])
        self.unweighted_graph = nx.Graph()
        self.unweighted_graph.add_nodes_from([0, 1, 2, 3, 4, 5])
        self.unweighted_graph.add_edges_from([(0, 1), (0, 2), (1, 2), (2, 3), (3, 4), (3, 5), (4, 5)])

    def test_paris(self):
        dendrogram_weighted = paris(self.weighted_graph)
        dendrogram_unweighted = paris(self.unweighted_graph)

        self.assertEqual(dendrogram_unweighted.any(), dendrogram_weighted.any())

max_stars_repo_name: Charpenb/python-paris

max_issues_repo_name: Charpenb/python-paris

max_forks_repo_name: Charpenb/python-paris

reason: import networkx


text:
import argparse
import os
import cv2
import numpy as np
import colorsys
from uma_mot.image_viewer import ImageViewer

class Track:
    def __init__(self, track_bbox, track_id):
        self.track_bbox = track_bbox
        self.track_id = track_id

def gather_sequence_info(sequence_dir, det_dir, result_dir):
    """Gather sequence information, such as image filenames, detections,
    groundtruth (if available).

    Parameters
    ----------
    sequence_dir : str
        Path to the MOTChallenge sequence directory.
    detection_file : str
        Path to the detection file.

    Returns
    -------
    Dict
        A dictionary of the following sequence information:

        * sequence_name: Name of the sequence
        * image_filenames: A dictionary that maps frame indices to image
          filenames.
        * detections: A numpy array of detections in MOTChallenge format.
        * groundtruth: A numpy array of ground truth in MOTChallenge format.
        * image_size: Image size (height, width).
        * min_frame_idx: Index of the first frame.
        * max_frame_idx: Index of the last frame.

    """
    image_dir = os.path.join(sequence_dir, "img1")
    image_filenames = {
        int(os.path.splitext(f)[0]): os.path.join(image_dir, f)
        for f in os.listdir(image_dir)}
    groundtruth_file = os.path.join(sequence_dir, "gt/gt.txt")

    # detection_dir = os.path.join(sequence_dir, "det", "det.txt")
    detection_dir = os.path.join(det_dir, '%s.txt' % os.path.basename(sequence_dir))
    # print(detection_dir)
    if os.path.exists(detection_dir):
        detections = np.loadtxt(detection_dir, delimiter=',')
    if os.path.exists(result_dir):
        results = np.loadtxt(result_dir, delimiter=',')
        # results = results[results[:, 6] == 1.0, :]
        # index = (results[:,-1]).tolist().index(min(results[:,-1]))
        # print(results[index])
    groundtruth = None
    if os.path.exists(groundtruth_file):
        groundtruth = np.loadtxt(groundtruth_file, delimiter=',')

    if len(image_filenames) > 0:
        image = cv2.imread(next(iter(image_filenames.values())),
                           cv2.IMREAD_GRAYSCALE)
        image_size = image.shape
    else:
        image_size = None

    if len(image_filenames) > 0:
        min_frame_idx = min(image_filenames.keys())
        max_frame_idx = max(image_filenames.keys())
    else:
        min_frame_idx = int(detections[:, 0].min())
        max_frame_idx = int(detections[:, 0].max())

    info_filename = os.path.join(sequence_dir, "seqinfo.ini")
    if os.path.exists(info_filename):
        with open(info_filename, "r") as f:
            line_splits = [l.split('=') for l in f.read().splitlines()[1:]]
            info_dict = dict(
                s for s in line_splits if isinstance(s, list) and len(s) == 2)

        update_ms = 1000 / int(info_dict["frameRate"])
    else:
        update_ms = None

    # feature_dim = detections.shape[1] - 10 if detections is not None else 0
    seq_info = {
        "sequence_name": os.path.basename(sequence_dir),
        "image_filenames": image_filenames,
        "detections": detections,
        "groundtruth": groundtruth,
        "image_size": image_size,
        "min_frame_idx": min_frame_idx,
        "max_frame_idx": max_frame_idx,
        "update_ms": update_ms,
        "results": results
    }
    return seq_info

def frame_callback(vis, frame_idx, results):  # 跟踪过程

    print("Processing frame %05d" % frame_idx)
    # if frame_idx>300:
    frame_indices = results[:, 0].astype(np.int)
    mask = frame_indices == frame_idx
    record_trks = []
    for row in results[mask]:
        # bbox, confidence = row[2:6], row[6]
        bbox, id = row[2:6], int(row[1])
        record_trks.append(Track(bbox, id))
    image = cv2.imread(
        seq_info["image_filenames"][frame_idx], cv2.IMREAD_COLOR)
    vis.set_image(image.copy())
    vis.draw_trackers(record_trks)

def create_unique_color_float(tag, hue_step=0.41):
    """Create a unique RGB color code for a given track id (tag).

    The color code is generated in HSV color space by moving along the
    hue angle and gradually changing the saturation.

    Parameters
    ----------
    tag : int
        The unique target identifying tag.
    hue_step : float
        Difference between two neighboring color codes in HSV space (more
        specifically, the distance in hue channel).

    Returns
    -------
    (float, float, float)
        RGB color code in range [0, 1]

    """
    h, v = (tag * hue_step) % 1, 1. - (int(tag * hue_step) % 4) / 5.
    r, g, b = colorsys.hsv_to_rgb(h, 1., v)
    return r, g, b

def create_unique_color_uchar(tag, hue_step=0.41):
    """Create a unique RGB color code for a given track id (tag).

    The color code is generated in HSV color space by moving along the
    hue angle and gradually changing the saturation.

    Parameters
    ----------
    tag : int
        The unique target identifying tag.
    hue_step : float
        Difference between two neighboring color codes in HSV space (more
        specifically, the distance in hue channel).

    Returns
    -------
    (int, int, int)
        RGB color code in range [0, 255]

    """
    r, g, b = create_unique_color_float(tag, hue_step)
    return int(255*r), int(255*g), int(255*b)

class Visualization(object):
    """
    This class shows tracking output in an OpenCV image viewer.
    """

    def __init__(self, seq_info, update_ms):
        image_shape = seq_info["image_size"][::-1]
        aspect_ratio = float(image_shape[1]) / image_shape[0]
        image_shape = 1024, int(aspect_ratio * 1024)
        if update_ms is None:
            update_ms = seq_info["update_ms"]
        self.viewer = ImageViewer(
            update_ms, image_shape, "Figure %s" % seq_info["sequence_name"])
        self.viewer.thickness = 2
        self.frame_idx = seq_info["min_frame_idx"]
        self.last_idx = seq_info["max_frame_idx"]
        self.results = seq_info["results"]
        # self.detections = detections
        # self.iou = iou

    def run(self, frame_callback):
        self.viewer.run(lambda: self._update_fun(frame_callback))

    def _update_fun(self, frame_callback):
        if self.frame_idx > self.last_idx:
            return False  # Terminate
        frame_callback(self, self.frame_idx, self.results)
        self.frame_idx += 1
        return True

    def set_image(self, image):
        self.viewer.image = image

    def draw_groundtruth(self, track_ids, boxes):
        self.viewer.thickness = 2
        for track_id, box in zip(track_ids, boxes):
            self.viewer.color = create_unique_color_uchar(track_id)
            self.viewer.rectangle(*box.astype(np.int), label=str(track_id))

    # def draw_detections(self, detections):
    #     self.viewer.thickness = 2
    #     self.viewer.color = 0, 0, 255
    #     for i, detection in enumerate(detections):
    #         self.viewer.rectangle(*detection.tlwh)

    def draw_trackers(self, tracks):
        self.viewer.thickness = 5
        for track in tracks:
            self.viewer.color = create_unique_color_uchar(track.track_id)
            self.viewer.rectangle(
                *track.track_bbox.astype(np.int), label=str(track.track_id))

def parse_args():

    parser = argparse.ArgumentParser(description="Siamese Tracking")
    parser.add_argument(
        "--sequence_dir", help="Path to the MOTChallenge sequence directory.",
        default=None, required=True)
    parser.add_argument(
        "--result_file", help="Tracking output in MOTChallenge file format.",
        default=None, required=True)
    parser.add_argument(
        "--detection_file", help="Path to custom detections (optional).",
        default=None)
    parser.add_argument(
        "--update_ms", help="Time between consecutive frames in milliseconds. "
        "Defaults to the frame_rate specified in seqinfo.ini, if available.",
        default=None)
    return parser.parse_args()

if __name__ == "__main__":
    args = parse_args()
    seq_info = gather_sequence_info(args.sequence_dir, args.detection_file, args.result_file)
    visualizer = Visualization(seq_info, update_ms=args.update_ms)
    visualizer.run(frame_callback)

max_stars_repo_name: yestinchen/UMA-MOT

max_issues_repo_name: yestinchen/UMA-MOT

max_forks_repo_name: yestinchen/UMA-MOT


text:
# -*- coding: utf-8 -*-
"""Generic stimulus generation functions."""

import warnings
import numpy as np
from scipy import signal
from threading import Timer

from ..io import read_wav
from .._sound_controllers import SoundPlayer
from .._utils import _wait_secs, string_types

def window_edges(sig, fs, dur=0.01, axis=-1, window='hann', edges='both'):
    """Window the edges of a signal (e.g., to prevent "pops")

    Parameters
    ----------
    sig : array-like
        The array to window.
    fs : float
        The sample rate.
    dur : float
        The duration to window on each edge. The default is 0.01 (10 ms).
    axis : int
        The axis to operate over.
    window : str
        The window to use. For a list of valid options, see
        ``scipy.signal.get_window()``, but can also be 'dpss'.
    edges : str
        Can be ``'leading'``, ``'trailing'``, or ``'both'`` (default).

    Returns
    -------
    windowed_sig : array-like
        The modified array (float64).
    """
    fs = float(fs)
    sig = np.array(sig, dtype=np.float64)  # this will make a copy
    sig_len = sig.shape[axis]
    win_len = int(dur * fs)
    if win_len > sig_len:
        raise RuntimeError('cannot create window of size {0} samples (dur={1})'
                           'for signal with length {2}'
                           ''.format(win_len, dur, sig_len))
    if window == 'dpss':
        from mne.time_frequency.multitaper import dpss_windows
        win = dpss_windows(2 * win_len + 1, 1, 1)[0][0][:win_len]
        win -= win[0]
        win /= win.max()
    else:
        win = signal.windows.get_window(window, 2 * win_len)[:win_len]
    valid_edges = ('leading', 'trailing', 'both')
    if edges not in valid_edges:
        raise ValueError('edges must be one of {0}, not "{1}"'
                         ''.format(valid_edges, edges))
    # now we can actually do the calculation
    flattop = np.ones(sig_len, dtype=np.float64)
    if edges in ('trailing', 'both'):  # eliminate trailing
        flattop[-win_len:] *= win[::-1]
    if edges in ('leading', 'both'):  # eliminate leading
        flattop[:win_len] *= win
    shape = np.ones_like(sig.shape)
    shape[axis] = sig.shape[axis]
    flattop.shape = shape
    sig *= flattop
    return sig

def rms(data, axis=-1, keepdims=False):
    """Calculate the RMS of a signal

    Parameters
    ----------
    data : array-like
        Data to operate on.
    axis : int | None
        Axis to operate over. None will operate over the flattened array.
    keepdims : bool
        Keep dimension operated over.
    """
    return np.sqrt(np.mean(data * data, axis=axis, keepdims=keepdims))

def play_sound(sound, fs=None, norm=True, wait=False, backend='auto'):
    """Play a sound

    Parameters
    ----------
    sound : array
        1D or 2D array of sound values.
    fs : int | None
        Sample rate. If None, the sample rate will be inferred from the sound
        file (if sound is array, it is assumed to be 44100).
    norm : bool
        If True, normalize sound to between -1 and +1.
    wait : bool
        If True, wait until the sound completes to return control.
    backend : str
        The backend to use.

    Returns
    -------
    snd : instance of SoundPlayer
        The object playing sound. Can use "stop" to stop playback. Note that
        the sound player will be cleared/deleted once the sound finishes
        playing.
    """
    sound = np.array(sound)
    fs_default = 44100
    if isinstance(sound, string_types):
        sound, fs_default = read_wav(sound)
    if fs is None:
        fs = fs_default
    if sound.ndim == 1:  # make it stereo
        sound = np.array((sound, sound))
    if sound.ndim != 2:
        raise ValueError('sound must be 1- or 2-dimensional')
    if norm:
        m = np.abs(sound).max() * 1.000001
        m = m if m != 0 else 1
        sound /= m
    if np.abs(sound).max() > 1.:
        warnings.warn('Sound exceeds +/-1, will clip')
    # For rtmixer it's possible this will fail on some configurations if
    # resampling isn't built in to the backend; when we hit this we can
    # try/except here and do the resampling ourselves.
    snd = SoundPlayer(sound, fs=fs, backend=backend)
    dur = sound.shape[1] / float(fs)
    snd.play()  # will clip as necessary
    del_wait = 0.5
    if wait:
        _wait_secs(dur)
    else:
        del_wait += dur
    if hasattr(snd, 'delete'):  # for backward compatibility
        Timer(del_wait, snd.delete).start()
    return snd

def add_pad(sounds, alignment='start'):
    """Add sounds of different lengths and channel counts together

    Parameters
    ----------
    sounds : list
        The sounds to add together. There is a maximum of two channels.
    alignment : str
        How to align the sounds. Either by ``'start'`` (add zeros at end),
        ``'center'`` (add zeros on both sides), or ``'end'`` (add zeros at
        start).

    Returns
    -------
    y : float
        The summed sounds. The number of channels will be equal to the maximum
        number of channels (by appending a copy), and the length will be equal
        to the maximum length (by appending zeros).

    Notes
    -----
        Even if the original sounds were all 0- or 1-dimensional, the output
        will be 2-dimensional (channels, samples).
    """
    if alignment not in ['start', 'center', 'end']:
        raise(ValueError("alignment must be either 'start', 'center', "
                         "or 'end'"))
    x = [np.atleast_2d(y) for y in sounds]
    if not np.all(y.ndim == 2 for y in x):
        raise ValueError('Sound data must have no more than 2 dimensions.')
    shapes = [y.shape for y in x]
    ch_max, len_max = np.max(shapes, axis=0)
    if ch_max > 2:
        raise ValueError('Only 1- and 2-channel sounds are supported.')
    for xi, (ch, length) in enumerate(shapes):
        if length < len_max:
            if alignment == 'start':
                n_pre = 0
                n_post = len_max - length
            elif alignment == 'center':
                n_pre = (len_max - length) // 2
                n_post = len_max - length - n_pre
            elif alignment == 'end':
                n_pre = len_max - length
                n_post = 0
            x[xi] = np.pad(x[xi], ((0, 0), (n_pre, n_post)), 'constant')
        if ch < ch_max:
            x[xi] = np.tile(x[xi], [ch_max, 1])
    return np.sum(x, 0)

max_stars_repo_name: drammock/expyfun

max_issues_repo_name: drammock/expyfun

max_forks_repo_name: drammock/expyfun


text:
import numpy as np
from basic_operations import choose_mode
from typing import List, AnyStr

X = List[List[float]]

def get_x_with_index(neighbors):
    return list(enumerate(neighbors))

def get_closest_neighbors(x: List, k: int):
    neighbors = get_x_with_index(x)
    neighbors.sort(key=lambda tup: tup[1])
    return neighbors[1:k + 1]

def calculate_density(closest_neighbors, k: int):
    values = list(map(lambda tup: tup[1], closest_neighbors))
    return 1 / ((1 / k) * sum(values))

def calculate_ard(i: int, x_string: AnyStr, k: int):
    x = np.array(choose_mode(x_string))
    closest_neighbors = get_closest_neighbors(x[i], k)
    density = calculate_density(closest_neighbors, k)
    ard_sum = 0
    for neighbor in closest_neighbors:
        (index, _) = neighbor
        closest_neighbors = get_closest_neighbors(x[index], k)
        ard_sum = ard_sum + calculate_density(closest_neighbors, k)

    return density / ((1 / k) * ard_sum)

# string_table = """O1 O2 O3 O4 O5 O6 O7 O8 O9 O10
# O1 0 8.55 0.43 1.25 1.14 3.73 2.72 1.63 1.68 1.28
# O2 8.55 0 8.23 8.13 8.49 6.84 8.23 8.28 8.13 7.66
# O3 0.43 8.23 0 1.09 1.10 3.55 2.68 1.50 1.52 1.05
# O4 1.25 8.13 1.09 0 1.23 3.21 2.17 1.29 1.33 0.56
# O5 1.14 8.49 1.10 1.23 0 3.20 2.68 1.56 1.50 1.28
# O6 3.73 6.84 3.55 3.21 3.20 0 2.98 2.66 2.50 3.00
# O7 2.72 8.23 2.68 2.17 2.68 2.98 0 2.28 2.30 2.31
# O8 1.63 8.28 1.50 1.29 1.56 2.66 2.28 0 0.25 1.46
# O9 1.68 8.13 1.52 1.33 1.50 2.50 2.30 0.25 0 1.44
# O10 1.28 7.66 1.05 0.56 1.28 3.00 2.31 1.46 1.44 0"""
#
# print(calculate_ard(i=1, x_string=string_table, k=2))

string_table2 = '5.7 6.0 6.2 6.3 6.4 6.6 6.7 6.9 7.0 7.4'

print(calculate_ard(i=9, x_string=string_table2, k=3))


text:
import argparse
import os
import numpy
import torch

from tqdm.auto import tqdm

from torchvision import transforms
from PIL import Image
import numpy as np

from pytorch_slim_cnn.slimnet import SlimNet

if __name__ == '__main__':
    torch.backends.cudnn.benchmark = True

    parser = argparse.ArgumentParser(description='seg')
    parser.add_argument('img_path', type=str)
    args = parser.parse_args()

    model = SlimNet.load_pretrained('./pytorch_slim_cnn/models/celeba_20.pth').cuda().eval()

    labels = np.array(['5_o_Clock_Shadow', 'Arched_Eyebrows', 'Attractive', 'Bags_Under_Eyes',
                       'Bald', 'Bangs', 'Big_Lips', 'Big_Nose', 'Black_Hair', 'Blond_Hair',
                       'Blurry', 'Brown_Hair', 'Bushy_Eyebrows', 'Chubby', 'Double_Chin',
                       'Eyeglasses', 'Goatee', 'Gray_Hair', 'Heavy_Makeup', 'High_Cheekbones',
                       'Male', 'Mouth_Slightly_Open', 'Mustache', 'Narrow_Eyes', 'No_Beard',
                       'Oval_Face', 'Pale_Skin', 'Pointy_Nose', 'Receding_Hairline',
                       'Rosy_Cheeks', 'Sideburns', 'Smiling', 'Straight_Hair', 'Wavy_Hair',
                       'Wearing_Earrings', 'Wearing_Hat', 'Wearing_Lipstick',
                       'Wearing_Necklace', 'Wearing_Necktie', 'Young'])

    transform = transforms.Compose([
        transforms.Resize((178, 218)),
        transforms.ToTensor(),
        transforms.Normalize((0.5, 0.5, 0.5), (0.5, 0.5, 0.5)),
    ])

    root = args.img_path
    count = 0
    total = 0
    for img_path in tqdm(os.listdir(root)):
        if img_path.lower().endswith(('.png', '.jpeg', '.jpg')):
            img = transform(Image.open(os.path.join(root, img_path))).cuda().unsqueeze(dim=0)
            logits = torch.sigmoid(model(img))
            predictions = (logits > 0.5).squeeze().cpu().numpy().astype(bool)
            attributes = labels[predictions]
            if 'Smiling' in attributes:
                count += 1
            total += 1

    print(args.img_path, count, total)

max_stars_repo_name: iviazovetskyi/rewriting

max_issues_repo_name: iviazovetskyi/rewriting

max_forks_repo_name: iviazovetskyi/rewriting


text:
"""
Plot Z500 comparisons between HIT and FIT experiments. These are
sea ice thickness perturbation experiments using WACCM4.

Notes
-----
    Author : Zachary Labe
    Date   : 13 August 2017
"""

### Import modules
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, addcyclic, shiftgrid
import nclcmaps as ncm
import datetime
import read_MonthlyOutput as MO
import calc_Utilities as UT

### Define directories
directorydata = '/surtsey/zlabe/simu/'
directoryfigure = '/home/zlabe/Desktop/'
#directoryfigure = '/home/zlabe/Documents/Research/SITperturb/Figures/'

### Define time
now = datetime.datetime.now()
currentmn = str(now.month)
currentdy = str(now.day)
currentyr = str(now.year)
currenttime = currentmn + '_' + currentdy + '_' + currentyr
titletime = currentmn + '/' + currentdy + '/' + currentyr
print('\n' '----Plotting Z - %s----' % titletime)

### Alott time series
year1 = 1900
year2 = 2000
years = np.arange(year1,year2+1,1)

### Call function for 500 mb height data
lat,lon,time,lev,Z500h = MO.readExperi(directorydata,'Z500','HIT','surface')
lat,lon,time,lev,Z500f = MO.readExperi(directorydata,'Z500','FIT','surface')

### Separate per periods (ON,DJ,FM)
Z500h_on = np.nanmean(Z500h[:,9:11,:,:],axis=1)
Z500f_on = np.nanmean(Z500f[:,9:11,:,:],axis=1)

Z500h_dj,Z500f_dj = UT.calcDecJan(Z500h,Z500f,lat,lon,'surface',1)

Z500h_fm = np.nanmean(Z500h[:,1:3,:,:],axis=1)
Z500f_fm = np.nanmean(Z500f[:,1:3,:,:],axis=1)

### Calculate period differenceds
diff_on = np.nanmean((Z500f_on-Z500h_on),axis=0)
diff_dj = np.nanmean((Z500f_dj-Z500h_dj),axis=0)
diff_fm = np.nanmean((Z500f_fm-Z500h_fm),axis=0)

### Calculate significance
stat_on,pvalue_on = UT.calc_indttest(Z500h_on,Z500f_on)
stat_dj,pvalue_dj = UT.calc_indttest(Z500h_dj,Z500f_dj)
stat_fm,pvalue_fm = UT.calc_indttest(Z500h_fm,Z500f_fm)

###########################################################################
###########################################################################
###########################################################################
### Plot 500 mb heights
plt.rc('text',usetex=True)
plt.rc('font',**{'family':'sans-serif','sans-serif':['Avant Garde']})

### Set limits for contours and colorbars
limit = np.arange(-60,60.1,1)
barlim = np.arange(-60,61,30)

### Begin plot
fig = plt.figure()
ax1 = plt.subplot(131)

m = Basemap(projection='ortho',lon_0=0,lat_0=89,resolution='l',
            area_thresh=10000.)

var, lons_cyclic = addcyclic(diff_on, lon)
var, lons_cyclic = shiftgrid(180., var, lons_cyclic, start=False)
lon2d, lat2d = np.meshgrid(lons_cyclic, lat)
x, y = m(lon2d, lat2d)

pvalue_onq,lons_cyclic = addcyclic(pvalue_on, lon)
pvalue_onq,lons_cyclic = shiftgrid(180.,pvalue_onq,lons_cyclic,start=False)

m.drawmapboundary(fill_color='white',color='dimgray',linewidth=0.7)
m.drawcoastlines(color='dimgray',linewidth=0.8)
parallels = np.arange(-90,90,45)
meridians = np.arange(-180,180,60)
#m.drawparallels(parallels,labels=[True,True,True,True],
#                linewidth=0.6,color='dimgray',fontsize=6)
#m.drawmeridians(meridians,labels=[True,True,True,True],
#                linewidth=0.6,color='dimgray',fontsize=6)
#m.drawlsmask(land_color='dimgray',ocean_color='mintcream')

cs = m.contourf(x,y,var,limit,extend='both')
cs1 = m.contourf(x,y,pvalue_onq,colors='None',hatches=['....'],
             linewidths=0.4)

ax1.annotate(r'\textbf{ON}',
            xy=(0, 0),xytext=(0.35,1.05),xycoords='axes fraction',
            fontsize=25,color='dimgrey',rotation=0)

cmap = ncm.cmap('nrl_sirkes')
cs.set_cmap(cmap)

###########################################################################

ax2 = plt.subplot(132)

m = Basemap(projection='ortho',lon_0=0,lat_0=89,resolution='l',
            area_thresh=10000.)

var, lons_cyclic = addcyclic(diff_dj, lon)
var, lons_cyclic = shiftgrid(180., var, lons_cyclic, start=False)
lon2d, lat2d = np.meshgrid(lons_cyclic, lat)
x, y = m(lon2d, lat2d)

pvalue_djq,lons_cyclic = addcyclic(pvalue_dj, lon)
pvalue_djq,lons_cyclic = shiftgrid(180.,pvalue_djq,lons_cyclic,start=False)

m.drawmapboundary(fill_color='white',color='dimgray',linewidth=0.7)
m.drawcoastlines(color='dimgray',linewidth=0.8)
parallels = np.arange(-90,90,45)
meridians = np.arange(-180,180,60)
#m.drawparallels(parallels,labels=[True,True,True,True],
#                linewidth=0.6,color='dimgray',fontsize=6)
#m.drawmeridians(meridians,labels=[True,True,True,True],
#                linewidth=0.6,color='dimgray',fontsize=6)
#m.drawlsmask(land_color='dimgray',ocean_color='mintcream')

cs = m.contourf(x,y,var,limit,extend='both')
cs1 = m.contourf(x,y,pvalue_djq,colors='None',hatches=['....'],
             linewidths=0.4)

ax2.annotate(r'\textbf{DJ}',
            xy=(0, 0),xytext=(0.35,1.05),xycoords='axes fraction',
            fontsize=25,color='dimgrey',rotation=0)

cmap = ncm.cmap('nrl_sirkes')
cs.set_cmap(cmap)

###########################################################################

ax3 = plt.subplot(133)

m = Basemap(projection='ortho',lon_0=0,lat_0=89,resolution='l',
            area_thresh=10000.)

var, lons_cyclic = addcyclic(diff_fm, lon)
var, lons_cyclic = shiftgrid(180., var, lons_cyclic, start=False)
lon2d, lat2d = np.meshgrid(lons_cyclic, lat)
x, y = m(lon2d, lat2d)

pvalue_fmq,lons_cyclic = addcyclic(pvalue_fm, lon)
pvalue_fmq,lons_cyclic = shiftgrid(180.,pvalue_fmq,lons_cyclic,start=False)

m.drawmapboundary(fill_color='white',color='dimgray',linewidth=0.7)
m.drawcoastlines(color='dimgray',linewidth=0.8)
parallels = np.arange(-90,90,45)
meridians = np.arange(-180,180,60)
#m.drawparallels(parallels,labels=[True,True,True,True],
#                linewidth=0.6,color='dimgray',fontsize=6)
#m.drawmeridians(meridians,labels=[True,True,True,True],
#                linewidth=0.6,color='dimgray',fontsize=6)
#m.drawlsmask(land_color='dimgray',ocean_color='mintcream')

cs = m.contourf(x,y,var,limit,extend='both')
cs1 = m.contourf(x,y,pvalue_fmq,colors='None',hatches=['....'],
             linewidths=0.4)

ax3.annotate(r'\textbf{FM}',
            xy=(0, 0),xytext=(0.35,1.05),xycoords='axes fraction',
            fontsize=25,color='dimgrey',rotation=0)

cmap = ncm.cmap('nrl_sirkes')
cs.set_cmap(cmap)

cbar_ax = fig.add_axes([0.312,0.23,0.4,0.03])
cbar = fig.colorbar(cs,cax=cbar_ax,orientation='horizontal',
                    extend='max',extendfrac=0.07,drawedges=False)
cbar.set_label(r'\textbf{m}',fontsize=11,color='dimgray')
cbar.set_ticks(barlim)
cbar.set_ticklabels(list(map(str,barlim)))
cbar.ax.tick_params(axis='x', size=.01)

plt.subplots_adjust(wspace=0.01)

plt.savefig(directoryfigure + 'Z500_diff_FIT-HIT.png',dpi=300)
print('Completed: Script done!')

max_stars_repo_name: zmlabe/ThicknessSensitivity

max_issues_repo_name: zmlabe/ThicknessSensitivity

max_forks_repo_name: zmlabe/ThicknessSensitivity


text:
import pyqtgraph as pg
import numpy as np
import sys
from PyQt5 import QtCore, QtGui
from .clickablePlotSettings_ui import Ui_LineSettingsDialog
from .PlotDataErrorItem import *
from ..packageSettings import config_options
from .fittingHelpers import *
from ..helperFuncs import *
from scipy.optimize import curve_fit as spcf

class signalBlocker(object):
    def __init__(self, toBlock):
        # print("Init. To block", toBlock)
        self.toBlock = toBlock
    def __call__(self, f):
        # print("called, f", f)
        def wrappedf(*args, **kwargs):
            args[0].ui.__getattribute__(self.toBlock).blockSignals(True)
            ret = f(*args, **kwargs)
            args[0].ui.__getattribute__(self.toBlock).blockSignals(False)
            return ret

        return wrappedf

class ClickablePlotWidget(pg.PlotWidget):
    # emitted when the window is closed
    sigPlotClosed = QtCore.pyqtSignal(object)
    sigMouseMoved = QtCore.pyqtSignal(object)

    # Emitted when an arrow key is pressed. Emits a
    # 2x1 np.array, first index is left/right (-1/+1) and
    # second is down/up (-1/+1, resp)
    sigArrowKeyPressed = QtCore.pyqtSignal(object)

    sigFitSettingsChange = QtCore.pyqtSignal(object)

    sigCrosshairsMoved = QtCore.pyqtSignal(object)

    def __init__(self, *args, **kwargs):
        super(ClickablePlotWidget, self).__init__(*args, **kwargs)

        self.opts = config_options.copy()
        self.updateOpts()

        self.selectedCurves = {}

        self.fitSettings = {"linearRegion": pg.LinearRegionItem(),
                            "button": None,
                            "curveToFit": None,
                            "fitFunc": None,
                            "p0": None,
                            "menu": None,
                            "fitCurve": pg.PlotDataItem(pen=pg.mkPen('k', width=3)),
                            "pText": pg.TextItem("Select a function", color=(0, 0, 0))}

        self.crosshairSettings = {"crosshairsVisible": False,
                                  "freeCrosshairs": False,
                                  "dataCrosshairs": False,
                                  "dataCrosshairsIdx": 0,
                                  "dataCrosshairsSource": 0}

        # Map these up here so I don't need to constantly figure out
        # how to add on to the contet menu of the view
        self.addContextMenu = self.plotItem.vb.menu.addMenu
        self.addContextAction = self.plotItem.vb.menu.addAction

        # set up a timer to tell when you've double
        # clicked a
        self.doubleClickTimer = QtCore.QTimer()
        self.doubleClickTimer.setInterval(250)
        self.doubleClickTimer.setSingleShot(True)

        self.sigFitSettingsChange.connect(self.updateFitSettings)

        self.scene().sigMouseMoved.connect(self.mousemove)
        self.scene().sigMouseClicked.connect(self.mouseClickedEvent)
        self.sigArrowKeyPressed.connect(self.updateFreeCrosshairs)
        self.sigArrowKeyPressed.connect(self.updateDataCrosshairs)

        self.doubleClickCurve = None

        p = pg.mkPen(color='r', width=3)
        self.crosshairs = [pg.InfiniteLine(pen=p), pg.InfiniteLine(angle=0, pen=p)]
        self.crosshairWindow = QtGui.QMessageBox(self)
        self.crosshairWindow.setWindowFlags(QtCore.Qt.WindowStaysOnTopHint)
        self.crosshairWindow.setWindowFlags(QtCore.Qt.Tool)
        self.crosshairWindow.setModal(False)
        self.plotItem.addItem(self.crosshairs[0])
        self.plotItem.addItem(self.crosshairs[1])
        self.sigCrosshairsMoved.connect(self.updateCrosshairWindow)
        # Hide these things
        self.removeFreeCrosshairs()

        self.plotItem.addItem(self.fitSettings["linearRegion"])
        self.fitSettings["button"] = QtGui.QToolButton()
        self.fitSettings["button"].setArrowType(QtCore.Qt.RightArrow)
        self.fitSettings["button"].clicked.connect(self.popupFitMenu)
        self.plotItem.vb.scene().addWidget(self.fitSettings["button"])

        self.fitSettings["linearRegion"].sigRegionChanged.connect(self.updateFitButtonPosition)
        self.fitSettings["linearRegion"].sigRegionChanged.connect(self.attemptFitting)
        self.plotItem.sigRangeChanged.connect(self.updateFitButtonPosition)
        self.plotItem.sigRangeChanged.connect(self.updateFitTextPosition)

        self.plotItem.addItem(self.fitSettings["pText"])
        self.fitSettings["pText"].hide()
        self.fitSettings["pText"].setFont(QtGui.QFont("", 24))
        self.fitSettings["pText"].setZValue(1001)

        # make it sit on top
        self.fitSettings["fitCurve"].setZValue(1000)
        # add it directly to the viewBox, NOT the plotItem, to avoid
        # adding it to the list of plotItem.curves, etc, since a lot
        # of my stuff depends heavily on plotItem.curves and I don't want
        # the fit curve to be seen with it.
        # Also, ignore the bounds so it can be ignored when calculating the autorange
        self.plotItem.vb.addItem(self.fitSettings["fitCurve"], ignoreBounds=True)

        self.plotItem.ctrl.logXCheck.toggled.connect(
            lambda: self.fitSettings["fitCurve"].setLogMode(
                self.plotItem.ctrl.logXCheck.isChecked(),
                self.plotItem.ctrl.logYCheck.isChecked()
            )
        )
        self.plotItem.ctrl.logYCheck.toggled.connect(
            lambda: self.fitSettings["fitCurve"].setLogMode(
                self.plotItem.ctrl.logXCheck.isChecked(),
                self.plotItem.ctrl.logYCheck.isChecked()
            )
        )
        self.fitSettings["fitCurve"].hide()

        self.removeFitRegion()

    def updateOpts(self, **opts):
        self.opts.update(opts)
        if self.opts["boundedAxes"]:
            self.plotItem.showAxis("top")

            self.plotItem.showAxis("right")
        for pos in ["top", "left", "right", "bottom"]:
            axis = self.plotItem.getAxis(pos)
            axis.setTickFont(QtGui.QFont("", self.opts["axisFontSize"]))
            axis.setLabel(**{"font-size":"{}pt".format(self.opts["axisLabelFontSize"])})
            axis.setPen(color=self.opts["foreground"],width=self.opts["axisThickness"])
            axis.setStyle(tickLength=self.opts["axisTickLength"])
            if axis.orientation in ["left", "right"]:
                # axis.setWidth(100)
                pass
                # axis.setStyle(tickTextOffset=20, tickTextWidth=20, autoExpandTextSpace=True)
            elif axis.orientation in ["bottom", "top"]:
                pass
            # axis.setMinimumWidth(300)
            # axis.setMinimumHeight(300)
            pg.ImageView

    def plot(self, *args, **kwargs):
        aargs, kkwargs = getPlotPens(pw=self, *args, **kwargs)
        kwargs.update(kkwargs)
        args = aargs
        p = self.plotItem.plot(*args, **kwargs)
        self.setItemClickable(p)
        return p

    def brazilPlot(self,*args, **kwargs):
        """
        Pass an nx3 as first arg, or kwargs data:
            Plot d[:,0]vs d[:,1]+d[:,2] :: d[:,1]-d[:,2]
            as d[:,1] is a data value, and d[:,2] is errors

        Pass two Nx2

        :param args:
        :param kwargs:
        :return:
        """
        data = kwargs.pop("data", None)
        data1 = kwargs.pop("data1", None)
        data2 = kwargs.pop("data2", None)
        args = list(args)

        if (None is data1 or None is data2) and None is data:
            # Assume it's in the *args... Figure out how
            # many numpy arrays there are.
            numnp = [ii for ii in args if type(ii) is np.ndarray]
            if len(numnp) == 0:
                raise RuntimeError("I don't know what to plot")
            if len(numnp) == 1:
                # Follow it down to below parsing
                if args[0].shape[1] in [3, ]:
                    data = args[0]
                    args.pop(0)
            if len(numnp) == 2:
                if args[0].shape[1] == 2 and args[1].shape[1] == 2:
                    data1 = args[0]
                    data2 = args[1]
                    args.pop(1)
                    args.pop(0)
                elif args[0].shape[1] == 1 and args[1].shape[1] == 2:
                    # args[0] is x,
                    # args[1] is [y1, y2]
                    data1 = np.column_stack((args[0], args[1][:, 0]))
                    data2 = np.column_stack((args[0], args[1][:, 1]))
                    args.pop(1)
                    args.pop(0)
            if len(numnp) == 3:
                if args[0].ndim == 1 and args[1].ndim == 1 and args[2].ndim == 1:
                    data = np.column_stack((args[0], args[1], args[2]))
                    args.pop(2)
                    args.pop(1)
                    args.pop(0)

        if data is not None:
            # either Nx3 or Nx4?
            # Is Nx3 [x, y1, y2], or [x, y, yerr]? need a flag?
            # I'm gonna assume the later cause that's how I'd use it...
            if data.shape[1] == 3:
                data1 = np.column_stack((data[:, 0], data[:, 1] + data[:, 2]))
                data2 = np.column_stack((data[:, 0], data[:, 1] - data[:, 2]))
            else:
                raise RuntimeError("Not sure how to parse this brazil input: {}".format(data.shape))

        aargs, kkwargs = getPlotPens(self, *args, **kwargs)
        kwargs.update(kkwargs)
        args = aargs

        c1 = pg.PlotDataItem(data1, **kwargs)
        c2 = pg.PlotDataItem(data2, **kwargs)
        c1.hide()
        c2.hide()
        # c1 = self.plot(data1, *args, **kwargs)
        # c2 = self.plot(data2, *args, **kwargs)
        kwargs.pop("symbolPen") #not used in fillbetween
        kwargs.pop("symbolBrush")  # not used in fillbetween
        kwargs.pop("symbol")  # not used in fillbetween
        brush = kwargs.pop("brush", None)
        if brush is None:
            color = kwargs.get("pen").color()

            alpha = kwargs.pop("alpha", None)
            if alpha is None:
                alpha = 0.5
            color.setAlphaF(alpha)

            brush = pg.mkBrush(color)
            kwargs["brush"] = brush
            # kwargs["pen"] = pg.mkPen(color)
        kwargs["pen"] = None
        fill = pg.FillBetweenItem(c1, c2, **kwargs)
        # p = self.plotItem.plot(*args, **kwargs)
        # self.setItemClickable(p)
        self.plotItem.addItem(fill)
        # self.setItemClickable(c1)
        # self.setItemClickable(c2)
        # self.plotItem.addItem(c1)
        # self.plotItem.addItem(c2)

        return fill

    def errorbars(self, *args, **kwargs):
        # """
        # create and add an errordataitem
        #
        # errorbars(<xdata>, <ydata>, <error>, *kwargs)
        # errorbars(<ydata>, <errors>, **kwargs) to assume
        #     x as an index
        # errorbars(x=<xdata>, y=<ydata>, errorbars=<errors>, **kwargs)
        # errorbars(Nx(2, 3), **kwargs) to unpack as
        #     x, y, yerr
        #
        # kwargs:
        # capWidth: widths of bar caps
        # """
        #
        # aargs, kkwargs = getPlotPens(self, *args, **kwargs)
        # kwargs.update(kkwargs)
        # args = aargs
        #
        #
        # if len(args)==2:
        #     kwargs.update(y=args[0])
        #     kwargs.update(errorbars=args[1])
        # elif len(args)==3:
        #     kwargs.update(x=args[0])
        #     kwargs.update(y=args[1])
        #     kwargs.update(errorbars=args[2])
        # elif len(args) == 1 and isinstance(args[0], np.ndarray):
        #     if args[0].shape[1] == 2:
        #         kwargs.update(y=args[0][:,0])
        #         kwargs.update(errorbars=args[0][:,1])
        #     elif args[0].shape[1] == 3:
        #         kwargs.update(x=args[0][:,0])
        #         kwargs.update(y=args[0][:,1])
        #         kwargs.update(errorbars=args[0][:,2])
        #     else:
        #         raise RuntimeError("I do not know how to parse this np.ndarray")
        # elif len(args)==1:
        #     raise RuntimeError("I do not know how to parse this argument", args)
        #
        # assert 'y' in kwargs
        # assert 'errorbars' in kwargs
        #
        # if 'capWidth' not in kwargs:
        #     kwargs['capWidth'] = config_options["errorBarCapWidth"]
        # kwargs['beam'] = kwargs.pop('capWidth')
        #
        # args, kwargs = getPlotPens(self, *args, **kwargs)
        erroritem = PlotDataErrorItem(pw = self, *args, **kwargs)

        self.addItem(erroritem)
        self.setItemClickable(erroritem)
        return erroritem

    def setItemClickable(self, item):
        """
        item: a data item which should be modified
        to add behavior for when it can be clicked
        and modified by the UI
        """

        item.curve.setClickable(True)
        item.sigClicked.connect(self._updateCurveClicked)
        item.opts['isEnabled'] = True
        item.opts['isSelected'] = False
        item.opts["plotWidget"] = self
        item.opts["xOffset"] = 0
        item.opts["yOffset"] = 0

    def updateLegendNames(self, curveItem):
        if curveItem.name() is None:
            return
        for sample, label in self.plotItem.legend.items:
            if sample.item is curveItem:
                label.setText(curveItem.opts["name"])
                # The legend will resize to fit new text,
                # but will not shrink back down after changing
                # to a shorter name.
                # For some reason. LegendItem.updateSize() doesn't
                # work properly.
                # For some reason, setting the geometry to 0 causes
                # it to resize properly.
                # Not sure what the miscommunication is from
                # PyQt/graph is.
                self.plotItem.legend.setGeometry(0, 0, 0, 0)
                break
        else:
            self.plotItem.legend.addItem(curveItem, name=curveItem.opts["name"])

    def _updateCurveClicked(self, plotDataItem):
        pen = plotDataItem.opts['pen']
        pen = pg.mkPen(pen)
        width = pen.width()
        # print("Curve clicked. Current width: {}".format(pen.width()))
        if self.doubleClickTimer.isActive():
            self.doubleClickTimer.stop()
            CurveItemSettings.getNewParameters(self, self.plotItem.curves, plotDataItem)
        elif plotDataItem.opts["isSelected"]:
            pen.setWidth(width - config_options["selectionThickness"])
            plotDataItem.setPen(pen)
            self.doubleClickTimer.start()
            plotDataItem.opts["isSelected"] = False
        else:
            pen.setWidth(width + config_options["selectionThickness"])
            plotDataItem.setPen(pen)
            plotDataItem.opts["isSelected"] = True
            self.doubleClickTimer.start()

    def closeEvent(self, *args, **kwargs):
        self.crosshairWindow.close()
        self.sigPlotClosed.emit(self)

    def mousemove(self, *args, **kwargs):
        self.sigMouseMoved.emit(self.plotItem.vb.mapSceneToView(args[0]))

    def keyPressEvent(self, ev):
        arr = [0, 0]
        if ev.key() == QtCore.Qt.Key_Left:
            arr[0] = -1
        elif ev.key() == QtCore.Qt.Key_Right:
            arr[0] = 1
        elif ev.key() == QtCore.Qt.Key_Up:
            arr[1] = 1
        elif ev.key() == QtCore.Qt.Key_Down:
            arr[1] = -1
        self.sigArrowKeyPressed.emit(np.array(arr))

        super(ClickablePlotWidget, self).keyPressEvent(ev)

    def updateFreeCrosshairs(self, newPos):
        if not self.crosshairSettings["freeCrosshairs"]:
            return
        drange = self.plotItem.vb.viewRange()
        xsp = np.abs(np.diff(drange[0])/100.)[0]
        ysp = np.abs(np.diff(drange[1])/100.)[0]

        change = np.array(newPos) * np.array([xsp, ysp])

        for i in [0, 1]:
            self.crosshairs[i].setPos(
                self.crosshairs[i].value() + change[i]
            )
        self.sigCrosshairsMoved.emit([i.value() for i in self.crosshairs])

    def updateDataCrosshairs(self, direction=np.array([0,0])):
        if not self.crosshairSettings["dataCrosshairs"]:
            return

        self.crosshairSettings["dataCrosshairsSource"] += direction[1]
        self.crosshairSettings["dataCrosshairsSource"] %= len(self.plotItem.curves)

        dataCurve = self.plotItem.curves[self.crosshairSettings["dataCrosshairsSource"]]
        xData, yData = dataCurve.getData()
        if xData is None or yData is None:
            print("it's nonE?")

        self.crosshairSettings["dataCrosshairsIdx"] += direction[0]
        self.crosshairSettings["dataCrosshairsIdx"] %= len(xData)

        self.crosshairs[0].setPos(
            xData[self.crosshairSettings["dataCrosshairsIdx"]]
        )
        self.crosshairs[1].setPos(
            yData[self.crosshairSettings["dataCrosshairsIdx"]]
        )
        self.sigCrosshairsMoved.emit([i.value() for i in self.crosshairs])

    def updateCrosshairWindow(self):
        x = self.crosshairs[0].value()
        y = self.crosshairs[1].value()
        if self.crosshairSettings["dataCrosshairs"]:
            curve = self.plotItem.curves[
                self.crosshairSettings["dataCrosshairsSource"]]
            idx =self.crosshairSettings["dataCrosshairsIdx"]
            xD, yD = curve.getData()
            name = curve.name()
            if name is None: name="Data"
            txt = "{}:\nx[{}]={}, y[{}]={}".format(name, idx, xD[idx], idx, yD[idx])
        else:
            txt = "x={}, y={}".format(x, y)
        self.crosshairWindow.setText(txt)

    def addFreeCrosshairs(self):
        self.crosshairSettings["crosshairsVisible"] = True
        self.crosshairSettings["freeCrosshairs"] = True
        drange = self.plotItem.vb.viewRange()
        self.crosshairs[0].setValue(np.mean(drange[0]))
        self.crosshairs[1].setValue(np.mean(drange[1]))

        # self.addItem(self.crosshairs[0])
        # self.addItem(self.crosshairs[1])
        self.crosshairs[0].show()
        self.crosshairs[1].show()

        self.crosshairWindow.show()
        self.updateCrosshairWindow()

    def removeFreeCrosshairs(self):
        self.crosshairSettings["crosshairsVisible"] = False
        self.crosshairSettings["freeCrosshairs"] = False
        self.crosshairWindow.hide()
        self.crosshairs[0].hide()
        self.crosshairs[1].hide()

    def addDataCrosshairs(self):
        self.crosshairSettings["crosshairsVisible"] = True
        self.crosshairSettings["dataCrosshairs"] = True

        self.updateDataCrosshairs()

        self.crosshairs[0].show()
        self.crosshairs[1].show()

        self.crosshairWindow.show()
        self.updateCrosshairWindow()

    def removeDataCrosshairs(self):
        self.crosshairSettings["crosshairsVisible"] = False
        self.crosshairSettings["dataCrosshairs"] = False
        self.crosshairWindow.hide()
        self.crosshairs[0].hide()
        self.crosshairs[1].hide()

    def addFitRegion(self):
        xrange = self.plotItem.vb.viewRange()[0]
        d = xrange[1]-xrange[0]
        self.fitSettings["linearRegion"].setRegion((xrange[0] + d*0.2, xrange[1] - d*0.2))
        self.updateFitButtonPosition()

        if self.fitSettings["curveToFit"] is None:
            self.fitSettings["curveToFit"] = self.plotItem.curves[0]

        self.fitSettings["linearRegion"].show()
        self.fitSettings["button"].show()
        self.fitSettings["fitCurve"].show()
        self.fitSettings["pText"].show()
        self.updateFitTextPosition()

    def updateFitButtonPosition(self):
        if not self.fitSettings["button"].isVisible(): return
        yrange = self.plotItem.vb.viewRange()[1]
        x = self.fitSettings["linearRegion"].getRegion()[1]
        p = self.plotItem.vb.mapViewToDevice(QtCore.QPointF(x, yrange[1]))

        self.fitSettings["button"].setGeometry(
            p.x(), p.y(), self.fitSettings["button"].geometry().width(), self.fitSettings["button"].geometry().height()
        )

    def updateFitTextPosition(self):
        if not self.fitSettings["pText"].isVisible(): return
        vrange = self.plotItem.vb.viewRange()
        dx = vrange[0][1] - vrange[0][0]
        dy = vrange[1][1] - vrange[1][0]

        x = np.sum(vrange[0])/2.
        self.fitSettings["pText"].setPos(vrange[0][0] + dx*0.05, vrange[1][1] - dy*0.05)

    def removeFitRegion(self):
        # self.plotItem.removeItem(self.fitSettings["linearRegion"])
        self.fitSettings["linearRegion"].hide()
        self.fitSettings["button"].hide()
        self.fitSettings["fitCurve"].hide()
        self.fitSettings["pText"].hide()

    def popupFitMenu(self):
        menu = self.getFitMenu()
        # Get the cursor position to poppup
        # the menu
        p = QtGui.QCursor.pos()
        menu.popup(p)

    def getFitMenu(self):
        menu = QtGui.QMenu()
        mData = menu.addMenu("Select Data...")
        mFit = menu.addMenu("Set Function...")

        for c in self.plotItem.curves:
            name = str(c.name())
            if name is None:
                name = str(c)
            a = mData.addAction(name)
            a.setCheckable(True)
            if c is self.fitSettings["curveToFit"]:
                a.setChecked(True)
            # Need the c=c to overcome a scoping thing in python
            # when you loop through lambdas like this
            a.triggered.connect(lambda x, c=c: self.sigFitSettingsChange.emit(c))

        funcs = list(getFuncs.keys())
        funcs.append("Other...")
        for n in funcs:
            a = mFit.addAction(n)
            a.setCheckable(True)
            if n == self.fitSettings["fitFunc"]:
                a.setChecked(True)
            a.triggered.connect(lambda x, n=n: self.sigFitSettingsChange.emit(n))

        self.fitSettings["menu"] = menu

        return menu

    def updateFitSettings(self, *args, **kwargs):
        if isinstance(args[0], pg.PlotDataItem):
            # Been told to change the data source
            self.fitSettings["curveToFit"] = args[0]
            self.fitSettings["p0"] = None
        elif isinstance(args[0], str):
            # been told to change the fit function
            if args[0] == "Other...": return
            self.fitSettings["fitFunc"] = args[0]
            self.fitSettings["p0"] = None
        self.attemptFitting()

    def attemptFitting(self):
        # get the full data
        try:
            x = self.fitSettings["curveToFit"].xData
            y = self.fitSettings["curveToFit"].yData
        except AttributeError:
            # no dataset has been selected, abort
            return
        if self.fitSettings["fitFunc"] is None: return

        # now need to slice it by the region of itnerest
        rg = self.fitSettings["linearRegion"].getRegion()
        xidx0 = np.where(x>rg[0])[0]
        xidx1 = np.where(x<rg[1])[0]
        xidx = list(set(xidx0) & set(xidx1))

        x = x[xidx]
        y = y[xidx]

        d = np.column_stack((x, y))
        d = d[d[:,0].argsort()]
        x = d[:,0]
        y = d[:,1]
        f = getFuncs[self.fitSettings["fitFunc"]]

        if self.fitSettings["p0"] is None:
            self.fitSettings["p0"] = f.guessParameters(x, y)
            # print("return from guess:", self.fitSettings["p0"])

        try:
            p, _ = spcf(f.__call__, x, y, p0 = self.fitSettings["p0"])
        except Exception as e:
            print("Error fitting function,",e)
            return
        self.fitSettings["p0"] = p

        # I really wanted to put this into a UI setting, but I can't
        # remember how things are laid out for this
        # This one will plot the full view range
        xfit = self.plotItem.vb.viewRange()[0]
        xfit = np.linspace(xfit[0], xfit[1])
        # This one will plot only within the linearregion
        ### xfit = np.linspace(x[0], x[-1], 250)

        y = f(xfit, *self.fitSettings["p0"])
        self.fitSettings["fitCurve"].setData(xfit, y)
        st = str(f).format(*p)

        self.fitSettings["pText"].setText(st, color=(0,0,0))

    def addFunctionLine(self):
        pass

    def addLegend(self, *args, **kwargs):
        if self.plotItem.legend is None:
            self.plotItem.addLegend()
        return self.plotItem.legend

    def mouseClickedEvent(self, ev):
        # Intercept if the crosshairs
        # are visible so we can move their position
        if self.crosshairSettings["crosshairsVisible"] and\
                        ev.button() == QtCore.Qt.LeftButton:

            ev.accept()
            p = self.plotItem.vb.mapSceneToView(ev.scenePos())

            # if free, set the crosshairs to the clicked position
            if self.crosshairSettings["freeCrosshairs"]:
                self.crosshairs[0].setPos(p.x())
                self.crosshairs[1].setPos(p.y())
            # otherwise, find the nearest datapoint and go to it.
            elif self.crosshairSettings["dataCrosshairs"]:
                xD, yD = self.plotItem.curves[
                    self.crosshairSettings["dataCrosshairsSource"]].getData()
                xD = np.nan_to_num(xD)
                yD = np.nan_to_num(yD)
                dist = (p.x()-xD)**2 + (p.y()-yD)**2
                self.crosshairSettings["dataCrosshairsIdx"] = np.argmin(dist)
                self.updateDataCrosshairs()
            self.updateCrosshairWindow()

class CurveItemSettings(QtGui.QDialog):
    # QListWidgets allows multiple selection of listItems
    # However, if you toggle the checkbox of a set of selected items,
    # it will de-select all but the checked one, which is frustrating
    # if you want to toggle it (i.e. to quickly see differences).
    # I put this in place to prevent this. If you check multiple
    # items, start a quick single-shot timer, which another
    # function will check to see if it's running to reselect those
    # items.
    _multiSelectHelper = {"timer": QtCore.QTimer(),
                          "selectionList": []}
    def __init__(self, *args, **kwargs):
        curveList = kwargs.pop('curves', None)
        clickedCurve = kwargs.pop('clickedCurve', None)
        super(CurveItemSettings, self).__init__(*args, **kwargs)
        if curveList is None:
            curveList = [pg.PlotDataItem()]
        self.listViewCurveItems = {}
        self.initUI(curveList, clickedCurve)

        self._multiSelectHelper["timer"].setInterval(100)
        self._multiSelectHelper["timer"].setSingleShot(True)

    def initUI(self, curveList, firstCurve):
        """

        :param curveItem: curveitim
        :type curveItem: pg.PlotDataItem
        :return:
        """
        self.ui = Ui_LineSettingsDialog()
        self.ui.setupUi(self)
        self.ui.sbLineWidth.setMinimum(0)
        self.ui.sbLineWidth.setSingleStep(1)
        self.ui.sbLineWidth.setOpts(int=True)

        self.ui.sbMarkerSize.setMinimum(0)
        self.ui.sbMarkerSize.setSingleStep(1)
        self.ui.sbMarkerSize.setOpts(int=True)
        self.originalSettings = {}

        initialListItem = None

        for curveItem in curveList:
            listItem = QtGui.QListWidgetItem(self.ui.lwCurves)
            listItem.setFlags(QtCore.Qt.ItemIsUserCheckable | QtCore.Qt.ItemIsSelectable |
                              QtCore.Qt.ItemIsEnabled | QtCore.Qt.ItemIsEditable)
            listItem.setCheckState(QtCore.Qt.Unchecked)
            listItem.setSelected(curveItem.opts["isSelected"])
            if curveItem is firstCurve:
                initialListItem = listItem

            name = curveItem.name()
            if name is None:
                name = curveItem
            name = str(name)
            listItem.setText(name)

            cs = QtCore.Qt.Checked if curveItem.opts["isEnabled"] else QtCore.Qt.Unchecked
            listItem.setCheckState(cs)

            self.ui.lwCurves.addItem(listItem)
            self.listViewCurveItems[curveItem] = listItem
            self.originalSettings[curveItem] = curveItem.opts.copy()

        self.ui.lwCurves.currentItemChanged.connect(self.handlecurrentItemChanged)
        self.ui.lwCurves.itemChanged.connect(self.handleitemChanged)
        self.ui.lwCurves.itemSelectionChanged.connect(self.handleitemSelectionChanged)

        self.ui.colLine.sigColorChanging.connect(self.updateLineColors)
        self.ui.cbLineStyle.currentIndexChanged.connect(self.updateLineStyle)
        self.ui.sbLineWidth.sigValueChanging.connect(self.updateLineWidth)

        self.ui.colMarker.sigColorChanging.connect(self.updateMarkerBrushColor)
        self.ui.cbMarkerStyle.currentIndexChanged.connect(self.updateMarkerStyle)
        self.ui.sbMarkerSize.sigValueChanging.connect(self.updateMarkerSize)

        self.show()
        # toggle it to set the interface to have the settings
        # of the clicked line
        self.handlecurrentItemChanged(initialListItem)

    def getCurveFromItem(self, item):
        for k, v in list(self.listViewCurveItems.items()):
            if v is item: return k

    def handlecurrentItemChanged(self, *args, **kwargs):
        item = args[0]
        curve = self.getCurveFromItem(item)
        pen = pg.mkPen(curve.opts["pen"])
        opts = curve.opts

        self.ui.colLine.blockSignals(True)
        self.ui.colLine.setColor(pen.color())
        self.ui.colLine.blockSignals(False)

        self.ui.cbLineStyle.blockSignals(True)
        self.ui.cbLineStyle.setCurrentIndex(pen.style())
        self.ui.cbLineStyle.blockSignals(False)

        self.ui.sbLineWidth.blockSignals(True)
        self.ui.sbLineWidth.setValue(pen.width() -
                                     config_options["selectionThickness"]*int(curve.opts["isSelected"]))
        self.ui.sbLineWidth.blockSignals(False)

        self.ui.cbMarkerStyle.blockSignals(True)
        self.ui.cbMarkerStyle.setCurrentIndex(
            self.ui.cbMarkerStyle.findText(str(opts["symbol"]))
        )
        self.ui.cbMarkerStyle.blockSignals(False)

        self.ui.sbMarkerSize.blockSignals(True)
        self.ui.sbMarkerSize.setValue(opts["symbolSize"])
        self.ui.sbMarkerSize.blockSignals(False)

        self.ui.colMarker.blockSignals(True)
        self.ui.colMarker.setColor(pg.mkBrush(opts["symbolBrush"]).color())
        self.ui.colMarker.blockSignals(False)

    @signalBlocker("lwCurves")
    def handleitemChanged(self, *args, **kwargs):
        clickedItem = args[0]
        if clickedItem not in self.ui.lwCurves.selectedItems():
            self.ui.lwCurves.setCurrentItem(clickedItem)
        for listItem in self.ui.lwCurves.selectedItems():
            curve = self.getCurveFromItem(listItem)
            if listItem is clickedItem:
                pass
            else:
                # Want to toggle other things if they're not in this list
                listItem.setCheckState(2*(not listItem.checkState()))
                listItem.setSelected(True)
            # if curve.opts["isEnabled"] ^ (listItem.checkState() == QtCore.Qt.Checked):
            if True:
                curve.opts["isEnabled"] = (listItem.checkState() == QtCore.Qt.Checked)
                if not curve.opts["isEnabled"]:
                    curve.curve.hide()
                    curve.scatter.hide()
                    if hasattr(curve, 'errorbars'):
                        curve.errorbars.hide()
                else:
                    curve.curve.show()
                    curve.scatter.show()
                    if hasattr(curve, 'errorbars'):
                        curve.errorbars.show()
        item = args[0]
        curve = self.getCurveFromItem(item)

        if curve.name() != str(item.text()):
            curve.opts["name"] = str(item.text())
            if curve.opts["plotWidget"].plotItem.legend is not None:
                curve.opts["plotWidget"].updateLegendNames(curve)

        self._multiSelectHelper["selectionList"] = self.ui.lwCurves.selectedItems()
        self._multiSelectHelper["timer"].start()

    @signalBlocker("lwCurves")
    def handleitemSelectionChanged(self, *args, **kwargs):
        if self._multiSelectHelper["timer"].isActive():
            for item in self._multiSelectHelper["selectionList"]:
                item.setSelected(True)
        # print("Timer active?", self._multiSelectHelper["timer"].isActive())
        # print("Selectionlist:", self._multiSelectHelper["selectionList"])

    def updateLineColors(self, colorButton):
        for listItem in self.ui.lwCurves.selectedItems():
            curve = self.getCurveFromItem(listItem)
            p = pg.mkPen(curve.opts["pen"])
            p.setColor(colorButton.color())
            curve.setPen(p)
        self.ui.colMarker.setColor(colorButton.color())
        self.updateMarkerBrushColor(colorButton)

    def updateLineStyle(self, newIdx):
        for listItem in self.ui.lwCurves.selectedItems():
            curve  = self.getCurveFromItem(listItem)
            p = pg.mkPen(curve.opts["pen"])
            p.setStyle(newIdx)
            curve.setPen(p)

    def updateLineWidth(self, sb):
        for listItem in self.ui.lwCurves.selectedItems():
            curve  = self.getCurveFromItem(listItem)
            p = pg.mkPen(curve.opts["pen"])
            p.setWidth(sb.value() +
                       config_options["selectionThickness"]*int(curve.opts["isSelected"]))
            curve.setPen(p)

    def updateMarkerStyle(self, newIdx):
        for listItem in self.ui.lwCurves.selectedItems():
            curve  = self.getCurveFromItem(listItem)
            if str(self.ui.cbMarkerStyle.currentText()) == "None":
                curve.setSymbol(None)
            else:
                curve.setSymbol(
                    str(self.ui.cbMarkerStyle.currentText())
                )

    def updateMarkerSize(self, sb):
        for listItem in self.ui.lwCurves.selectedItems():
            curve  = self.getCurveFromItem(listItem)
            curve.setSymbolSize(sb.value())

    def updateMarkerBrushColor(self, colorButton):
        for listItem in self.ui.lwCurves.selectedItems():
            curve  = self.getCurveFromItem(listItem)
            brush = pg.mkBrush(curve.opts["symbolBrush"])
            pen = pg.mkPen(curve.opts["symbolPen"])
            brush.setColor(colorButton.color())
            pen.setColor(colorButton.color())
            curve.setSymbolBrush(brush)
            curve.setSymbolPen(pen)

    @staticmethod
    def getNewParameters(parent, curves, clickedCurve = None):
        dialog = CurveItemSettings(curves=curves, clickedCurve = clickedCurve)
        parent.sigPlotClosed.connect(dialog.close)
        ok = dialog.exec_()
        if not ok:
            for curve in curves:
                # curve.opts = dialog.originalSettings[curve].copy()
                curve.opts.update(dialog.originalSettings[curve])
                curve.updateItems()
                if not curve.opts["isEnabled"]:
                    curve.curve.hide()
                    curve.scatter.hide()
                    if hasattr(curve, 'errorbars'):
                        curve.errorbars.hide()
                if curve.opts["plotWidget"].plotItem.legend is not None:
                    curve.opts["plotWidget"].updateLegendNames(curve)
        else:
            for curve in curves:
                if not curve.opts["isEnabled"]:
                    curve.opts["isSelected"] = False

if __name__=='__main__':
    ex = QtGui.QApplication([])
    wid = ClickablePlotWidget()
    wid.plotItem.addLegend()
    wid.plot([0, 1, 2, 3, 4], pen='g', width=5, name='green')
    p = wid.plot([1, -2, 4, 8, 2], pen='r', name='r')
    wid.plot([1]*5, pen='y')
    x = np.arange(5)
    y = np.random.normal(1, 1, 5)
    err = np.ones(5)

    err = wid.errorbars(np.column_stack((x, y, err)), pen='c', name='errors')

    wid.show()
    # sys.exit()
    sys.exit(ex.exec_())

max_stars_repo_name: SherwinGroup/Stele

max_issues_repo_name: SherwinGroup/Stele

max_forks_repo_name: SherwinGroup/Stele


text:
#!/usr/bin/env python
# -*- coding=utf-8 -*-

import numpy as np
import joblib
import argparse
import os

class GloveVec:
    def __init__(self, vocabulary, embed_dim=200):
        self.vocabulary = vocabulary
        vocab_size = len(self.vocabulary) + 1
        self.embed_mat = np.zeros((vocab_size, embed_dim))
        self.embed_idx = {}

        self._gen_word_idx_dict()
        self._pretrained_vec()
        self._gen_embed_mat()

    def _gen_word_idx_dict(self):
        self.idx2word = {}
        self.word2idx = {}
        idx = 1 # 0 is reserved
        for word in self.vocabulary:
            self.word2idx[word] = idx
            self.idx2word[idx] = word
            idx += 1

    def _pretrained_vec(self):
        with open('glove/glove.6B.200d.txt', 'r') as f:
            for line in f:
                values = line.split()
                word = values[0]
                coefs = np.asarray(values[1:], dtype='float32')
                self.embed_idx[word] = coefs
        assert len(self.embed_idx) == 400000

    def _gen_embed_mat(self):
        for word, idx in self.word2idx.items():
            embed_vec = self.embed_idx.get(word)
            if embed_vec is not None:
                self.embed_mat[idx] = embed_vec

def parse_args():
    parser = argparse.ArgumentParser()
    parser.add_argument('--dataset_root', type=str, required=True, help='Path to dataset root directory (flickr8k or COCO)')
    return parser.parse_args()

if __name__ == "__main__":
    args = parse_args()

    # Load vocabulary
    vocab = np.loadtxt(os.path.join(args.dataset_root, 'vocabulary.txt'), dtype=str)

    glove_vec = GloveVec(vocab)
    print('embedding_matrix:', glove_vec.embed_mat.shape)

    with open(os.path.join(args.dataset_root, 'word2idx.pkl'), 'wb') as f:
        joblib.dump(glove_vec.word2idx, f, compress=3)

    with open(os.path.join(args.dataset_root, 'idx2word.pkl'), 'wb') as f:
        joblib.dump(glove_vec.idx2word, f, compress=3)

max_stars_repo_name: satojkovic/DeepCaptioning


text:
import numpy as np
import pandas as pd

import galsim
from galsim import degrees, arcsec

from flux_conversion import npho_from_mag, integrated_throughput_in_band

'''
Given a table of galaxies, select those within a rectangular RA-Dec patch
spanning a given angular width on each side, centered at a given RA,Dec.
'''
def select_gals(center, halfwidth, gal_table):
    # Define a window around the image center
    # Be as generous with the RA boundary as possible
    ra_up = center.deproject(-halfwidth, -halfwidth).ra
    dec_up = center.deproject(0*arcsec, halfwidth).dec
    ra_down = center.ra - (ra_up - center.ra)
    dec_down = center.dec - (dec_up - center.dec)

    # Select the galaxies in this window
    ra_filter = gal_table.ra > ra_down/degrees
    ra_filter &= gal_table.ra < ra_up/degrees
    dec_filter = gal_table.dec > dec_down/degrees
    dec_filter &= gal_table.dec < dec_up/degrees
    return gal_table[ra_filter & dec_filter]

'''
Construct the light model for a single galaxy
'''
def draw_gal(gal, npho, dx, dy):
    '''
    First, construct the undistorted galaxy light profile:
    Model the disk and bulge separately, then add together.
    From the galsim docs: "The shear() method precisely preserves the area."
    So start with a radius (_eff_rad) of a circle with the same area as the final desired ellipse,
    and then apply an area-preserving shear to reach the desired ellipticities.
    '''
    # Bulge
    bulge_eff_rad = np.sqrt(gal.size_bulge_true * gal.size_minor_bulge_true)
    bulge = galsim.Sersic(n=gal.sersic_bulge, half_light_radius=bulge_eff_rad)
    bulge = bulge.shear(g1=gal.ellipticity_1_bulge_true, g2=gal.ellipticity_2_bulge_true)
    bulge *= gal.bulge_to_total_ratio_i / bulge.flux
    # Disk
    disk_eff_rad = np.sqrt(gal.size_disk_true * gal.size_minor_disk_true)
    disk = galsim.Sersic(n=gal.sersic_disk, half_light_radius=disk_eff_rad)
    disk = disk.shear(g1=gal.ellipticity_1_disk_true, g2=gal.ellipticity_2_disk_true)
    disk *= (1.-gal.bulge_to_total_ratio_i) / disk.flux
    # And finally sum
    lightmodel = disk + bulge

    '''
    Apply weak lensing distortions, now that all "true" (unlensed) quantities have been set.
    Need to first convert from "theoretical" shear and convergence to "reduced" shear and magnification.
    Using the same transformations as implemented here:
    http://galsim-developers.github.io/GalSim/_build/html/_modules/galsim/lensing_ps.html#theoryToObserved
    '''
    g1 = gal.shear_1/(1-gal.convergence)
    g2 = gal.shear_2/(1-gal.convergence)
    mu = 1/((1-gal.convergence)**2 - (gal.shear_1**2 + gal.shear_2**2))
    lightmodel = lightmodel.lens(g1, g2, mu)

    '''
    Overall flux normalization. To be compatible with the photon shooting procedure
    specified below, this should be the expected number of photons reaching the camera.
    '''
    lightmodel *= npho / lightmodel.flux

    lightmodel = lightmodel.shift(dx, dy)

    return lightmodel

'''
Combined galaxy light profile.
In expectation, i.e. without random exposure-dependent fluctuations.

Mean atmospheric throughput is modeled here, but the atmospheric PSF (since it will change from exposure to exposure) is not.
Similarly, the mean number of photons per galaxy is modeled here, but the actual number collected in any given exposure (which will be Poisson-distributed) is not.
Finally, mean galactic positions on the image are modeled here, but random shifts due to telescope dithering are not.
'''
def build_galaxy_composite(gals, image_center, params, verbose=False):
    # Truth information for all drawn galaxies
    true_pos = pd.DataFrame(np.zeros((len(gals.index), 2)), columns=['x','y'])
    npho = pd.DataFrame(data=np.zeros((len(gals.index), len(params.filters))), columns=params.filters)

    # Combined galaxy light profile
    # One for every filter
    gal_composite_drawer = [galsim.Sersic(n=1, scale_radius=1, flux=0.0) # Empty profile
                            for _ in range(len(params.filters))]

    integrated_total_throughput = {band: integrated_throughput_in_band(band, 'total')
                                    for band in params.filters}

    for g_iloc in range(len(gals.index)):
        if verbose:
            if (g_iloc+1)%500 == 0:
                print(g_iloc+1, '/', len(gals.index))

        gal = gals.iloc[g_iloc]

        '''
        The values returned by project() should have an orientation that is exactly
        compatible with shift(). However, project() returns Angle instances, while
        shift() expects bare numbers that are compatible with the pixel scale.
        Docs:
        http://galsim-developers.github.io/GalSim/_build/html/gsobject.html?highlight=drawimage#galsim.GSObject.shift
        http://galsim-developers.github.io/GalSim/_build/html/wcs.html?highlight=project#galsim.CelestialCoord.project
        '''
        gal_coord = galsim.CelestialCoord(gal.ra*degrees, gal.dec*degrees)
        dx, dy = image_center.project(gal_coord)
        # Since the pixel scale is given in arcsec, give dx and dy in arcsec as well
        dx, dy = dx/arcsec, dy/arcsec

        # The true position of this galaxy, in grid coordinates
        true_pos.x[g_iloc] = params.ngrid_x/2 - 1 + (1/params.scale)*dx
        true_pos.y[g_iloc] = params.ngrid_y/2 - 1 + (1/params.scale)*dy

        for band_idx, band in enumerate(params.filters):
            mag = gal[f'mag_{band}_noMW']
            # The expected number of photons collected from this galaxy in this band
            npho_gal = npho_from_mag(mag, integrated_total_throughput[band],
                                     params.exposure_time, params.area)
            npho[band][g_iloc] = npho_gal

            gal_drawer = draw_gal(gal, npho_gal, dx, dy)
            gal_composite_drawer[band_idx] += gal_drawer

    return gal_composite_drawer, true_pos, npho