text:
import os
import random
import json
import sys
import itertools
import functools
import warnings
import imageio
import shutil
from seq2seq import hlog

import numpy as np
import torch

from absl import app, flags
from tqdm import tqdm

from torch import optim
from torch.utils import data as torch_data
from torch.nn.utils.rnn import pad_sequence
from torch.utils.data import DataLoader
from torch.utils.data.distributed import DistributedSampler
from torch.nn.parallel import DistributedDataParallel

import options
from src import utils
from src.lex import FilterModel
from src.vqvae import VectorQuantizedVAE
from src.cvqvae import CVQVAE
from src.cvae import CVAE
from src.dae import DAE
from src import parallel
from src.datasets import get_data

import torchvision
from torchvision.utils import make_grid
from torchvision.utils import save_image

FLAGS = flags.FLAGS

flags.DEFINE_string("lex_path", default='',
                    help='A prelearned lexicon path to be used in text-image '
                         'vqvae models')

flags.DEFINE_string('vae_path', default='',
                    help='A pretrained vae path for conditional vae models.')

def evaluate_cvae(model,loader):
    val_recon_error = 0.0
    val_loss = 0.0
    cnt = 0
    for (cmd, img, _) in loader:
        img = img.to(utils.device)
        cmd = cmd.to(utils.device)
        cnt += img.shape[0]
        loss = model(img, cmd)
        x_tilde, *_ = model.predict(cmd)
        val_recon_error += (x_tilde-img).pow(2).sum().item()
        val_loss += loss.item() * img.shape[0]
        if cnt > 100:
            break
    return val_recon_error / cnt, val_loss / cnt

def train_cvae():
    assert FLAGS.vae_path != ""

    train, val, test = get_data(datatype=FLAGS.datatype,
                                dataroot=FLAGS.dataroot)
    vis_folder = utils.flags_to_path(FLAGS)
    os.makedirs(vis_folder, exist_ok=True)
    hlog.log("vis folder: %s" % vis_folder)

    if FLAGS.test:
        model = torch.load(FLAGS.model_path)
        model = model.to(utils.device)

        loader = DataLoader(train,
                            batch_size=FLAGS.n_batch,
                            shuffle=True,
                            collate_fn=train.collate,
                            num_workers=FLAGS.n_workers)

        test_loader = DataLoader(test,
                                 batch_size=36,
                                 shuffle=True,
                                 collate_fn=train.collate,
                                 num_workers=FLAGS.n_workers)
    else:
        if FLAGS.modeltype == "CVQVAE":
            model = CVQVAE(3,
                           FLAGS.h_dim,
                           FLAGS.n_latent,
                           train.vocab,
                           rnn_dim=FLAGS.rnn_dim,
                           n_codes=FLAGS.n_codes,
                           cc=FLAGS.commitment_cost,
                           decay=FLAGS.decay,
                           epsilon=FLAGS.epsilon,
                           beta=FLAGS.beta,
                           cmdproc=False,
                           size=train.size,
                           ).to(utils.device)
            model.vqvae.load_state_dict(torch.load(FLAGS.vae_path).state_dict())
            for p in model.vqvae.parameters():
                p.requires_grad = False
            if FLAGS.lex_path != "":
                model.set_src_copy(model.get_src_lexicon(FLAGS.lex_path))

            model = model.to(utils.device)

            loader = DataLoader(train,
                                batch_size=FLAGS.n_batch,
                                shuffle=True, collate_fn=train.collate,
                                num_workers=FLAGS.n_workers)
            test_loader = DataLoader(test,
                                     batch_size=32,
                                     shuffle=True,
                                     collate_fn=train.collate,
                                     num_workers=FLAGS.n_workers)

            model.eval()
            T = torchvision.transforms.ToPILImage(mode=train.color)
            with torch.no_grad():
                sample, _ = model.vqvae.sample(B=32)
                sample = sample.cpu().data * train.std[None, :, None, None] + train.mean[None, :, None, None]
                T(make_grid(sample.clip_(0, 1))).convert("RGB").save(os.path.join(vis_folder, "initial_samples.png"))
                generator = iter(test_loader)
                cmd, img, _ = next(generator)
                cmd = cmd.to(utils.device)
                img = img.to(utils.device)
                _, recon, *_ = model.vqvae(img, cmd)
                recon = recon.cpu().data * train.std[None, :, None, None] + train.mean[None, :, None, None]
                img = img.cpu().data * train.std[None, :, None, None] + train.mean[None, :, None, None]
                res = torch.cat((recon, img), 0).clip_(0, 1)
                T(make_grid(res)).convert("RGB").save(os.path.join(vis_folder, f"initial_recons.png"))

        elif FLAGS.modeltype == "CVAE":
            model = CVAE(3,
                         FLAGS.h_dim,
                         FLAGS.n_latent,
                         train.vocab,
                         beta=FLAGS.beta,
                         rnn_dim=FLAGS.rnn_dim,
                         size=train.size)

            model.vae.load_state_dict(torch.load(FLAGS.vae_path).state_dict())
            model = model.to(utils.device)
            for p in model.vae.parameters():
                p.requires_grad = False
        else:
            raise ValueError(f"Unknown model type {FLAGS.modeltype}")

        optimizer = optim.Adam(model.parameters(), lr=FLAGS.lr)
        train_loss = 0.
        cnt = 0.

        generator = iter(loader)
        model.train()
        model.vqvae.eval()
        for i in tqdm(range(FLAGS.n_iter)):
            try:
                cmd, img, _ = next(generator)
            except StopIteration:
                generator = iter(loader)
                cmd, img, _ = next(generator)

            cmd = cmd.to(utils.device)
            img = img.to(utils.device)

            optimizer.zero_grad()
            loss = model(img, cmd)
            loss.backward()
            optimizer.step()

            train_loss += loss * img.shape[0]
            cnt += img.shape[0]

            if i == 10 or (i+1) % 100 == 0:
                model.eval()
                with torch.no_grad():
                    hlog.log('%d iterations' % (i+1))
                    hlog.log('%.6f train loss' % (train_loss / cnt))
                    val_recon, val_loss = evaluate_cvae(model, test_loader)
                    hlog.log('val_recon_error: %.6f' % val_recon)
                    hlog.log('val_loss: %.6f' % val_loss)
                    T = torchvision.transforms.ToPILImage(mode=train.color)
                    test_iter = iter(test_loader)
                    for j in range(5):
                        cmd, img, _ = next(test_iter)
                        cmd = cmd.to(utils.device)
                        img = img.to(utils.device)
                        recon, pred_encodings, copy_probs = model.predict(cmd, sample=True, top_k=10)
                        pred_encodings = model.make_number_grid(pred_encodings.cpu())
                        copy_heatmap = model.make_copy_grid(copy_probs.cpu())
                        _, encodings, *_ = model.vqvae.encode(img)
                        encodings = model.make_number_grid(encodings.cpu())
                        encodings = torch.cat((copy_heatmap, pred_encodings, encodings), 0)
                        recon = recon.cpu().data * train.std[None, :, None, None] + train.mean[None, :, None, None]
                        img = img.cpu().data * train.std[None, :, None, None] + train.mean[None, :, None, None]
                        res = torch.cat((recon, img), 0).clip_(0, 1)
                        T(make_grid(encodings, nrow=encodings.shape[0]//3)).convert("RGB").save(os.path.join(vis_folder, f"{i}_{j}_encodings.png"))
                        T(make_grid(res, nrow=res.shape[0]//2)).convert("RGB").save(os.path.join(vis_folder, f"{i}_{j}.png"))
                        hlog.log("saved")
                model.train()
                model.vqvae.eval()
        torch.save(model, os.path.join(vis_folder, f"model.pt"))

    model.eval()
    with torch.no_grad():
        val_recon, val_loss = evaluate_cvae(model, test_loader)
        hlog.log('val_recon_error: %.6f' % val_recon)
        hlog.log('val_loss: %.6f' % val_loss)
        T = torchvision.transforms.ToPILImage(mode=train.color)
        test_iter = iter(test_loader)
        for j in range(5):
            cmd, img, _ = next(test_iter)
            cmd = cmd.to(utils.device)
            img = img.to(utils.device)
            hlog.log("sampling")
            recon, *_ = model.predict(cmd, top_k=10, sample=True)
            hlog.log("sampled")
            recon = recon.cpu().data * train.std[None, :, None, None] + train.mean[None, :, None, None]
            img = img.cpu().data * train.std[None, :, None, None] + train.mean[None, :, None, None]
            res = torch.cat((recon, img), 0).clip_(0, 1)
            hlog.log("saving samples")
            T(make_grid(res, nrow=res.shape[0]//2)).convert("RGB").save(os.path.join(vis_folder, f"eval_{j}.png"))
            hlog.log("saved")

def main(_):
    train_cvae()

if __name__ == "__main__":
    app.run(main)

max_stars_repo_name: ekinakyurek/lexsym

max_issues_repo_name: ekinakyurek/lexsym

max_forks_repo_name: ekinakyurek/lexsym


text:
from pyquil import Program
from pyquil.gates import X
from openfermion import QubitOperator, qubit_operator_sparse, IsingOperator
import numpy as np
import pytest

from .interfaces.estimator_test import EstimatorTests
from .interfaces.mock_objects import MockQuantumBackend, MockQuantumSimulator
from .estimator import (
    BasicEstimator,
    ExactEstimator,
    get_context_selection_circuit,
    get_context_selection_circuit_for_group,
)
from .measurement import ExpectationValues
from .circuit import Circuit

class TestEstimatorUtils:
    def test_get_context_selection_circuit_offdiagonal(self):
        term = ((0, "X"), (1, "Y"))
        circuit, ising_operator = get_context_selection_circuit(term)

        # Need to convert to QubitOperator in order to get matrix representation
        qubit_operator = QubitOperator()
        for ising_term in ising_operator.terms:
            qubit_operator += QubitOperator(
                ising_term, ising_operator.terms[ising_term]
            )

        target_unitary = qubit_operator_sparse(QubitOperator(term))
        transformed_unitary = (
            circuit.to_unitary().conj().T
            @ qubit_operator_sparse(qubit_operator)
            @ circuit.to_unitary()
        )

        assert np.allclose(target_unitary.todense(), transformed_unitary)

    def test_get_context_selection_circuit_diagonal(self):
        term = ((4, "Z"), (2, "Z"))
        circuit, ising_operator = get_context_selection_circuit(term)
        assert len(circuit.gates) == 0
        assert ising_operator == IsingOperator(term)

    def test_get_context_selection_circuit_for_group(self):
        group = QubitOperator(((0, "X"), (1, "Y"))) - 0.5 * QubitOperator(((1, "Y"),))
        circuit, ising_operator = get_context_selection_circuit_for_group(group)

        # Need to convert to QubitOperator in order to get matrix representation
        qubit_operator = QubitOperator()
        for ising_term in ising_operator.terms:
            qubit_operator += QubitOperator(
                ising_term, ising_operator.terms[ising_term]
            )

        target_unitary = qubit_operator_sparse(group)
        transformed_unitary = (
            circuit.to_unitary().conj().T
            @ qubit_operator_sparse(qubit_operator)
            @ circuit.to_unitary()
        )

        assert np.allclose(target_unitary.todense(), transformed_unitary)

class TestBasicEstimator(EstimatorTests):
    @pytest.fixture()
    def estimator(self, request):
        return BasicEstimator()

    @pytest.fixture()
    def target_operator(self, request):
        return QubitOperator("Z0")

    @pytest.fixture()
    def circuit(self, request):
        return Circuit(Program(X(0)))

    @pytest.fixture()
    def backend(self, request):
        return MockQuantumBackend(n_samples=20)

    @pytest.fixture()
    def n_samples(self, request):
        return 10

    def test_get_estimated_expectation_values(
        self, estimator, backend, circuit, target_operator, n_samples
    ):
        # When
        values = estimator.get_estimated_expectation_values(
            backend=backend,
            circuit=circuit,
            target_operator=target_operator,
            n_samples=n_samples,
        ).values
        value = values[0]
        # Then
        assert len(values) == 1
        assert value >= -1
        assert value <= 1

    def test_get_estimated_expectation_values_samples_from_backend(
        self,
        estimator,
        backend,
        circuit,
        target_operator,
    ):
        # Given
        # When
        values = estimator.get_estimated_expectation_values(
            backend=backend,
            circuit=circuit,
            target_operator=target_operator,
        ).values
        value = values[0]
        # Then
        assert len(values) == 1
        assert value >= -1
        assert value <= 1

    def test_n_samples_is_restored(self, estimator, backend, circuit, target_operator):
        # Given
        backend.n_samples = 5
        # When
        values = estimator.get_estimated_expectation_values(
            backend, circuit, target_operator, n_samples=10
        )
        # Then
        assert backend.n_samples == 5

    def test_get_estimated_expectation_values_with_constant(
        self, estimator, backend, circuit, n_samples
    ):
        # Given
        coefficient = -2
        constant_qubit_operator = QubitOperator((), coefficient) + QubitOperator(
            (0, "X")
        )

        # When
        values = estimator.get_estimated_expectation_values(
            backend=backend,
            circuit=circuit,
            target_operator=constant_qubit_operator,
            n_samples=n_samples,
        ).values
        value = values[1]
        # Then
        assert len(values) == 2
        assert coefficient == value

    def test_get_estimated_expectation_values_optimal_shot_allocation(
        self, estimator, backend, circuit, target_operator
    ):
        # TODO: After a deterministic testing backend is imlemented, this test
        # should be updated to actually check that shots are being correctly
        # allocated and the expectation values correctly estimated.

        # Given
        # When
        values = estimator.get_estimated_expectation_values(
            backend=backend,
            circuit=circuit,
            target_operator=target_operator,
            shot_allocation_strategy="optimal",
            n_total_samples=100,
        ).values
        value = values[0]
        # Then
        assert len(values) == 1
        assert value >= -1
        assert value <= 1

    def test_get_estimated_expectation_values_optimal_shot_allocation_with_prior(
        self, estimator, backend, circuit, target_operator
    ):
        # TODO: After a deterministic testing backend is imlemented, this test
        # should be updated to actually check that shots are being correctly
        # allocated and the expectation values correctly estimated.

        # Given
        # When
        estimator.prior_expectation_values = ExpectationValues(
            np.array([0 for _ in target_operator.terms])
        )
        values = estimator.get_estimated_expectation_values(
            backend=backend,
            circuit=circuit,
            target_operator=target_operator,
            shot_allocation_strategy="optimal",
            n_total_samples=100,
        ).values
        value = values[0]
        # Then
        assert len(values) == 1
        assert value >= -1
        assert value <= 1

    @pytest.mark.parametrize(
        "n_samples,n_total_samples,shot_allocation_strategy",
        [
            (None, 100, "uniform"),
            (100, None, "optimal"),
            (100, 100, "optimal"),
            (100, 100, "uniform"),
            (100, None, "foo"),
        ],
    )
    def test_get_estimated_expectation_values_invalid_options(
        self,
        estimator,
        backend,
        circuit,
        target_operator,
        n_samples,
        n_total_samples,
        shot_allocation_strategy,
    ):
        with pytest.raises(ValueError):
            estimator.get_estimated_expectation_values(
                backend=backend,
                circuit=circuit,
                target_operator=target_operator,
                shot_allocation_strategy=shot_allocation_strategy,
                n_total_samples=n_total_samples,
                n_samples=n_samples,
            )

class TestExactEstimator(EstimatorTests):
    @pytest.fixture()
    def estimator(self, request):
        return ExactEstimator()

    @pytest.fixture()
    def target_operator(self, request):
        return QubitOperator("Z0")

    @pytest.fixture()
    def circuit(self, request):
        return Circuit(Program(X(0)))

    @pytest.fixture()
    def backend(self, request):
        return MockQuantumSimulator()

    @pytest.fixture()
    def n_samples(self, request):
        return None

    def test_require_quantum_simulator(
        self, estimator, backend, circuit, target_operator
    ):
        backend = MockQuantumBackend()
        with pytest.raises(AttributeError):
            value = estimator.get_estimated_expectation_values(
                backend=backend,
                circuit=circuit,
                target_operator=target_operator,
            ).values

    def test_get_estimated_expectation_values(
        self, estimator, backend, circuit, target_operator
    ):
        # Given
        # When
        values = estimator.get_estimated_expectation_values(
            backend=backend,
            circuit=circuit,
            target_operator=target_operator,
            n_samples=None,
        ).values
        value = values[0]
        # Then
        assert len(values) == 1
        assert value >= -1
        assert value <= 1

ext: py

max_stars_repo_name: kottmanj/z-quantum-core

max_issues_repo_name: kottmanj/z-quantum-core

max_forks_repo_name: kottmanj/z-quantum-core

reason: import numpy


text:
import pandas as pd
import numpy as np
import scipy as s
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import preprocessing
plt.style.use('fivethirtyeight')
pd.set_option('display.max_columns', 500)
from pylab import cm
import warnings
warnings.filterwarnings('ignore')
from pywaffle import Waffle
from copy import deepcopy
from sklearn import feature_selection

def prepare():

    dataf = pd.read_csv(r"./books.csv")
    dataf.columns = ["Timestamp", "Format", "Beverage", "Music","Target", "Frequency","Connect","Variety"]
    dataf["VarCount"] = dataf.Variety.apply(lambda stri: len(stri.split(",")))   #Count of variety of books read by people
    format_count = dataf.groupby('Format')['Format'].count()
    return dataf, format_count

def plot1(dataf, format_count):

    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(22,8))
    sns.countplot(dataf.Format, palette="summer", ax=ax[0])
    ax[0].set_xlabel("Book format")
    ax[0].set_ylabel("Preferred choice of readers")
    sns.boxplot(dataf.Format, dataf.VarCount, palette='summer', ax=ax[1])
    ax[1].set_xlabel("Book format")
    ax[1].set_ylabel("No. of different genres read")
    plt.suptitle("Readership analysis", fontsize=20)
    ax[0].text(-0.05,format_count["Paperback"]+1, s=format_count["Paperback"])
    ax[0].text(0.95,format_count["Hard cover"]+1, s=format_count["Hard cover"])
    ax[0].text(1.95,format_count["Ebook"]+1, s=format_count["Ebook"])
    ax[0].text(2.95,format_count["Audiobook"]+1, s=format_count["Audiobook"])

def postprocess(dataf):

    dataf = dataf[(dataf.Format=="Paperback") | (dataf.Format=="Hard cover")]
    d = dataf[["Variety"]]                                                                      #Only the variety column
    d.Variety = d.Variety.apply(lambda s:s.strip(' ').split(","))                               #Picking up the varieties of genres
    genres = []
    for i in d.Variety:
        for j in i:
            genres.append(j.strip(' '))                                                         #Storing all the genre varieties in genres

    genres = pd.DataFrame(genres).groupby(0)[0].count()
    genres.drop(['In fiction- fantasy fiction to be exact. Adventure books. Enid Blyton( a little kiddish ik). I also would like to read all the Vedas one by one.','None other than those required for my coursework','Something else'],axis=0,inplace=True) #Dropping single ultra specific entrie(s)
    genres=genres.sort_values(0)                                                                #Sorting the genres dataframe

    dataf.Connect.replace({'Yes! no. Well maybe....':'May connect','Yes':'Want to connect','No':'Do not want to connect'},inplace=True)
    d = dataf.groupby('Connect')['Connect'].count()                                             #Grouping preferences of "Want to connect" attribute
    pphc = pd.DataFrame(dataf.Target)                                                           #Main column of interest
    pphc.Target.replace({"I wish I could read them without buying a lot of books":"Want to read without buying",
    "I want to build a library duh!":"Want to build a library"}, inplace=True)
    pphc = pphc.groupby('Target')['Target'].count()                                             #Grouping by interest to read books without buying
    dataf.Beverage.replace({'No drink necessary':'I do not drink but I know things','None':'I do not drink but I know things','No drink necessary':'I do not drink but I know things','I drink but not with books':'I do not drink but I know things'},inplace=True)
    bv =pd.DataFrame(dataf.groupby('Beverage')['Beverage'].count())                             #Grouping by interest of beverage consumption
    bv.drop(['There is no connection between books and beverage'],axis=0,inplace=True)
    bv.columns=['Count']

    return dataf, genres, pphc, d, bv

def plot2(dataf, genres):
    fig,ax=plt.subplots(nrows=1, ncols=2,figsize=(20,8))
    genres.plot(kind='barh', color="dodgerblue", ax=ax[0])
    #sns.catplot(genres, ax=ax[0])
    plt.suptitle('Diversity among readers', fontsize=20)
    ax[0].set_ylabel('Genre')
    ax[0].set_xlabel('Genre read by different people')
    sns.lineplot(x="VarCount", y="Frequency",data=dataf, ax=ax[1],color='royalblue')
    ax[1].set_xlabel('No. of different genres read by a reader')
    ax[1].set_ylabel('No. of books read by the reader in the last six months')
    #ax[0].axhline(genres.mean(), color='red')
    #ax[0].text(40,genres.mean()+0.5,"Mean", fontsize=18)

def plot3(pphc,d):
    plt.figure(figsize=(20,8))
    plt.pie(pphc, autopct='%2.1f%%',colors=['lavender','royalblue'], explode=[0.02,0.02], pctdistance=1.05,labels=pphc.index, labeldistance=1.1)
    plt.title("Interest among readers about a book cafe", fontsize=22)
    #draw circle
    centre_circle = plt.Circle((0,0),0.70,fc='white')
    fig = plt.gcf()
    fig.gca().add_artist(centre_circle)
    # Equal aspect ratio ensures that pie is drawn as a circle
    plt.axis('equal')
    plt.pie(d, radius=0.85, autopct='%1.1f%%', explode=[0.02,0.02,0.02], pctdistance=0.7, labels=d.index, labeldistance=0.55, colors=['mintcream','palegreen','mediumspringgreen'])
    plt.tight_layout()
    plt.show()

def plot4(bv):
    # To plot the waffle Chart
    fig = plt.figure(FigureClass = Waffle, rows = 3, values = bv. Count, labels = list(bv.index) , figsize=(20,8))
    plt.title('Beverage preference of readers',fontsize=20)

def plot5(dataf):
    da = dataf[dataf.Beverage=='I do not drink but I know things']
    x = pd.DataFrame(da.groupby('Connect')['Connect'].count())
    x.columns=['Count']
    fig = plt.figure(FigureClass = Waffle, rows = 1, values = x.Count, labels = list(x.index) , figsize=(20,4))
    plt.title('Interest to meet fellow book readers among people who do not drink while reading')

def plot6(dmod):
    dmod.Beverage.replace({'I do not drink but I know things':'No drink required','None':'No drink required','No drink necessary':'No drink required','I drink but not with books':'No drink required','Depends upon mood and time of day':'No drink required','There is no connection between books and beverage':'No drink required'},inplace=True)
    dmod.Connect.replace({'Yes! no. Well maybe....':'Maybe'},inplace=True)
    dmod.Music.replace({'Yes, like my life':'Yes to music','Nope':'No to music'},inplace=True)
    dk = dmod[["Beverage","Connect","Music"]].groupby(['Connect','Music','Beverage'])[["Beverage"]].count()
    dk.columns=["Count"]
    dk.sort_values(["Connect","Music"], ascending=[False,True], inplace=True)
    dk = pd.concat([dk.iloc[0:7,:],dk.iloc[13:20,:]])
    return dk

def modprocess(datax):

    datax = datax[datax.Beverage!='There is no connection between books and beverage']
    datax.Beverage.replace({'I do not drink but I know things':'No drink required','None':'No drink required','No drink necessary':'No drink required','I drink but not with books':'No drink required','Depends upon mood and time of day':'No drink required'},inplace=True)
    datax.Connect.replace({'Yes':'Yes to connect','No':'No to connect','Yes! no. Well maybe....':'Open to connect'}, inplace=True)
    datax.Music.replace({'Nope':'No to music','Yes, like my life':'Yes to music'},inplace=True)
    datax[list(pd.DataFrame(datax["Format"].unique())[0].sort_values())] = pd.get_dummies(datax.Format)
    datax[list(pd.DataFrame(datax["Beverage"].unique())[0].sort_values())] = pd.get_dummies(datax.Beverage)
    datax[list(pd.DataFrame(datax["Music"].unique())[0].sort_values())] = pd.get_dummies(datax.Music)
    datax.Target.replace({"I wish I could read them without buying a lot of books":1,"I want to build a library duh!":0}, inplace=True)
    datax[list(pd.DataFrame(datax["Connect"].unique())[0].sort_values())] = pd.get_dummies(datax.Connect)
    datax.drop(['Format','Beverage','Music','Connect','Timestamp','Variety'], axis=1, inplace=True)
    datax.drop(['Alcohol','No to connect','Audiobook','Ebook'],axis=1, inplace=True)
    return datax

def plot7(dframe):
    x = dframe.drop(['Target'],axis=1)
    y = dframe.Target
    from sklearn.feature_selection import SelectKBest, chi2, f_classif
    fs = SelectKBest(f_classif, k="all")
    fs.fit(x,y)
    sc = pd.concat([pd.DataFrame(x.columns),pd.DataFrame(fs.scores_)], axis=1)
    sc.columns = ['Feature','Score']
    sc.sort_values('Score', inplace=True, ascending=False)
    plt.figure(figsize=(20,10))
    sns.barplot(sc.Score, sc.Feature, color='seagreen')
    plt.title('Feature importance among people wanting to read books without buying')

def modeller(dmod):

    dmod['Fiction'] = dmod.Variety.apply(lambda x:'Fiction' in x).astype('int')
    dmod['Science fiction'] = dmod.Variety.apply(lambda x:'Science fiction' in x).astype('int')
    dmod['Humour'] = dmod.Variety.apply(lambda x:'Humour' in x).astype('int')
    dmod['Philosophy'] = dmod.Variety.apply(lambda x:'Philosophy' in x).astype('int')
    dmod['History'] = dmod.Variety.apply(lambda x:'History' in x).astype('int')
    dmod['Business'] = dmod.Variety.apply(lambda x:'Business' in x).astype('int')
    dmod['Literature'] = dmod.Variety.apply(lambda x:'Literature' in x).astype('int')
    dmod['Biography'] = dmod.Variety.apply(lambda x:'Biography' in x).astype('int')
    dmod['Travel'] = dmod.Variety.apply(lambda x:'Travel' in x).astype('int')
    dmod['Religion'] = dmod.Variety.apply(lambda x:'Religion' in x).astype('int')
    y = dmod.Target
    dmod.Beverage.replace({'Alcohol':'No drink required','Juice':'No drink required'}, inplace=True)
    dmod =pd.concat([dmod,pd.get_dummies(dmod.Format, drop_first=True),pd.get_dummies(dmod.Beverage, prefix='Beverage_', drop_first=True),
    pd.get_dummies(dmod.Music, drop_first=True), pd.get_dummies(dmod.Connect, prefix='Connect', drop_first=True)],axis=1)
    X=dmod.drop(['Timestamp','Variety','Format','Connect','Beverage','Music','Target'],axis=1)
    y=y.replace({'I wish I could read them without buying a lot of books':1 , 'I want to build a library duh!':0})
    return X,y

def get_palette(pal,n):
    a=[]
    cmap = cm.get_cmap(pal, n)    # PiYG

    for i in range(cmap.N):
        rgba = cmap(i)
        # rgb2hex accepts rgb or rgba
        a.append(matplotlib.colors.rgb2hex(rgba))
    return a

reason: import numpy,import scipy


text:
from . import BaseVideo

from pathlib import Path
import os
import numpy as np

# directory structure:
# /VSUMM
#   /database
#   /UserSummary
#       /user1
#       ...
#       /user5

class VSUMMVideo(BaseVideo):
    NUM_USERS = 5
    SUMMARY_NAME_PREFIX = "Frame"

    def __init__(self, name, dataset_dir="./VSUMM"):
        file_name = f"{name}.mpg"
        video_dir = Path(dataset_dir, "database")
        super(VSUMMVideo, self).__init__(file_name, video_dir)

    def ground_truth(self):
        data_dir = self.video_dir / f"../UserSummary/{self.name}"

        user_score = np.zeros((self.NUM_USERS, self.nframes), dtype=np.uint8)
        for user_index in range(0, self.NUM_USERS):
            user_dir = data_dir / f"user{user_index + 1}"
            for summary in user_dir.iterdir():
                summary_name = summary.stem
                if not summary_name.startswith(self.SUMMARY_NAME_PREFIX):
                    print("Not a summary file:", summary)
                    continue
                frame_index = int(summary_name[len(self.SUMMARY_NAME_PREFIX):])
                user_score[user_index, frame_index] = 1

        return user_score

max_stars_repo_name: jongwookyi/AC-SUM-GAN

max_issues_repo_name: jongwookyi/AC-SUM-GAN

max_forks_repo_name: jongwookyi/AC-SUM-GAN


text:
import pylab as pl
import numpy as np
#from sklearn import datasets, linear_model
import pandas as pd
import statsmodels.api as sm

# import the cleaned up dataset
df = pd.read_csv('../datasets/loanf.csv')

intrate = df['Interest.Rate']
loanamt = df['Loan.Amount']
fico = df['FICO.Score']

# prepping data for the normalize equantion
y = np.matrix(intrate).transpose()
x1 = np.matrix(fico).transpose()
x2 = np.matrix(loanamt).transpose()

#joins the two matrices
x = np.column_stack([x1, x2])

X = sm.add_constant(x)
model = sm.OLS(y, X)
f = model.fit()

print 'Coefficients: ', f.params[0:2]
print 'Intercept: ', f.params[2]
print 'P-Values: ', f.pvalues
print 'R-Squared: ', f.rsquared

max_stars_repo_name: blaiseyuri/LearnDataScience

reason: import numpy,import statsmodels


text:
import curses
from typing import Text
import numpy as np
import math

text = """
                 !#########       #
               !########!          ##!
            !########!               ###
         !##########                  ####
       ######### #####                ######
        !###!      !####!              ######
          !           #####            ######!
                        !####!         #######
                           #####       #######
                             !####!   #######!
                                ####!########
             ##                   ##########
           ,######!          !#############
         ,#### ########################!####!
       ,####'     ##################!'    #####
     ,####'            #######              !####!
    ####'                                      #####
    ~##                                          ##~ """

brightness = ".:-=+*#%@"
s = 30
#points = [[2*s*(x-0.5), 2*s*(y-0.5), 2*s*(z-0.5)]
#          for x in range(2) for y in range(2) for z in range(2)]

#points += [[x, -s, -s] for x in range(-s, s+1)]
#points += [[x, s, -s] for x in range(-s, s+1)]
#points += [[s, x, -s] for x in range(-s, s+1)]
#points += [[-s, -x, -s] for x in range(-s, s+1)]

#points += [[s, s, x] for x in range(-s, s+1)]
#points += [[s, -s, x] for x in range(-s, s+1)]
#points += [[-s, s, x] for x in range(-s, s+1)]
#points += [[-s, -s, x] for x in range(-s, s+1)]

#points += [[-s, x, s] for x in range(-s, s+1)]
#points += [[s, x, s] for x in range(-s, s+1)]
#points += [[x, s, s] for x in range(-s, s+1)]
#points += [[x, -s, s] for x in range(-s, s+1)]
screen = curses.initscr()
w, h = screen.getmaxyx()
curses.start_color()
curses.init_pair(1,curses.COLOR_YELLOW, curses.COLOR_RED)
curses.init_pair(2,curses.COLOR_YELLOW, curses.COLOR_YELLOW)

text = text.splitlines()

points=[]
for i, line in enumerate(text):
    for j, char in enumerate(line):

        if char != " ":
            points.append([2 * s * (i - len(text)/2)/len(text),
                          2 * s*(j - len(line)/2)/len(line), 0])
            points.append([2 * s * (i - len(text)/2)/len(text),
                          2 * s*(j - len(line)/2)/len(line), -1])
            points.append([2 * s * (i - len(text)/2)/len(text),
                          2 * s*(j - len(line)/2)/len(line), 1])
            points.append([2 * s * (i - len(text)/2)/len(text),
                          2 * s*(j - len(line)/2)/len(line), -2])
            points.append([2 * s * (i - len(text)/2)/len(text),
                          2 * s*(j - len(line)/2)/len(line), 2])

curses.curs_set(0)
curses.noecho()
curses.cbreak()
screen.keypad(True)
screen.clear()
height, width = screen.getmaxyx()
theta = math.pi/72

Rx = [[1.0, 0.0, 0.0],
      [0.0, math.cos(theta), -math.sin(theta)],
      [0.0, math.sin(theta), math.cos(theta)]]

Ry = [[math.cos(theta), 0,  math.sin(theta)],
      [0, 1, 0],
      [-math.sin(theta), 0, math.cos(theta)]]

rx = [[1.0, 0.0, 0.0],
      [0.0, math.cos(-theta), -math.sin(-theta)],
      [0.0, math.sin(-theta), math.cos(-theta)]]

ry = [[math.cos(-theta), 0,  math.sin(-theta)],
      [0, 1, 0],
      [-math.sin(-theta), 0, math.cos(-theta)]]

while True:
    screen.refresh()
    screen.clear()
    curses.napms(30)
    screen.bkgd(' ',curses.color_pair(1))
    for point in points:
        screen.addch(int(w/2 + w/2 / s * point[0]*np.interp(point[2], [-math.sqrt(3)*s, math.sqrt(3)*s], [0.5, 1])),
                     int(h/2 + h/2 / s * point[1]*np.interp(point[2], [-math.sqrt(3)*s, math.sqrt(3)*s], [0.5, 1])), brightness[int(np.interp(point[2], [-math.sqrt(3)*s, math.sqrt(3)*s], [0, len(brightness)-1]))],curses.color_pair(2))

    key = screen.getch()
    if key == curses.KEY_RIGHT:
        for i in range(len(points)):
            points[i] = [[sum(a*b for a, b in zip(X_row, Y_col))
                          for Y_col in zip(*Rx)] for X_row in [points[i]]][0]
    elif key == curses.KEY_UP:
        for i in range(len(points)):
            points[i] = [[sum(a*b for a, b in zip(X_row, Y_col))
                          for Y_col in zip(*Ry)] for X_row in [points[i]]][0]
    elif key == curses.KEY_LEFT:
        for i in range(len(points)):
            points[i] = [[sum(a*b for a, b in zip(X_row, Y_col))
                          for Y_col in zip(*rx)] for X_row in [points[i]]][0]
    elif key == curses.KEY_DOWN:
        for i in range(len(points)):
            points[i] = [[sum(a*b for a, b in zip(X_row, Y_col))
                          for Y_col in zip(*ry)] for X_row in [points[i]]][0]


text:
import pandas as pd
import numpy as np
import os

from bs4 import BeautifulSoup
from selenium import webdriver
from selenium.webdriver.chrome.options import Options
from sqlalchemy import create_engine, exc
from selenium.common.exceptions import NoSuchElementException
from selenium.webdriver.support.ui import WebDriverWait
from selenium.common.exceptions import TimeoutException
from selenium.webdriver.support.ui import WebDriverWait
from selenium.webdriver.support import expected_conditions as EC
from selenium.webdriver.common.by import By

class SpendingHistory():

    def __init__(self, email, password):
        self.email = email
        self.password = password

    def spending_history(self):

        try:

            options = Options()

            # --- heroku --- #
            options.binary_location = os.environ['GOOGLE_CHROME_BIN']

            options.add_argument('--disable-gpu')
            options.add_argument('--no-sandbox')
            options.add_argument('--headless')
            options.add_argument("--disable-extensions")

            # --- local testing --- #
            # driver = webdriver.Chrome(executable_path='chromedriver', chrome_options=options)

            # --- heroku --- #
            driver = webdriver.Chrome(executable_path=os.environ['CHROMEDRIVER_PATH'], chrome_options=options)

            print("[INFO] Opening browser...")

            driver.get("https://get.cbord.com/gwu/full/login.php")
            driver.find_element_by_id(
                "login_username_text").send_keys(self.email)
            driver.find_element_by_id(
                "login_password_text").send_keys(self.password)
            driver.find_element_by_name("submit").click()
            driver.get("https://get.cbord.com/gwu/full/history.php")

            soup = BeautifulSoup(driver.page_source, "lxml")
            driver.close()

            print("[INFO] Closing browser, found data...")

            try:

                table = soup.find('div', class_='history_table table-responsive')

                print("[INFO] Creating dataframe...")

                df = pd.read_html(table.prettify())[0]
                df['Amount ($ / Meals)'] = df['Amount ($ / Meals)'].replace({'\$': '', '\+ ': '', '- ': '-'}, regex=True)
                df['Date & Time'] = pd.to_datetime(df['Date & Time'])

                print('[INFO] Cleaning data...')
                df.dropna(inplace=True)
                df.columns = ['account', 'datetime', 'vendor', 'price']
                df['price'] = df['price'].astype(np.float64)
                df = df.reset_index(drop=True)
                df.index.names = ['id']
                df['email'] = self.email
                df.sort_values(by='datetime', inplace=True, ascending=True)

                # # --- local testing --- #
                # engine = create_engine("postgresql://suraj:password@localhost/gworld")

                # --- heroku --- #
                engine = create_engine("postgresql+psycopg2://gefimacpqdfbke:07654d71914a5e524f41fb4b9cb31e8f1e9ba4ba9f7676bcd4afadf1caa9e5c5@ec2-50-16-217-122.compute-1.amazonaws.com:5432/d6ncrb2kqoer99")

                print('[INFO] Storing in database...')

                for i in range(len(df)):
                    try:
                        df.iloc[i:i + 1].to_sql(name="history",
                                                if_exists='append', con=engine)
                    except exc.IntegrityError:
                        pass

                print("[INFO] Finished collecting data...")

            except AttributeError as e:
                print("Wrong username/password")
        except NoSuchElementException as e:
            print("No internet connection")


text:
from math import *
import numpy as np
import tqdm
from get_csection import cross_section
import matplotlib.pyplot as plt

def get_Halo_cs_relation(M200_min=1e13,M200_max=1e15):
    Mstar = 10**11.5
    Re = 3
    z1 = .3
    z2 = 1.5
    Halo_mass = np.logspace(np.log10(M200_min),np.log10(M200_max),200)
    cs_area = np.zeros(len(Halo_mass))

    for i in tqdm.tqdm(range(len(Halo_mass))):
        hmass = Halo_mass[i]
        csection = cross_section(M200=hmass)
        beta_caus = csection.get_loc()[1]
        cs_area[i] = pi*beta_caus**2

    plt.figure(figsize=(7,5))
    plt.loglog(Halo_mass,cs_area)
    plt.title(r'$M_*=$'+"%.3e"%Mstar+r'$M_{\odot},R_e =$'+str(Re)+r'$,z_{lens}=$'+str(z1)+r'$,z_{source}=$'+str(z2))
    plt.xlabel(r'$M_h(M_{\odot})$')
    plt.ylabel(r'$\sigma_{tot}(arcsec^2)$')
    plt.savefig("./plots/Halo_mass.png")

def get_Mstar_cs_relation(Mstar_min =10**11,Mstar_max=10**11.5):
    M200 = 1e13
    Re = 3
    z1 = .3
    z2 = 1.5
    Mstar_mass = np.logspace(np.log10(Mstar_min),np.log10(Mstar_max),200)
    cs_area = np.zeros(len(Mstar_mass))

    for i in tqdm.tqdm(range(len(Mstar_mass))):
        Mstar = Mstar_mass[i]
        csection = cross_section(Mstar=Mstar)
        beta_caus = csection.get_loc()[1]
        cs_area[i] = pi*beta_caus**2

    plt.figure(figsize=(7,5))
    plt.loglog(Mstar_mass,cs_area)
    plt.title(r'$M_h=$'+"%.3e"%M200+r'$M_{\odot},R_e =$'+str(Re)+r'$,z_{lens}=$'+str(z1)+r'$,z_{source}=$'+str(z2))
    plt.xlabel(r'$M_*(M_{\odot})$')
    plt.ylabel(r'$\sigma_{tot}(arcsec^2)$')
    plt.savefig("./plots/Stellar_mass.png")

def get_Re_cs_relation(Re_min=2,Re_max=5):
    M200 = 1e13
    Mstar = 10**11.5
    z1 = .3
    z2 = 1.5
    Reff = np.linspace(Re_min,Re_max,200)
    cs_area = np.zeros(len(Reff))

    for i in tqdm.tqdm(range(len(Reff))):
        Re = Reff[i]
        csection = cross_section(Re=Re)
        beta_caus = csection.get_loc()[1]
        cs_area[i] = pi*beta_caus**2

    plt.figure(figsize=(7,5))
    plt.loglog(Reff,cs_area)
    plt.title(r'$M_h=$'+"%.3e"%M200+r'$M_{\odot},M_*=$'+"%.3e"%Mstar+r'$M_{\odot},z_{lens}=$'+str(z1)+r'$,z_{source}=$'+str(z2))
    plt.xlabel(r'$R_e(kpc)$')
    plt.ylabel(r'$\sigma_{tot}(arcsec^2)$')
    plt.savefig("./plots/Re.png")

def get_zlens_cs_relation(zlens_min=.3,zlens_max=1):
    M200 = 1e13
    Mstar = 10**11.5
    Re = 3
    z2 = 1.5

    Zlens = np.logspace(np.log10(zlens_min),np.log10(zlens_max),200)
    cs_area = np.zeros(len(Zlens))

    for i in tqdm.tqdm(range(len(Zlens))):
        zlens = Zlens[i]
        csection = cross_section(z1=zlens)
        beta_caus = csection.get_loc()[1]
        cs_area[i] = pi*beta_caus**2

    plt.figure(figsize=(7,5))
    plt.loglog(Zlens,cs_area)
    plt.title(r'$M_h=$'+"%.3e"%M200+r'$M_{\odot},M_*=$'+"%.3e"%Mstar+r'$M_{\odot},R_e =$'+str(Re)+r'$,z_{source}=$'+str(z2))
    plt.xlabel(r'$z_{lens}$')
    plt.ylabel(r'$\sigma_{tot}(arcsec^2)$')
    plt.savefig("./plots/zlens.png")

def get_zsource_cs_relation(zsource_min=1.5,zsource_max=4):
    M200 = 1e13
    Mstar = 10**11.5
    Re = 3
    z1 = .3

    Z_source = np.logspace(np.log10(zsource_min),np.log10(zsource_max),200)
    cs_area = np.zeros(len(Z_source))

    for i in tqdm.tqdm(range(len(Z_source))):
        z2 = Z_source[i]
        csection = cross_section(z2=z2)
        beta_caus = csection.get_loc()[1]
        cs_area[i] = pi*beta_caus**2

    plt.figure(figsize=(7,5))
    plt.loglog(Z_source,cs_area)
    plt.title(r'$M_h=$'+"%.3e"%M200+r'$M_{\odot},M_*=$'+"%.3e"%Mstar+r'$M_{\odot},R_e =$'+str(Re)+r'$,z_{lens}=$'+str(z1))
    plt.xlabel(r'$z_{source}$')
    plt.ylabel(r'$\sigma_{tot}(arcsec^2)$')
    plt.savefig("./plots/zsource.png")

if __name__=="__main__":
    get_Halo_cs_relation()
    get_Mstar_cs_relation()
    get_Re_cs_relation()
    get_zlens_cs_relation()
    get_zsource_cs_relation()


text:
import scxx as sc
import scxx.preprocessing as pp
import scxx.plotting as pl
import sys
import numpy as np
import pandas as pd
from keras import optimizers

mode=str(sys.argv[1])
lambda_regulizer=float(sys.argv[2])
batch_input_size=int(sys.argv[3])
batch_input_size2=int(sys.argv[4])
test_input_file1="./data/atacseq/snareseq/12A/ExpressionMatrix_cDNA/"
test_input_file2="./data/atacseq/snareseq/12A/GeneActivity_chromatin/"
test_result_path="./results/atacseq/12A/scxx_%s/scmb_%s_%s_%s_%s/" % (str(sys.argv[1]).lower(),str(sys.argv[1]).lower(),str(sys.argv[2]),str(sys.argv[3]),str(sys.argv[4]))

ann = pd.read_csv("./data/atacseq/snareseq/12A/annotation.csv",sep=" ")
adata1 = pp.read_sc_data(test_input_file1,batch_name="ATAC",fmt="10x_mtx")
adata1.obs=adata1.obs.join(ann)
adata2 = pp.read_sc_data(test_input_file2,batch_name="RNA",fmt="10x_mtx")
adata2.obs=adata2.obs.join(ann)

# Run scxx. Here, rawdata is optional. It will be generated from datasets if not provided.
handle=sc.SCXX(
							datasets=[adata1, adata2],
							res_path=test_result_path,
							mode=mode,
							split_by="_batch",
							patience_es=100,
							patience_lr=50,
							mt_ratio=0.8,
							lambda_regulizer=lambda_regulizer,
							batch_input_size=batch_input_size,
							batch_input_size2=batch_input_size2,
							hidden_layers=[128,64,32,16],
							index_unique="-",
							#optimizer=optimizers.RMSprop(lr=0.01),
							model_file="model.h5"
							)
adata_transform=handle.fit_transform()
adata_transform.obs.to_csv("./data/atacseq/snareseq/12A/annotation_vipcca.csv")
#
import numpy as np
np.savetxt("./data/atacseq/snareseq/12A/genelist_python.txt",adata_transform.var_names.values,fmt="%s")

# pl.plotPrediction2(adata_transform.raw.X,adata_transform.X,result_path=test_result_path,lim=20)
# pl.run_embedding(adata_transform,path=test_result_path,method="umap")
# pl.plotEmbedding(adata_transform, path=test_result_path, method='umap', group_by="_batch",legend_loc="right margin")
#pl.plotEmbedding(adata_transform, path=test_result_path, method='umap', group_by="celltype",legend_loc="on data")
# pl.plotEmbedding(adata_transform, path=test_result_path, method='umap', group_by="louvain",legend_loc="right margin")


text:
#!/usr/bin/env python3

import cv2
import tensorflow as tf
import os
from imutils.video import WebcamVideoStream  # For more performant non-blocking multi-threaded OpenCV Web Camera Stream
from scipy.misc import imread
from lib.mtcnn import detect_face  # for MTCNN face detection
from flask import Flask, request, render_template
from werkzeug.utils import secure_filename
from waitress import serve
from utils import (
    load_model,
    get_face,
    get_faces_live,
    forward_pass,
    save_embedding,
    load_embeddings,
    identify_face,
    allowed_file,
    remove_file_extension,
    save_image
)

app = Flask(__name__)
app.secret_key = os.urandom(24)
APP_ROOT = os.path.dirname(os.path.abspath(__file__))
uploads_path = os.path.join(APP_ROOT, 'uploads')
embeddings_path = os.path.join(APP_ROOT, 'embeddings')
allowed_set = set(['png', 'jpg', 'jpeg'])  # allowed image formats for upload

@app.route('/upload', methods=['POST', 'GET'])
def get_image():
    """Gets an image file via POST request, feeds the image to the FaceNet model then saves both the original image
     and its resulting embedding from the FaceNet model in their designated folders.

        'uploads' folder: for image files
        'embeddings' folder: for embedding numpy files.
    """

    if request.method == 'POST':
        if 'file' not in request.files:
            return render_template(
                template_name_or_list="warning.html",
                status="No 'file' field in POST request!"
            )

        file = request.files['file']
        filename = file.filename

        if filename == "":
            return render_template(
                template_name_or_list="warning.html",
                status="No selected file!"
            )

        if file and allowed_file(filename=filename, allowed_set=allowed_set):
            filename = secure_filename(filename=filename)
            # Read image file as numpy array of RGB dimension
            img = imread(name=file, mode='RGB')

            # Detect and crop a 160 x 160 image containing a human face in the image file
            img = get_face(
                img=img,
                pnet=pnet,
                rnet=rnet,
                onet=onet,
                image_size=image_size
            )

            # If a human face is detected
            if img is not None:

                embedding = forward_pass(
                    img=img,
                    session=facenet_persistent_session,
                    images_placeholder=images_placeholder,
                    embeddings=embeddings,
                    phase_train_placeholder=phase_train_placeholder,
                    image_size=image_size
                )
                # Save cropped face image to 'uploads/' folder
                save_image(img=img, filename=filename, uploads_path=uploads_path)

                # Remove file extension from image filename for numpy file storage being based on image filename
                filename = remove_file_extension(filename=filename)

                # Save embedding to 'embeddings/' folder
                save_embedding(
                    embedding=embedding,
                    filename=filename,
                    embeddings_path=embeddings_path
                )

                return render_template(
                    template_name_or_list="upload_result.html",
                    status="Image uploaded and embedded successfully!"
                )

            else:
                return render_template(
                    template_name_or_list="upload_result.html",
                    status="Image upload was unsuccessful! No human face was detected!"
                )

    else:
        return render_template(
            template_name_or_list="warning.html",
            status="POST HTTP method required!"
        )

@app.route('/predictImage', methods=['POST', 'GET'])
def predict_image():
    """Gets an image file via POST request, feeds the image to the FaceNet model, the resulting embedding is then
    sent to be compared with the embeddings database. The image file is not stored.

    An html page is then rendered showing the prediction result.
    """
    if request.method == 'POST':
        if 'file' not in request.files:
            return render_template(
                template_name_or_list="warning.html",
                status="No 'file' field in POST request!"
            )

        file = request.files['file']
        filename = file.filename

        if filename == "":
            return render_template(
                template_name_or_list="warning.html",
                status="No selected file!"
            )

        if file and allowed_file(filename=filename, allowed_set=allowed_set):
            # Read image file as numpy array of RGB dimension
            img = imread(name=file, mode='RGB')

            # Detect and crop a 160 x 160 image containing a human face in the image file
            img = get_face(
                img=img,
                pnet=pnet,
                rnet=rnet,
                onet=onet,
                image_size=image_size
            )

            # If a human face is detected
            if img is not None:

                embedding = forward_pass(
                    img=img,
                    session=facenet_persistent_session,
                    images_placeholder=images_placeholder,
                    embeddings=embeddings,
                    phase_train_placeholder=phase_train_placeholder,
                    image_size=image_size
                )

                embedding_dict = load_embeddings()
                if embedding_dict:
                    # Compare euclidean distance between this embedding and the embeddings in 'embeddings/'
                    identity = identify_face(
                        embedding=embedding,
                        embedding_dict=embedding_dict
                    )

                    return render_template(
                        template_name_or_list='predict_result.html',
                        identity=identity
                    )

                else:
                    return render_template(
                        template_name_or_list='predict_result.html',
                        identity="No embedding files detected! Please upload image files for embedding!"
                    )

            else:
                return render_template(
                    template_name_or_list='predict_result.html',
                    identity="Operation was unsuccessful! No human face was detected!"
                )
    else:
        return render_template(
            template_name_or_list="warning.html",
            status="POST HTTP method required!"
        )

@app.route("/live", methods=['GET'])
def face_detect_live():
    """Detects faces in real-time via Web Camera."""

    embedding_dict = load_embeddings()
    if embedding_dict:
        try:
            # Start non-blocking multi-threaded OpenCV video stream
            cap = WebcamVideoStream(src=0).start()

            while True:
                frame_orig = cap.read()  # Read frame

                # Resize frame to half its size for faster computation
                frame = cv2.resize(src=frame_orig, dsize=(0, 0), fx=0.5, fy=0.5)

                # Convert the image from BGR color (which OpenCV uses) to RGB color
                frame = frame[:, :, ::-1]

                if cv2.waitKey(1) & 0xFF == ord('q'):
                    break

                if frame.size > 0:
                    faces, rects = get_faces_live(
                        img=frame,
                        pnet=pnet,
                        rnet=rnet,
                        onet=onet,
                        image_size=image_size
                    )

                    # If there are human faces detected
                    if faces:
                        for i in range(len(faces)):
                            face_img = faces[i]
                            rect = rects[i]

                            # Scale coordinates of face locations by the resize ratio
                            rect = [coordinate * 2 for coordinate in rect]

                            face_embedding = forward_pass(
                                img=face_img,
                                session=facenet_persistent_session,
                                images_placeholder=images_placeholder,
                                embeddings=embeddings,
                                phase_train_placeholder=phase_train_placeholder,
                                image_size=image_size
                            )

                            # Compare euclidean distance between this embedding and the embeddings in 'embeddings/'
                            identity = identify_face(
                                embedding=face_embedding,
                                embedding_dict=embedding_dict
                            )

                            cv2.rectangle(
                                img=frame_orig,
                                pt1=(rect[0], rect[1]),
                                pt2=(rect[2], rect[3]),
                                color=(0, 0, 255),  # Display Red (BGR)
                                thickness=2
                            )

                            W = int(rect[2] - rect[0]) // 2

                            cv2.putText(
                                img=frame_orig,
                                text=identity,
                                org=(rect[0] + W - (W // 2), rect[1]-7),
                                fontFace=cv2.FONT_HERSHEY_SIMPLEX,
                                fontScale=0.5,
                                color=(0, 0, 255),  # Display Red (BGR)
                                thickness=1,
                                lineType=cv2.LINE_AA
                            )

                        cv2.imshow(winname='Video', mat=frame_orig)
                    # Keep showing camera stream even if no human faces are detected
                    cv2.imshow(winname='Video', mat=frame_orig)
                else:
                    continue

            cap.stop()  # Stop multi-threaded Video Stream
            cv2.destroyAllWindows()

            return render_template(template_name_or_list='index.html')

        except Exception as e:
            print(e)

    else:
        return render_template(
            template_name_or_list="warning.html",
            status="No embedding files detected! Please upload image files for embedding!"
        )

@app.route("/")
def index_page():
    """Renders the 'index.html' page for manual image file uploads."""
    return render_template(template_name_or_list="index.html")

@app.route("/predict")
def predict_page():
    """Renders the 'predict.html' page for manual image file uploads for prediction."""
    return render_template(template_name_or_list="predict.html")

if __name__ == '__main__':
    """Server and FaceNet Tensorflow configuration."""

    # Load FaceNet model and configure placeholders for forward pass into the FaceNet model to calculate embeddings
    model_path = 'model/20170512-110547/20170512-110547.pb'
    facenet_model = load_model(model_path)
    config = tf.ConfigProto()
    config.gpu_options.allow_growth = True
    image_size = 160
    images_placeholder = tf.get_default_graph().get_tensor_by_name("input:0")
    embeddings = tf.get_default_graph().get_tensor_by_name("embeddings:0")
    phase_train_placeholder = tf.get_default_graph().get_tensor_by_name("phase_train:0")

    # Initiate persistent FaceNet model in memory
    facenet_persistent_session = tf.Session(graph=facenet_model, config=config)

    # Create Multi-Task Cascading Convolutional (MTCNN) neural networks for Face Detection
    pnet, rnet, onet = detect_face.create_mtcnn(sess=facenet_persistent_session, model_path=None)

    # Start flask application on waitress WSGI server
    serve(app=app, host='0.0.0.0', port=5000)

reason: from scipy


text:
# model.py
from mesa import Agent, Model
from mesa.time import RandomActivation
from mesa.space import MultiGrid
from mesa.datacollection import DataCollector

# import matplotlib.pyplot as plt
import random
import sys,time
import numpy as np
import operator

from draw_plane import DrawPlane
SAVE_FIG = False
RUN_LEARNING = True

START_POINT = 168
END_POINT = 190
TICK_RESOLUTION = 10
AI_CELL = 'Cpaaa'

PLANE_THRESHOLD = 3			#draw the cells in the [-2,+2] planes (5 total)
PLANE_DRAW = 7
CANVAS_WIDTH = 310
CANVAS_HEIGHT = 370
CANVAS_WIDTH_OFFSET = 140
CANVAS_HEIGHT_OFFSET = 180
CANVAS_DISPLAY_SCALE_FACTOR = 2		### zoom in/out the size of canvas
FRESH_TIME = 0.02
FRESH_PERIOD = 1

#Here used manually adjusted data
EMBRYO_LONG_AXIS_1 = [162, 194]
EMBRYO_LONG_AXIS_2 = [293, 351]
EMBRYO_SHORT_AXIS_1 = [178, 317]
EMBRYO_SHORT_AXIS_2 = [274, 232]

PLANE7_LONG_AXIS_1 = [170, 205]
PLANE7_LONG_AXIS_2 = [283, 333]
PLANE7_SHORT_AXIS_1 = [190, 303]
PLANE7_SHORT_AXIS_2 = [263, 237]

PLANE7_CENTER = [225, 275]
PLANE7_VEC0 = [1, 1]
PLANE7_VEC1 = [1, -1]
PLANE7_E0 = np.linalg.norm(np.array(PLANE7_LONG_AXIS_2) - np.array(PLANE7_CENTER))
PLANE7_E1 = np.linalg.norm(np.array(PLANE7_SHORT_AXIS_2) - np.array(PLANE7_CENTER))

EMBRYO_TOTAL_PLANE = 30
PLANE_RESOLUTION = 5

RADIUS_SCALE_FACTOR = 1.0     ####if manual adjustment needed, because radii are estimated.

#State cell list find from the Cell Neighbor Dertermination model
STATE_CELL_LIST = ['ABarpppap', 'ABarppppa', 'ABarppppp', 'Caaaa', 'ABprapapp', 'Epra', 'ABprapaaa', \
				'ABprapaap', 'Cpaap', 'ABprapapa', 'ABarppapp', 'Caaap', 'Eprp', 'ABarpppaa', 'Eplp', \
				'ABarppapa', 'Epla', 'ABarppaap']

AI_CELL_SPEED_PER_MIN = 2

AI_TARGET = [248, 299, 7.0*5]
AI_CELL_TARGET_TOLERANCE = 1
AI_CELL_BEGIN_REWARD = 20
class SeqRosModel(Model):
	def __init__(self):
		self.file_path = './nuclei/t%03d-nuclei'
		self.start_point = START_POINT
		self.end_point = END_POINT
		self.ai_cell = AI_CELL
		self.ai_cell_target = np.array(AI_TARGET)
		self.ai_cell_target_tolerance = AI_CELL_TARGET_TOLERANCE
		self.ticks = 0
		self.tick_resolution = TICK_RESOLUTION
		self.end_tick = (self.end_point - self.start_point) * self.tick_resolution
		self.stage_destination_point = self.start_point

		self.embryo_long_axis = np.linalg.norm(np.array(EMBRYO_LONG_AXIS_2) - np.array(EMBRYO_LONG_AXIS_1))
		self.embryo_short_axis = np.linalg.norm(np.array(EMBRYO_LONG_AXIS_2) - np.array(EMBRYO_LONG_AXIS_1))
		self.embryo_total_plane = EMBRYO_TOTAL_PLANE
		self.plane_resolution = PLANE_RESOLUTION
		self.embryo_volume = 4 / 3.0 * np.pi \
						* (0.5 * self.embryo_long_axis) \
						* (0.5 * self.embryo_short_axis) \
						* (0.5 * self.embryo_total_plane * self.plane_resolution)
		self.radius_scale_factor = RADIUS_SCALE_FACTOR

		self.current_cell_list = []
		self.dividing_cell_overall = []
		self.next_stage_destination_list = {}
		self.state_cell_list = STATE_CELL_LIST
		self.state_value_dict = {}
		self.schedule = RandomActivation(self)

		self.init_env()
		self.update_stage_destination()

		self.plane = DrawPlane(width=CANVAS_WIDTH,
							height=CANVAS_HEIGHT,
							w_offset = CANVAS_WIDTH_OFFSET,
							h_offset = CANVAS_HEIGHT_OFFSET,
							scale_factor = CANVAS_DISPLAY_SCALE_FACTOR)
		self.canvas = self.plane.canvas
		self.plane_draw = PLANE_DRAW
		self.draw(self.plane_draw)

		########## Reinforcement Learning related ##############
		self.ai_cell_speed = AI_CELL_SPEED_PER_MIN / float(TICK_RESOLUTION)
		self.n_observations = (len(self.state_cell_list) + 1) * 3
		self.actions = ['n','s','w','e','nw','ne','sw','se']
		self.n_actions = len(self.actions)

	def dist_point_ellipse(self, old_point, origin=PLANE7_CENTER, direc_vec0=PLANE7_VEC0, direc_vec1=PLANE7_VEC1):
		origin = np.array(origin)
		old_point = np.array(old_point)
		direc_vec0 = np.array(direc_vec0)
		direc_vec1 = np.array(direc_vec1)

		vec = old_point - origin
		vec_norm = np.linalg.norm(vec)
		direc_vec0_norm = np.linalg.norm(direc_vec0)
		direc_vec1_norm = np.linalg.norm(direc_vec1)
		cos_vec0 = vec.dot(direc_vec0) / vec_norm / direc_vec0_norm
		cos_vec1 = vec.dot(direc_vec1) / vec_norm / direc_vec1_norm
		x0 = vec_norm * cos_vec0
		x1 = vec_norm * cos_vec1
		new_point = np.array((x0,x1))

		new_point = np.abs(new_point)

		distance = -1.0
		if new_point[1] > 0:
			if new_point[0] > 0:
				###compute the unique root tbar of F(t) on (-e1*e1, +inf)
				t1 = - PLANE7_E1 * PLANE7_E1
				t2 = 1000000
				while True:
					ft_mid = (PLANE7_E0 * new_point[0] / (0.5 * (t1 + t2) + PLANE7_E0 ** 2)) ** 2 \
							+(PLANE7_E1 * new_point[1] / (0.5 * (t1 + t2) + PLANE7_E1 ** 2)) ** 2 - 1
					if ft_mid > 0:
						t1 = 0.5 * (t1 + t2)
					elif ft_mid < 0:
						t2 = 0.5 * (t1 + t2)
					else:
						break
					if np.abs(t1-t2)<0.001:
						break
				tbar = 0.5 * (t1 + t2)
				x0 = PLANE7_E0 * PLANE7_E0 * new_point[0] / (tbar + PLANE7_E0 * PLANE7_E0)
				x1 = PLANE7_E1 * PLANE7_E1 * new_point[1] / (tbar + PLANE7_E1 * PLANE7_E1)
				distance = np.sqrt((x0 - new_point[0]) ** 2 + (x1 - new_point[1]) ** 2)
			else:
				x0 = 0
				x1 = PLANE7_E1
				distance = np.abs(new_point[1] - PLANE7_E1)
		else:
			if new_point[0] < (PLANE7_E0 ** 2 - PLANE7_E1 ** 2) / PLANE7_E0:
				x0 = PLANE7_E0 * PLANE7_E0 * new_point[0] / (PLANE7_E0 ** 2 - PLANE7_E1 ** 2)
				x1 = PLANE7_E1 * np.sqrt(1 - (x0 / PLANE7_E0) ** 2)
				distance = np.sqrt((x0 - new_point[0]) ** 2 + x1 ** 2)
			else:
				x0 = PLANE7_E0
				x1 = 0
				distance = np.abs(new_point[0] - PLANE7_E0)

		return distance

	def draw(self,n_plane):
		self.canvas.delete("all")
		draw_range = np.arange(n_plane-PLANE_THRESHOLD, n_plane+PLANE_THRESHOLD+1, 1)
		draw_range = draw_range.tolist()
		draw_range.reverse()

		for n in draw_range:
			angle = np.pi *0.5 / (PLANE_THRESHOLD + 1) * np.abs(n - n_plane)
			for cell in self.schedule.agents:
				if cell.cell_name == self.ai_cell:
					type = 'AI'
					self.plane.draw_target(self.ai_cell_target[0:2], cell.diameter/2.0, self.ai_cell_target_tolerance)
				elif cell.cell_name in self.state_cell_list:
					type = 'STATE'
				else:
					type = 'NUMB'
				if round(cell.location[2]) == n:
					self.plane.draw_cell(center=cell.location[0:2],
										radius=cell.diameter/2.0*np.cos(angle),
										type=type)

		self.canvas.pack()
		self.canvas.update()
		time.sleep(FRESH_TIME)

	def get_radius(self, cell_name):
		if cell_name[0:2]=="AB":
			v=0.55*(0.5**(len(cell_name)-2))
		elif cell_name=="P1":
			v=0.45
		elif cell_name=="EMS":
			v=0.45*0.54
		elif cell_name=="P2":
			v=0.45*0.46
		elif cell_name[0:2]=="MS":
			v=0.45*0.54*0.5*(0.5**(len(cell_name)-2))
		elif cell_name=="E":
			v=0.45*0.54*0.5
		elif cell_name[0]=="E" and len(cell_name)>=2 and cell_name[1] != "M":
			v=0.45*0.54*0.5*(0.5**(len(cell_name)-1))
		elif cell_name[0]=="C":
			v=0.45*0.46*0.53*(0.5**(len(cell_name)-1))
		elif cell_name=="P3":
			v=0.45*0.46*0.47
		elif cell_name[0]=="D":
			v=0.45*0.46*0.47*0.52*(0.5**(len(cell_name)-1))
		elif cell_name=="P4":
			v=0.45*0.46*0.47*0.48
		else:
			print('ERROR!!!!! CELL NOT FOUND IN CALCULATING HER RADIUS!!!!')

		radius = pow(self.embryo_volume * v / (4 / 3.0 * np.pi), 1/3.0)
		radius = radius * self.radius_scale_factor

		return radius

	def get_cell_daughter(self, cell_name, cell_dict):
		daughter = []
		if cell_name == 'P0':
			daughter = ['AB', 'P1']
		elif cell_name == 'P1':
			daughter = ['EMS', 'P2']
		elif cell_name == 'P2':
			daughter = ['C', 'P3']
		elif cell_name == 'P3':
			daughter = ['D', 'P4']
		elif cell_name == 'P4':
			daughter = ['Z2', 'Z3']
		elif cell_name == 'EMS':
			daughter = ['MS', 'E']
		## standard name ###
		else:
			for cell in cell_dict.keys():
				if cell.startswith(cell_name) and len(cell) == len(cell_name) + 1:
					daughter.append(cell)
			daughter = sorted(daughter)
		if daughter == []:
			daughter = ['', '']
		return daughter

	def init_env(self):
		with open(self.file_path % self.start_point) as file:
			for line in file:
				line = line[:len(line)-1]
				vec = line.split(', ')
				id = int(vec[0])
				location = np.array((float(vec[5]), float(vec[6]), float(vec[7])))
				diameter = float(vec[8])
				cell_name = vec[9]
				if cell_name != '':
					self.current_cell_list.append(cell_name)
					a = CellAgent(id, self, cell_name, location, diameter)
					self.schedule.add(a)

	def set_cell_next_location(self,ai_action):
		for cell in self.schedule.agents:
			if cell.cell_name in self.next_stage_destination_list:
				cell.next_location = (self.next_stage_destination_list[cell.cell_name][0:3] - cell.location) \
								/ (self.tick_resolution - self.ticks % self.tick_resolution) + cell.location
				cell.diameter = self.next_stage_destination_list[cell.cell_name][3]
			else:
				### new cell born ###
				mother = cell.cell_name
				daughter = self.get_cell_daughter(cell.cell_name, self.next_stage_destination_list)
				if daughter[0] == '':
					print('ERROR!!!!! NO DAUGHTER FOUND!!!!!')
				cell.cell_name = daughter[0]
				cell.diameter = self.next_stage_destination_list[daughter[0]][3]
				cell.next_location = (self.next_stage_destination_list[daughter[0]][0:3] - cell.location) \
								/ (self.tick_resolution - self.ticks % self.tick_resolution) + cell.location
				new_id = len(self.schedule.agents) + 1
				new_diameter = self.next_stage_destination_list[daughter[1]][3]
				a = CellAgent(new_id, self, daughter[1], cell.location, new_diameter)
				self.schedule.add(a)
				a.next_location = (self.next_stage_destination_list[daughter[1]][0:3] - a.location) \
								/ (self.tick_resolution - self.ticks % self.tick_resolution) + a.location

				self.dividing_cell_overall.append(mother)

			if RUN_LEARNING and cell.cell_name == self.ai_cell:
				offset = self.get_ai_next_location_offset(ai_action)
				cell.next_location[0:2] = cell.location[0:2] + offset[0:2]

	def update_stage_destination(self):
		current_stage_destination_point = self.start_point + 1 + int(self.ticks / self.tick_resolution)
		if self.stage_destination_point == current_stage_destination_point:
			return
		else:
			self.stage_destination_point = current_stage_destination_point
			self.next_stage_destination_list.clear()
			with open(self.file_path % self.stage_destination_point) as file:
				for line in file:
					line = line[:len(line)-1]
					vec = line.split(', ')
					id = int(vec[0])
					loc_and_dia = np.array((float(vec[5]), float(vec[6]), float(vec[7]), float(vec[8])))
					cell_name = vec[9]
					if cell_name != '':
						self.next_stage_destination_list[cell_name] = loc_and_dia

	def render(self):
		if self.ticks % FRESH_PERIOD == 0:
			self.draw(self.plane_draw)

	def reset(self):
		self.ticks = 0
		self.stage_destination_point = self.start_point
		self.current_cell_list = []
		self.dividing_cell_overall = []
		self.next_stage_destination_list = {}
		self.state_value_dict = {}
		del self.schedule.agents[:]

		self.init_env()
		self.update_stage_destination()

		s = self.get_state()

		return s

	def get_state(self):
		s = np.round(self.state_value_dict[self.ai_cell], decimals=4).tolist()
		sorted_state_value_dict = sorted(self.state_value_dict.items(), key=operator.itemgetter(0))
		for item in sorted_state_value_dict:
			if item[0] != self.ai_cell:
				s = s + np.round(item[1], decimals=4).tolist()
		for i in range(0,len(s)):
			if i % 3 == 2:
				s[i] *= self.plane_resolution

		return s

	def get_reward(self):
		r = 0
		done = False
		ai_location = np.copy(self.state_value_dict[self.ai_cell])
		ai_location[2] = ai_location[2] * self.plane_resolution
		ai2bdr_dist = self.dist_point_ellipse(ai_location[0:2])
		ai_radius = self.get_radius(self.ai_cell)

		dist2target = np.linalg.norm(self.ai_cell_target - ai_location)
		if dist2target < AI_CELL_BEGIN_REWARD:
			r = (AI_CELL_BEGIN_REWARD -  dist2target)
			if dist2target < self.ai_cell_target_tolerance:
				r = 20

				return r, done

		#########boundary control rule: (1) dist>0.8*radius, ok  (2)0.5*r<dist<0.8*r, bad    (3)dist<0.5*r, dead
		if ai2bdr_dist > 0.8 * ai_radius:
			r += 0
		elif ai2bdr_dist > 0.5 * ai_radius and ai2bdr_dist <= 0.8 * ai_radius:
			r += (0.8 - float(ai2bdr_dist) / ai_radius) / (0.5 - 0.8)		## 0.5->-1, 0.8->0
		elif ai2bdr_dist < 0.5 * ai_radius:
			r = -10
			done = True
			print('hit boundary')
			return r, done

		######### pressure with other cells control rule ################
		######### (1) dist>0.6*r, ok  (2)0.3*r<dist<0.6*r, bad    (3)dist<0.3*r, dead
		for item in self.state_value_dict.keys():
			# print(r)
			if item != self.ai_cell:
				cell_location = np.copy(self.state_value_dict[item])
				cell_location[2] = cell_location[2] * self.plane_resolution
				dist = np.linalg.norm(cell_location - ai_location)
				cell_radius = self.get_radius(item)
				sum_radius = cell_radius + ai_radius
				dead_factor = 0.4
				ok_factor = 0.7
				if dist > ok_factor * sum_radius:
					r += 0
				elif dist > dead_factor * sum_radius and dist <= ok_factor * sum_radius:
					r += (ok_factor - float(dist) / sum_radius) / (dead_factor - ok_factor)		## 0.4->-1, 0.6->0

				elif dist < dead_factor * sum_radius:
					print('hit other cell:', item)
					r = -10
					done = True
					return r, done

		return r, done

	def get_ai_next_location_offset(self,a):
		offset_45 = self.ai_cell_speed * pow(2,0.5) / 2.0
		if a == 0:
			offset = np.array((0, - self.ai_cell_speed, 0))
		elif a == 1:
			offset = np.array((0, self.ai_cell_speed, 0))
		elif a == 2:
			offset = np.array((- self.ai_cell_speed, 0, 0))
		elif a == 3:
			offset = np.array((self.ai_cell_speed, 0, 0))

		elif a == 4:
			offset = np.array((- offset_45, -offset_45, 0))
		elif a == 5:
			offset = np.array((offset_45, -offset_45, 0))
		elif a == 6:
			offset = np.array((- offset_45, offset_45, 0))
		elif a == 7:
			offset = np.array((offset_45, offset_45, 0))
		elif a == 8:
			offset = np.array((0, 0, 0))
		return offset

	def step(self, a):
		done = False
		if self.ticks > 0 and self.ticks % self.tick_resolution == 0:
			self.update_stage_destination()
			if SAVE_FIG:
				name = 'screenshots_' + str(self.stage_destination_point-1) + '.ps'
				self.canvas.postscript(file=name, colormode='color')

		self.set_cell_next_location(ai_action=a)
		self.schedule.step()
		self.ticks += 1

		s_ = self.get_state()
		r, done = self.get_reward()
		if self.ticks == self.end_tick:
			done = True
		return s_, r, done

class CellAgent(Agent):
	def __init__(self, unique_id, model, name, location, diameter):
		super().__init__(unique_id, model)
		self.cell_name = name
		self.location = location
		self.diameter = diameter
		self.min_dist_list = []
		self.next_location = None

		self.set_state()

	def set_state(self):
		if self.cell_name == self.model.ai_cell:
			self.model.state_value_dict[self.cell_name] = self.location
		elif self.cell_name in self.model.state_cell_list:
			self.model.state_value_dict[self.cell_name] = self.location

	def move(self):
		self.location = self.next_location
		self.next_location = None

	def step(self):
		self.move()
		self.set_state()

if __name__ == '__main__':
	model = SeqRosModel()
	for i_episode in range(10):
		s = model.reset()
		counter = 0
		r_overall = 0
		while True:
			model.render()
			a = 2
			s_, r, done = model.step(a)
			counter += 1
			r_overall += r
			if done:
				break
		print('Episode:', i_episode, 'Done in', counter, 'steps. Reward:',r_overall)


text:
"""

mweigert@mpi-cbg.de

"""
from __future__ import absolute_import
from __future__ import print_function
import numpy as np
from biobeam.simlsm.simlsm import SimLSM_Base
from six.moves import range

class SimLSM_DSLM(SimLSM_Base):

    def _prepare_u0_illum(self, zfoc):
        self.u0_illum = self._bpm_illum.u0_beam(NA = self.NA_illum, zfoc = zfoc)

    def propagate_illum_single(self,cz = 0, **bpm_kwargs):
        bpm_kwargs.update({"return_comp":"intens"})
        offset = int(cz/self._bpm_illum.dy)
        assert abs(offset)<= self.u0_illum.shape[0]//2

        print("offset: ",offset)
        u0 = np.roll(self.u0_illum, offset ,axis=0)
        u = self._bpm_illum.propagate(u0,**bpm_kwargs)
        return self._trans_illum(u, inv = True)

    def propagate_illum(self,cz = 0,  dx_parallel = None,**bpm_kwargs):
        bpm = self._bpm_illum
        bpm_kwargs.update({"return_comp":"intens"})
        offset = int(cz/bpm.dy)
        assert abs(offset)<= self.u0_illum.shape[0]//2
        print("offset",  offset)

        u0_base = np.roll(self.u0_illum, offset ,axis=0)

        #prepare the parallel scheme
        max_NA = self.NA_illum if np.isscalar(self.NA_illum) else max(self.NA_illum)

        if dx_parallel is None:
            dx_parallel = 2*bpm.lam/max_NA

        print("dslm prop with parallelize beams of distance dx = %s mu"%dx_parallel)

        # the beamlet centers in the simul_xy coordinates
        ind_step = int(np.ceil(1.*dx_parallel/bpm.dx))
        #make sure its divisible by the grid size dimension
        ind_step = [i for i in range(ind_step, bpm.simul_xy[0]+1) if bpm.simul_xy[0]%i==0][0]

        if ind_step<=0:
            raise ValueError("dx resolution too coarse to propagate in parallel pleae increase dx_parallel")

        inds = np.arange(0,bpm.simul_xy[0],ind_step)

        print(inds,ind_step, bpm.simul_xy)

        u0 = np.sum([np.roll(u0_base,i,axis=1) for i in inds],axis=0)

        u = None
        # now scan
        for i in range(ind_step):
            print("propagating beamlets %s/%s"%(i+1,ind_step))
            u_part = bpm.propagate(u0,**bpm_kwargs)
            u = u_part if u is None else u+u_part

            u0 = np.roll(u0,1,axis=1)

        return self._trans_illum(u, inv = True)

if __name__ == '__main__':

    dn = np.zeros((256,512,256))

    signal = np.zeros_like(dn)

    #some point sources
    np.random.seed(0)
    for _ in range(4000):
        k,j,i = np.random.randint(dn.shape[0]),np.random.randint(dn.shape[1]),np.random.randint(dn.shape[2])
        signal[k,j,i] = 1.

    if not "m" in locals():
        m = SimLSM_DSLM(dn = dn,
                        signal = signal,
                        NA_illum= .4,
                        NA_detect=.7,
                        units = (.4,)*3,
                        #simul_xy_detect=(512,512),
                        #simul_xy_illum=(512,1024),
                         )

    u1 = m.propagate_illum(cz = -10)
    u2 = m.propagate_illum(cz = 10)

    h = m.psf((0.,0,0))

    #
    # im = m.simulate_image_z(cz=0, zslice=16,
    #                         psf_grid_dim=(16,16),
    #                         conv_sub_blocks=(8,8),
    #                         conv_pad_factor=3,
    #                         )
    #
    #

max_stars_repo_name: imagejan/biobeam

max_issues_repo_name: imagejan/biobeam

max_forks_repo_name: imagejan/biobeam


text:
from __future__ import unicode_literals

import numpy as np
import pytest

import properties
import properties.extras
from lfview.resources import files

def test_empty_init():
    arr = files.Array()
    with pytest.raises(properties.ValidationError):
        arr.validate()

def test_object_without_array():
    arr = files.Array()
    arr.content_length = 64
    arr.dtype = 'Int8Array'
    arr.shape = [64]
    assert arr.validate()

@pytest.mark.parametrize(
    'input_arr', [
        [1, 2, 3],
        (1, 2, 3),
        np.array([1, 2, 3]),
    ]
)
def test_array_init(input_arr):
    arr = files.Array(input_arr)
    assert arr.validate()
    assert isinstance(arr.array, np.ndarray)
    assert np.array_equal(arr.array, np.array([1, 2, 3]))
    assert arr.dtype.startswith('Int')
    assert arr.shape == [3]
    assert arr.content_type == 'application/octet-stream'
    assert arr.content_length == np.dtype(arr.array.dtype).itemsize * 3
    assert arr.is_1d()

@pytest.mark.parametrize(
    'input_arr', [
        ['a', 'b', 'c'],
        ('a', 'b', 'c'),
        np.array(['a', 'b', 'c']),
    ]
)
def test_array_assignment(input_arr):
    arr = files.Array()
    with pytest.raises(ValueError):
        arr.array = input_arr

@pytest.mark.parametrize(
    'input_arr', [
        [[1., 2], [3, 4]],
        ((1., 2), (3, 4)),
        np.array([[1., 2], [3, 4]]),
    ]
)
def test_hasarray_assignment(input_arr):
    class HasArray(properties.extras.HasUID):
        arr = properties.extras.Pointer('array', files.Array)

    has_array = HasArray()
    has_array.arr = input_arr
    assert has_array.validate()
    assert isinstance(has_array.arr.array, np.ndarray)
    assert np.allclose(has_array.arr.array, np.array([[1., 2], [3, 4]]))
    assert has_array.arr.dtype.startswith('Float')
    assert has_array.arr.shape == [2, 2]
    assert has_array.arr.content_type == 'application/octet-stream'
    dtype = np.dtype(has_array.arr.array.dtype)
    assert has_array.arr.content_length == dtype.itemsize * 4
    assert not has_array.arr.is_1d()

def test_array_serialize():
    arr = files.Array([[1]])
    arr.validate()
    serial_arr = arr.serialize()
    assert serial_arr == {
        '__class__': 'Array',
        'content_length': np.dtype(arr.array.dtype).itemsize,
        'content_type': 'application/octet-stream',
        'dtype': 'Int32Array',
        'shape': [1, 1],
    }

def test_unsupported_array():
    with pytest.raises(ValueError):
        files.Array('bad array')
    int64array = np.array([100000000000])
    with pytest.raises(ValueError):
        files.Array(int64array)

def test_bad_input():
    with pytest.raises(properties.ValidationError):
        files.Array(dtype='bad')
    arr = files.Array(dtype='Int32Array', shape=[1])
    with pytest.raises(properties.ValidationError):
        arr.content_length = 100

@pytest.mark.parametrize(
    ('input', 'is_1d'),
    [(None, False), ([1, 2, 3], True), ([[1, 2], [3, 4]], False)]
)
def test_is_1d(input, is_1d):
    assert files.Array(input).is_1d() == is_1d

def test_types():
    assert files.Array.BASE_TYPE == 'files'
    assert files.Array.SUB_TYPE == 'array'

if __name__ == '__main__':
    pytest.main()

max_stars_repo_name: seequent/lfview-resources-files

max_issues_repo_name: seequent/lfview-resources-files

max_forks_repo_name: seequent/lfview-resources-files


text:
# -*- coding: mbcs -*-
"""
Functions for read Abaqus odb files.

Generate plain data about:

* Parts
* Instances
* Materials
* Steps
* Nodes coordinates
* Deformed shapes
* Max. Von Mises Stresses
* Max. Nominal Strains
* Max. Plastic EQV. Strains
* Max. Reactions forces
"""

import_error_str = """
Please try run in Abaqus/CAE script mode, as follows:

>> abaqus cae noGUI="test_script.py"

or verify that you have Abaqus installed.
"""

import numpy as np
import Tkinter
import tkMessageBox
try:
    from odbAccess import *
    from abaqusConstants import *
except ImportError:
    print import_error_str
    pass

def get_max_values(dbpath,var,fname):
    """
    Get Max value from field variable (all steps - all nodes)

    Parameters
    ----------
    dbpath  : str
        Abaqus ODB file
    fname   : str
        Filename for output data
    """
    odb = openOdb(path=dbpath)
    if not odb.steps[find_last_step(dbpath)].frames[-1].fieldOutputs.has_key(var):
        root = Tkinter.Tk().withdraw()
        tkMessageBox.showerror("PyQus 0.1.0", "Undefined field variable '%s'"%(var))
        return False
    f = open(fname,"w")
    k = 0
    for step in odb.steps.keys():
        tbystep = odb.steps[step].frames[-1].frameValue
        for frame in odb.steps[step].frames:
            temp_list = []
            for td in frame.fieldOutputs[var].values:
                data = {"VECTOR_2D":td.magnitude,
                        "VECTOR_3D":td.magnitude,
                        "TENSOR_2D":td.maxPrincipal,
                        "TENSOR_3D":td.maxPrincipal,
                        "DEFAULT": 0}
                if td.type == VECTOR and len(td.data)==3:
                    tp = "VECTOR_3D"
                elif td.type == VECTOR:
                    tp = "VECTOR_2D"
                elif td.type == TENSOR_2D_SURFACE:
                    tp = "TENSOR_2D"
                elif td.type == TENSOR_3D_FULL:
                    tp = "TENSOR_3D"
                else:
                    tp = "DEFAULT"
                temp_list.append(data[tp])
            maxv = max(temp_list) # Max value
            f.write("%g,%g\n"%(frame.frameValue + k*tbystep, maxv))
        k += 1
    f.close()
    odb.close()

def get_field_var(dbpath,var,inst,stepname,nframe=-1,fname="field_out.txt"):
    """
    Get field output

    Parameters
    ----------
    dbpath   : str
        Abaqus ODB file
    inst     : str
        Instance name
    stepname : str
        Step name
    nframe   : int
        Number of frame
    fname    : str
        Filename for output data
    """

    odb = openOdb(path=dbpath)
    _inst = odb.rootAssembly.instances[inst]
    if not odb.steps[stepname].frames[nframe].fieldOutputs.has_key(var):
        root = Tkinter.Tk().withdraw()
        tkMessageBox.showerror("PyQus 0.1.0", "Undefined field variable '%s'"%(var))
        return False
    f = open(fname,"w")
    valofinst = odb.steps[stepname].frames[nframe].fieldOutputs[var].getSubset(region=_inst).values
    for k,jj in enumerate(valofinst):
        data = {"VECTOR_2D":["%g,%g,%g,%g\n",(jj.nodeLabel,jj.data[0],jj.data[1],jj.magnitude)],
                "VECTOR_3D":["%g,%g,%g,%g,%g\n",(jj.nodeLabel,jj.data[0],jj.data[1],jj.data[2],jj.magnitude)],
                "TENSOR_2D":["%g,%g,%g\n",(jj.elementLabel,jj.maxPrincipal,jj.mises)],
                "TENSOR_3D":["%g,%g,%g\n",(jj.elementLabel,jj.maxPrincipal,jj.mises)],
                "DEFAULT":["%s\n",(jj)]}

        if jj.type == VECTOR and len(jj.data)==2:
            tp = "VECTOR_2D"
        elif jj.type == VECTOR and len(jj.data)==3:
            tp = "VECTOR_3D"
        elif jj.type == TENSOR_2D_SURFACE:
            tp = "TENSOR_2D"
        elif jj.type == TENSOR_3D_FULL:
            tp = "TENSOR_3D"
        else:
            tp = "DEFAULT"
        f.write(data[tp][0]%data[tp][1])

    f.close()
    odb.close()

def get_nodes_coordinates(dbpath,inst,stepname,nframe=-1,fname="nodes.txt"):
    """
    Get nodes coordinates from meshed instance // last frame by default

    Parameters
    ----------
    dbpath   : str
        Abaqus ODB file
    inst     : str
        Instance name
    stepname : str
        Step name
    nframe   : int
        Number of frame
    fname    : str
        Filename for output data
    """
    odb = openOdb(path=dbpath)
    f = open(fname,"w")
    _inst = odb.rootAssembly.instances[inst]
    ic = odb.rootAssembly.instances[inst].nodes # Initial coordinates
    valofinst = odb.steps[stepname].frames[nframe].fieldOutputs['U'].getSubset(region=_inst).values
    for k,jj in enumerate(valofinst):
        if _inst.embeddedSpace == THREE_D: # 3D case
            #f.write("%g,"%(ic[k].label)) # Write node label ?
            f.write("%0.4f,%0.4f,%0.4f\n"%(ic[k].coordinates[0]+jj.data[0],
                                           ic[k].coordinates[1]+jj.data[1],
                                           ic[k].coordinates[2]+jj.data[2]))
        else: # 2D case
            f.write("%0.4f,%0.4f\n"%(ic[k].coordinates[0]+jj.data[0],
                                     ic[k].coordinates[1]+jj.data[1]))
    f.close()

def get_nodes_distance(dbpath,node1,node2,inst,stepname,nframe=-1):
    """
    Calculate distance between two nodes

    Parameters
    ----------
    dbpath   : str
        Abaqus ODB file
    node1    : int
        Number of node 1
    node2    : int
        Number of node 2
    inst     : str
        Instance name
    stepname : str
        Step name
    nframe   : int
        Number of frame, -1 for last frame (See Python List Indexing)
    """
    odb = openOdb(path=dbpath)
    _inst = odb.rootAssembly.instances[inst]
    ic = odb.rootAssembly.instances[inst].nodes
    us = odb.steps[stepname].frames[nframe].fieldOutputs['U'].getSubset(region=_inst).values
    xx1 = ic[node1-1].coordinates[0]+us[node1-1].data[0]
    yy1 = ic[node1-1].coordinates[1]+us[node1-1].data[1]
    xx2 = ic[node2-1].coordinates[0]+us[node2-1].data[0]
    yy2 = ic[node2-1].coordinates[1]+us[node2-1].data[1]
    if _inst.embeddedSpace == THREE_D:
        zz1 = ic[node1-1].coordinates[2]+us[node1-1].data[2]
        zz2 = ic[node2-1].coordinates[2]+us[node2-1].data[2]
        d = np.sqrt((xx2-xx1)**2 + (yy2-yy1)**2 + (zz2-zz1)**2)
    else:
        d = np.sqrt((xx2-xx1)**2+(yy2-yy1)**2)
    return d

def get_distances(dbpath,M1,M2,inst,stepname,nframe=-1,filename="distances.txt"):
    """
    Get distances between two arrays of nodes. For example, if M1 and M2
    are two arrays of nodes defined by:

    ::

        M1 = (1,2,3)
        M2 = (4,5,6)

    it returns the distances between 1-4, 2-5 and 3-6 nodes.

    Parameters
    ----------
    dbpath  : str
        Abaqus ODB file
    M1      : tuple, list
        Arrays of nodes, i.e.  M1 = (3,9,31)
    M2      : tuple, list
        Arrays of nodes, i.e.  M2 = (10,18,38)
    inst    : str
        Instance name
    stepname: str
        Step name
    nframe  : int
        Number of frame, -1 for last frame (See Python List Indexing)
    filename: str
        Name of output file
    """
    f = open(filename,"w")
    for ii,jj in enumerate(M1):
        f.write("%4.6f\n"%(get_nodes_distance(dbpath,jj,M2[ii],inst,stepname,nframe),))
    f.close()

def get_elements(dbpath,inst,fname="elements.txt"):
    """
    Get connectivity elements.

    Parameters
    ----------
    dbpath : str
        Abaqus ODB file
    inst   : str
        Instance name
    fname  : str
        Name of output file
    """
    odb = openOdb(path=dbpath)
    f = open(fname,"w")
    _inst = odb.rootAssembly.instances[inst]
    for el in _inst.elements:
        econn = el.connectivity
        nel = len(econn)
        str_flag = (nel*"%g,")[:-1]+"\n"
        elmts = tuple([econn[k] for k in range(nel)])
        f.write(str_flag%elmts)
    f.close()

def find_deformable_body(dbpath):
    """
    Return name of first deformable body in the analysis.

    Parameters
    ----------
    dbpath : str
        Abaqus ODB file
    """
    odb = openOdb(path=dbpath)
    for name,_inst in odb.rootAssembly.instances.items():
        if _inst.type == DEFORMABLE_BODY and name!="ASSEMBLY":
            return name
    return None

def find_last_step(dbpath):
    """
    Return last step name from odb file

    Parameters
    ----------
    dbpath : str
        Abaqus ODB file
    """
    #~ if not isinstance(dbpath,OdbType):
    odb = openOdb(path=dbpath)
    _steps = odb.steps.keys()
    last_step = _steps[-1]
    # odb.close() <- This failed
    return last_step

def get_steps(dbpath):
    """
    Return all steps names from odb file, also include
    additional information as time duration, number of frames
    per step and procedure type.

    Parameters
    ----------
    dbpath : str
        Abaqus ODB file
    """
    odb = openOdb(path=dbpath)
    _steps = []
    for _name,_stp in odb.steps.items():
        _time = _stp.timePeriod
        _procedure = _stp.procedure
        _nframes = len(_stp.frames)
        _steps.append((_name,_time,_nframes,_procedure))
    odb.close()
    return _steps

def get_materials(dbpath):
    """
    Return all materials names from odb file

    Parameters
    ----------
    dbpath : str
        Abaqus ODB file
    """
    odb = openOdb(path=dbpath)
    _materials = []
    for _name in odb.materials.items():
        _materials.append(_name)
    odb.close()
    return _materials

def get_materials_properties(dbpath): #<un-named>nook
    """
    Get elastic/plastic properties for all materials in the odb file.
    Return a list of tuples, one tuple for every material.

    Parameters
    ----------
    dbpath : str
        Abaqus ODB file
    """
    odb = openOdb(path=dbpath)
    data = []
    for _name,_mat in odb.materials.items():
        _elastic_mod = _mat.elastic.table[0][0]
        _poisson = _mat.elastic.table[0][1]
        if hasattr(_mat,"plastic"):
            _plastic = _mat.plastic.table
        else:
            _plastic = []
        data.append((_name,_elastic_mod,_poisson,_plastic))
    odb.close()
    return data

def get_parts(dbpath):
    """
    Get parts names and types. Return a list of tuples.

    Parameters
    ----------
    dbpath : str
        Abaqus ODB file
    """
    odb = openOdb(path=dbpath)
    _parts = []
    for _name,_prt in odb.parts.items():
        _type = _prt.type
        _parts.append((_name,_type))
    return _parts

def get_instances(dbpath):
    """
    Get instances names. Return a list of tuples.

    Parameters
    ----------
    dbpath : str
        Abaqus ODB file
    """
    odb = openOdb(path=dbpath)
    _instances = []
    for _name,_inst in odb.rootAssembly.instances.items():
        _nodes = len(_inst.nodes)
        _elements = len(_inst.elements)
        _instances.append((_name,_nodes,_elements))
    return _instances

def get_job(dbpath):
    """
    Get job info:

    - Analysis code (type of analysis i.e. Dynamic/Explicit)
    - Creation time (date + time string)
    - Precision (double/single)

    Parameters
    ----------
    dbpath : str
        Abaqus ODB file
    """
    odb = openOdb(path=dbpath)
    analysis_code = odb.jobData.analysisCode
    creation_time = odb.jobData.creationTime
    precision = odb.jobData.precision
    return [(analysis_code,creation_time,precision)]

if __name__ == '__main__':
    pass

max_stars_repo_name: JorgeDeLosSantos/pyqus

max_issues_repo_name: JorgeDeLosSantos/pyqus

max_forks_repo_name: JorgeDeLosSantos/pyqus


text:
# -*- coding: utf-8 -*-
from six.moves import range
import numpy as np

from ..general.neighbours import neighbours
from ..general.coord2K import coord2K
from ..stest.designmatrix import designmatrix

def kriging(ck,ch,zh,model,param,nhmax,dmax,order=None,options=None):
  '''
% kriging                   - prediction using kriging methods
%
% Standard linear kriging algorithm for processing hard data
% that can reasonably be assumed as Gaussian distributed. The
% function is intendend to be as general as possible, covering
% various situations, like non-stationarity of the mean,
% multivariate cases, nested models, space-time estimations,
% etc. Depending on the case, specific format are needed for
% the input variables. This function is a special case of the more
% general BMEprobaMoments.mfunction, which processes both hard and
% soft data.
%
% SYNTAX :
%
% [zk,vk]=kriging(ck,ch,zh,model,param,nhmax,dmax,order,options);
%
% INPUT :
%
% ck        nk by d   matrix of coordinates for the estimation locations.
%                     A line corresponds to the vector of coordinates at
%                     an estimation location, so the number of columns
%                     corresponds to the dimension of the space. There is
%                     no restriction on the dimension of the space.
% ch        nh by d   matrix of coordinates for the hard data locations,
%                     with the same convention as for ck.
% zh        nh by 1   vector of values for the hard data at the coordinates
%                     specified in ch.
  covmodel  list      list of m nested covariance models in each of
                      which the spatial and temporal components are put
                      in a list as [covmodelS,covmodelT]
  covparam0 list      list of intial values covariance parameters for
                      m covmodels in each component the parameters are
                      listed as [sill,[covparamS1,covparamS2,..],
                      [covparamT1,covparamT2,..]]
% nhmax     scalar    maximum number of hard data values that are considered
%                     for the estimation at the locations specified in ck.
% dmax      scalar    maximum distance between an estimation location and
%                     existing hard data locations. All hard data locations
%                     separated by a distance smaller than dmax from an
%                     estimation location will be included in the estimation
%                     process for that location, whereas other hard data
%                     locations are neglected.
% order     scalar    order of the polynomial mean along the spatial axes at
%                     the estimation locations. For the zero-mean case, NaN
%                     (Not-a-Number) is used. Note that order=NaN can only be
%                     used with covariance models and not with variogram models.
% options   scalar    optional parameter that can be used if the default value
%                     is not satisfactory (otherwise it can simply be omitted
%                     from the input list of variables). options(1) is taking
%                     the value 1 or 0 if the user wants or does not want to
%                     display the order number of the location which is
%                     currently processed, respectively.
%
% OUTPUT :
%
% zk        nk by 1   vector of estimated values at the estimation locations. A
%                     value coded as NaN means that no estimation has been performed
%                     at that location due to the lack of available data.
% vk        nk by 1   vector of estimation (kriging) variances at the estimation
%                     locations. As for zk, a value coded as NaN means that no
%                     estimation has been performed at the corresponding location.
%

    input
    ck: nk by d float
    ch: nh by d float
    zh: nh by 1 float
    model: mn by 1 string
    param: mn by 3 float
    nhmax: int
    dmax: 1 by 3 float
    order: None or 0
    options: not no use here

    return
    zk: nk by 1 float
    vk: nk by 1 float
  '''

    #initial return variable
  nk = ck.shape[0]
  zk = np.zeros( ( nk, 1 ) )
  vk = np.zeros( ( nk, 1 ) )
  zk[:] = np.nan
  vk[:] = np.nan

  if order is None:
    order = np.nan

  for i in range( ck.shape[0] ):
    ci = ck[i:i+1]
    ci_nebr,zi_nebr,di_nebr,ni_nebr,idxi_nebr=neighbours(ci,ch,zh,nhmax,dmax)
    if ni_nebr > 0:
      K, dummyKK = coord2K( ci_nebr, ci_nebr, model, param )
      k, dummykk = coord2K( ci_nebr, ci, model, param ) #kk is a list contain ki
      k0, dummykk0 = coord2K(ci, ci, model, param)
      X=designmatrix(ci_nebr,order)[0]
      x=designmatrix(ci,order)[0]
#            unit = np.ones(k.shape) # n by 1, X in matlab
#            unit_t_add = np.append( unit.T, [[0]], axis = 1 ) # 1 by n+1, x in matlab

            #change shape for kriging
#            Kadd = np.hstack( ( K, unit ) )
#            Kadd = np.vstack( ( Kadd, unit_t_add ) )
      Kadd = K
      Xm,Xn = X.shape
      Kadd = np.vstack([np.hstack([K,X]),
                           np.hstack([X.T,np.zeros((Xn,Xn))])])
#            kadd = np.append( k, [[0]], axis = 0 )
      kadd = np.vstack([k,x.T])
      weight = np.linalg.solve(Kadd,kadd)[0:Xm]
      weight_t = weight.T
            #compute zk, vk
      zk[i] = weight_t.dot( zi_nebr )
      vk[i] = (k0-2*weight_t.dot(k)+weight_t.dot(K).dot(weight))[0]

    else:
      pass #already give NaN

  return zk, vk

def krigconstr(c0,c,order=None):
  '''
  order is None or 0 for now
  '''

  return X,x

def krigingFac( ck, ch, zh, model, param, nhmax, dmax, order = 0, options = None):

    '''
    input
    ck: nk by d float
    ch: nh by d float
    zh: nh by 1 float
    model: mn by 1 string ( ex. "exponential/gaussian" )
    param: mn by 3 float (ex. [ c, s, t ] )
    nhmax: int
    dmax: 1 by 3 float
    order: always 0, equals NaN in matlab
    options: not no use here, 0 or 1 in matlab for display echo

    return
    zk: mn+1 by 1 float
    vk: mn+1 by 1 float
    '''

    #initial return variable
    nk = ck.shape[0]
    mn = len(model) #model.shape[0]
    zk = np.zeros( ( nk, mn+1 ) )
    vk = np.zeros( ( nk, mn+1 ) )
    zk[:] = np.nan
    vk[:] = np.nan

    for i in range( ck.shape[0] ):
        ci = ck[i:i+1]
        ci_nebr, zi_nebr, di_nebr, ni_nebr, idxi_nebr = neighbours(ci, ch, zh, nhmax, dmax )
        if ni_nebr > 0:
            K, dummyKK = coord2K( ci_nebr, ci_nebr, model, param )
            dummyk, kk = coord2K( ci_nebr, ci, model, param ) #kk is a list contain ki
            dummyk0, kk0 = coord2K(ci, ci, model, param)
            unit = np.ones(kk[0].shape) # n by 1
            unit_t_add = np.append( unit.T, [[0]], axis = 1 ) # 1 by n+1
            for idx_k, (ki,k0i) in enumerate( zip( kk, kk0 ) ):
                #change shape for kriging
                Kadd = np.hstack( ( K, unit ) )
                Kadd = np.vstack( ( Kadd, unit_t_add ) )
                kkadd = np.append( ki, [[0]], axis = 0 )
                weight = np.dot( np.linalg.inv( Kadd ), kkadd )[:-1,:]
                weight_t = weight.T
                #compute zk, vk
                zk[i:i+1,idx_k:idx_k+1] = weight_t.dot( zi_nebr )
                vk[i:i+1,idx_k:idx_k+1] = ( k0i - 2 * weight_t.dot( ki ) + weight_t.dot( K ).dot( weight ) )[0]
            #compute local mean trend
            kkadd = np.append( np.zeros( kk[0].shape ), [[1]], axis = 0 )
            weight = np.dot( np.linalg.inv( Kadd ), kkadd )[:-1,:]
            weight_t = weight.T
            zk[i:i+1,-1:] = weight_t.dot( zi_nebr )
            vk[i:i+1,-1:] = weight_t.dot( K ).dot( weight )[0]
        else:
            pass #already give NaN
    return zk, vk