text:
import os

import numpy as np
import pandas as pd

class ExplainerMixin(object):

    def __init__(
            self,
            path,
            data,
            features
    ):
        if not os.path.exists(path):
            os.makedirs(path)

        self.path = path
        self.data = data
        self.features = features

    @property
    def data_is_2d(self):
        if isinstance(self.data, np.ndarray) and self.data.ndim == 2:
            return True
        elif isinstance(self.data, pd.DataFrame):
            return True
        else:
            return False

    @property
    def data_is_3d(self):
        if isinstance(self.data, np.ndarray) and self.data.ndim == 3:
            return True
        return False

    @property
    def single_source(self):
        if isinstance(self.data, list) and len(self.data) > 1:
            return False
        else:
            return True

    @property
    def features(self):
        return self._features

    @features.setter
    def features(self, features):
        if self.data_is_2d:
            if type(self.data) == pd.DataFrame:

                features = self.data.columns.to_list()
            elif features is None:
                features = [f"Feature {i}" for i in range(self.data.shape[-1])]
            else:
                assert isinstance(features, list) and len(features) == self.data.shape[-1], f"""
                    features must be given as list of length {self.data.shape[-1]}
                    but are of len {len(features)}
                    """

                features = features
        elif not self.single_source and features is None:
            features = []
            for data in self.data:
                if isinstance(data, pd.DataFrame):
                    _features = data.columns.to_list()
                else:
                    _features = [f"Feature {i}" for i in range(data.shape[-1])]

                features.append(_features)

        elif self.data_is_3d and features is None:
            features = [f"Feature {i}" for i in range(self.data.shape[-1])]

        self._features = features

    @property
    def unrolled_features(self):
        # returns the possible names of features if unrolled over time
        if not self.data_is_2d and self.single_source:
            features = self.features
            if features is None:
                features = [f"Feature_{i}" for i in range(self.data.shape[-1])]

            lookback = self.data.shape[1]
            features = [[f"{f}_{i}" for f in features] for i in range(lookback)]
            features = [item for sublist in features for item in sublist]
        else:
            features = None
        return features


text:
from pathlib import Path

import numpy as np

from PIL import Image, ImageFile
import cv2

from torchvision import transforms
from torch.utils.data import Dataset

from util.kaggle_util import is_kagglekernel
ImageFile.LOAD_TRUNCATED_IMAGES = True

class BencolorTrainDataset(Dataset):
    def __init__(self, img_df, transform=transforms.ToTensor(), binary=False):
        self.img_df = img_df
        self.transform = transform
        self.binary = binary

        if is_kagglekernel():
            self.data_path = Path(__file__).parents[4] / "aptos2019-blindness-detection" / "train_images"
            self.cache_path = Path(__file__).parents[5] / "working" / "data" / "bencolor"
        else:
            self.data_path = Path(__file__).parents[2] / "input" / "train_images"
            self.cache_path = Path(__file__).parents[2] / "data" / "bencolor"

        if not self.cache_path.exists():
            self.cache_path.mkdir(parents=True)

    def __len__(self):
        return len(self.img_df)

    def __getitem__(self, idx):
        cache_filepath = self.cache_path / "{}.png".format(self.img_df.loc[idx, "id_code"])

        if not cache_filepath.exists():
            self.__save_cache(idx)

        img = Image.open(str(cache_filepath))
        img = self.transform(img)

        label = self.img_df.loc[idx, "diagnosis"]
        if self.binary:
            label = 1 if label > 0 else 0

        return {"image": img, "label": label}

    def __save_cache(self, idx):
        img_filepath = self.data_path / "{}.png".format(self.img_df.loc[idx, "id_code"])
        img = cv2.imread(str(img_filepath.resolve()))

        cropper = BenColorCrop()
        img = cropper(img)

        cache_filepath = self.cache_path / "{}.png".format(self.img_df.loc[idx, "id_code"])
        cv2.imwrite(str(cache_filepath), img)

class BencolorTestDataset(Dataset):
    def __init__(self, img_df, transform=transforms.ToTensor()):
        self.img_df = img_df
        self.transform = transform

        if is_kagglekernel():
            self.data_path = Path(__file__).parents[4] / "aptos2019-blindness-detection" / "test_images"
            self.cache_path = Path(__file__).parents[5] / "working" / "data" / "bencolor"
        else:
            self.data_path = Path(__file__).parents[2] / "input" / "test_images"
            self.cache_path = Path(__file__).parents[2] / "data" / "bencolor"

        if not self.cache_path.exists():
            self.cache_path.mkdir(parents=True)

    def __len__(self):
        return len(self.img_df)

    def __getitem__(self, idx):
        cache_filepath = self.cache_path / "{}.png".format(self.img_df.loc[idx, "id_code"])

        if not cache_filepath.exists():
            self.__save_cache(idx)

        img = Image.open(str(cache_filepath.resolve()))
        img = self.transform(img)

        return {"image": img}

    def __save_cache(self, idx):
        img_filepath = self.data_path / "{}.png".format(self.img_df.loc[idx, "id_code"])
        img = cv2.imread(str(img_filepath.resolve()))

        cropper = BenColorCrop()
        img = cropper(img)

        cache_filepath = self.cache_path / "{}.png".format(self.img_df.loc[idx, "id_code"])
        cv2.imwrite(str(cache_filepath), img)

class BenColorCrop(object):
    def __init__(self, img_size=512, sigmaX=10, tol=7):
        self.img_size = img_size
        self.sigmaX = sigmaX
        self.tol = tol

    def __call__(self, img):
        img = self._load_ben_color(img)
        return img

    def _load_ben_color(self, img):
        img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
        img = self._crop_image_from_gray(img)
        img = cv2.resize(img, (self.img_size, self.img_size))
        img = cv2.addWeighted(img, 4, cv2.GaussianBlur(img, (0, 0), self.sigmaX), -4, 128)
        img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
        return img

    def _crop_image_from_gray(self, img):
        if img.ndim == 2:
            mask = img > self.tol
            return img[np.ix_(mask.any(1), mask.any(0))]
        elif img.ndim == 3:
            gray_img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY)
            mask = gray_img > self.tol

            check_shape = img[:, :, 0][np.ix_(mask.any(1), mask.any(0))].shape[0]
            if check_shape == 0:
                return img
            else:
                img1 = img[:, :, 0][np.ix_(mask.any(1), mask.any(0))]
                img2 = img[:, :, 1][np.ix_(mask.any(1), mask.any(0))]
                img3 = img[:, :, 2][np.ix_(mask.any(1), mask.any(0))]
                img = np.stack([img1, img2, img3], axis=-1)
            return img

ext: py

reason: import numpy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Tools for working with the FSPS filter set.

This module uses filter information shipped with FSPS itself in
``$SPS_HOME/data``.
"""

from __future__ import (division, print_function, absolute_import,
                        unicode_literals)

__all__ = ["find_filter", "FILTERS", "get_filter", "list_filters"]

import os
import numpy as np
from pkg_resources import resource_stream, resource_exists

# Cache for $SPS_HOME/data/magsun.dat parsed by numpy
MSUN_TABLE = None

# Cache for $SPS_HOME/data/filter_lambda_eff.dat parsed by numpy
LAMBDA_EFF_TABLE = None

# Cache for bandpass transmission listings: a dictionary keyed by bandpass
# name with values of wavelength, transmission tuples.
TRANS_CACHE = None

class Filter(object):

    def __init__(self, index, name, fullname):
        self.index = index - 1
        self.name = name.lower()
        self.fullname = fullname

    def __str__(self):
        return "<Filter({0})>".format(self.name)

    def __repr__(self):
        return "<Filter({0})>".format(self.name)

    @property
    def msun_ab(self):
        """Solar absolute magnitude in Filter, AB zeropoint."""
        # if self._msun_ab is None:
        if MSUN_TABLE is None:
            self._load_msun_table()
        if self.index < MSUN_TABLE.shape[0]:
            assert MSUN_TABLE[self.index, 0] == self.index + 1
            return float(MSUN_TABLE[self.index, 1])
        else:
            return np.nan

    @property
    def msun_vega(self):
        """Solar absolute magnitude in Filter, VEGAMAG zeropoint."""
        if MSUN_TABLE is None:
            self._load_msun_table()
        if self.index < MSUN_TABLE.shape[0]:
            assert MSUN_TABLE[self.index, 0] == self.index + 1
            return float(MSUN_TABLE[self.index, 2])
        else:
            return np.nan

    @property
    def lambda_eff(self):
        """Effective wavelength of Filter, in Angstroms."""
        if LAMBDA_EFF_TABLE is None:
            self._load_lambda_eff_table()
        if self.index < LAMBDA_EFF_TABLE.shape[0]:
            assert LAMBDA_EFF_TABLE[self.index, 0] == self.index + 1
            return float(LAMBDA_EFF_TABLE[self.index, 1])
        else:
            return np.nan

    @property
    def transmission(self):
        """Returns filter transmission: a tuple of wavelength (Angstroms) and
        an un-normalized transmission arrays.
        """
        if TRANS_CACHE is None:
            # Load the cache for all filters.
            self._load_transmission_cache()
        try:
            return TRANS_CACHE[self.name]
        except KeyError as e:
            e.args += ("Could not find transmission data "
                       "for {0}".format(self.name))
            raise

    def _load_msun_table(self):
        global MSUN_TABLE
        MSUN_TABLE = np.loadtxt(
            os.path.expandvars("$SPS_HOME/data/magsun.dat"))

    def _load_lambda_eff_table(self):
        global LAMBDA_EFF_TABLE
        LAMBDA_EFF_TABLE = np.loadtxt(
            os.path.expandvars("$SPS_HOME/data/filter_lambda_eff.dat"))

    def _load_transmission_cache(self):
        """Parse the allfilters.dat file into the TRANS_CACHE."""
        global TRANS_CACHE
        path = os.path.expandvars("$SPS_HOME/data/allfilters.dat")
        names = list_filters()
        TRANS_CACHE = {}
        filter_index = -1
        lambdas, trans = [], []
        with open(path) as f:
            for line in f:
                line.strip()
                if line[0].startswith("#"):
                    # Close out filter
                    if filter_index > -1:
                        TRANS_CACHE[names[filter_index]] = (
                            np.array(lambdas), np.array(trans))
                    # Start new filter
                    filter_index += 1
                    lambdas, trans = [], []
                else:
                    try:
                        l, t = line.split()
                        lambdas.append(float(l))
                        trans.append(float(t))
                    except(ValueError):
                        pass

def _load_filter_dict():
    """
    Load the filter list, creating a dictionary of :class:`Filter` instances.
    """
    # Load filter table from FSPS
    filter_list_path = os.path.expandvars(
        os.path.join("$SPS_HOME", "data", "FILTER_LIST"))
    filters = {}
    with open(filter_list_path) as f:
        for line in f:
            columns = line.strip().split()
            if len(columns) > 1:
                fsps_id, key = columns[:2]
            else:
                continue
            comment = ' '.join(columns[2:])
            filters[key.lower()] = Filter(int(fsps_id), key, comment)

    return filters

FILTERS = _load_filter_dict()

def find_filter(band):
    """
    Find the FSPS name for a filter.

    Usage:

    ::

        >>> import fsps
        >>> fsps.find_filter("F555W")
        ['wfpc2_f555w', 'wfc_acs_f555w']

    :param band:
        Something like the name of the band.

    """
    b = band.lower()
    possible = []
    for k in FILTERS.keys():
        if b in k:
            possible.append(k)
    return possible

def get_filter(name):
    """Returns the :class:`fsps.filters.Filter` instance associated with the
    filter name.

    :param name:
        Name of the filter, as found with :func:`find_filter`.
    """
    try:
        return FILTERS[name.lower()]
    except KeyError as e:
        e.args += ("Filter {0} does not exist. "
                   "Try using fsps.find_filter('{0}').".format(name),)
        raise

def list_filters():
    """Returns a list of all FSPS filter names.

    Filters are sorted by their FSPS index.
    """
    lst = [(name, f.index) for name, f in FILTERS.items()]
    lst.sort(key=lambda x: x[1])
    return [l[0] for l in lst]

max_stars_repo_name: weiguangcui/python-fsps

max_issues_repo_name: weiguangcui/python-fsps

max_forks_repo_name: weiguangcui/python-fsps


text:
from time import sleep, time
import cv2
from scipy.ndimage.filters import convolve as filter2, gaussian_filter
import numpy as np

#
windowAvg =np.array([[1 / 12, 1 / 6, 1 / 12],
                     [1/6,    0, 1/6],
                     [1/12, 1/6, 1/12]], float)

windowX = np.array([[-1, 1],
                    [-1, 1]]) * .25 #kernel for computing d/dx

# windowX = np.array([[-1, 0,1],
#                     [-1, 0,1],
#                     [-1,0,1]]) * .25 #kernel for computing d/dx
#

windowY = np.array([[-1, -1],
                    [ 1, 1]]) * .25 #kernel for computing d/dy
#
# windowY = np.array([[-1, -1,-1],
#                     [0, 0,0],
#                     [1,1,1]]) * .25 #kernel for computing d/dx

windowT = np.ones((2, 2)) * .25

class HornSchunck:

    def process(self,frame1, frame2, alpha=0.001, NumOfIter=8):
        """
        frame1: frame at t=0
        frame2: frame at t=1
        alpha: regularization constant
        NumOfIter: number of iteration
        """
        #if the frame is integers then we need to convert it to floats
        frame1 = frame1.astype(np.float32)
        frame2 = frame2.astype(np.float32)

        # making the shape of horizontal and vertical change
        # Set initial value for the flow vectors
        H = np.zeros([frame1.shape[0], frame1.shape[1]])
        V = np.zeros([frame1.shape[0], frame1.shape[1]])

        # Estimate derivatives
        [fx, fy, ft] = self.derivatives(frame1, frame2)

        # Iteration to reduce error
        for i in range(NumOfIter):
            # avrageing the flow vectors
            hAvg = cv2.filter2D(H, -1, windowAvg)
            vAvg = cv2.filter2D(V, -1, windowAvg)
            # common part of update step
            top = fx*hAvg + fy*vAvg + ft
            down = alpha**2 + fx**2 + fy**2
            der = top/down

            # iterative step
            H = hAvg - fx * der
            V = vAvg - fy * der

        M = pow(pow(H, 2) + pow(V, 2), 0.5)

        return H,V, M

    def derivatives(self,frame1, frame2):
        t = time()
        fx = filter2(frame1, windowX) + filter2(frame2, windowX)
        fy = filter2(frame1, windowY) + filter2(frame2, windowY)
       # ft = im2 - im1
        ft = filter2(frame1, windowT) + filter2(frame2, -windowT)
        # print(time() - t)
        return fx,fy,ft

    def draw_vectors_hs(self,im1, im2, step = 10):
        print("drawing vectors")
        t = time()
        im1_gray = cv2.cvtColor(im1, cv2.COLOR_BGR2GRAY)
        im2_gray = cv2.cvtColor(im2, cv2.COLOR_BGR2GRAY)
        # im1_gray = cv2.GaussianBlur(im1_gray,(5,5),0)
        # im2_gray = cv2.GaussianBlur(im2_gray,(5,5),0)

        U, V, M = self.process(im1_gray, im2_gray)

        print(time() - t)
        rows, cols = im2_gray.shape
        # print(rows, cols, range(0, rows, step))
        for i in range(0, rows, step):
            for j in range(0, cols, step):
                x = int(U[i, j]*2)
                y = int(V[i, j] *2)
                cv2.arrowedLine(im2, (j, i), (j + x, i + y), (255, 0, 0))
        return im2

if __name__ == "__main__":
    cap = cv2.VideoCapture("2.mkv")
    i =0
    while(i < 200):
        ret, frame = cap.read()  # get first frame
        i+=1
    ret,old = cap.read()
    hs = HornSchunck()
    while(True):
        ret, new = cap.read()
        ret, new = cap.read()
        if not ret:
            break

        # kernel = np.ones((3, 3), np.float32) / 9
        # new  = cv2.filter2D(new, -1, kernel)
        ret = hs.draw_vectors_hs(old,new)
        old = new
        cv2.imshow("frame", ret)
        cv2.waitKey(1)

    cv2.destroyAllWindows()

def ccw(A,B,C):
	return (C[1]-A[1])*(B[0]-A[0]) > (B[1]-A[1])*(C[0]-A[0])

def intersect(A,B,C,D):
	return ccw(A,C,D) != ccw(B,C,D) and ccw(A,B,C) != ccw(A,B,D)

# buscamos si hay vectores de flujo optico que se intersectan, a fuerza bruta
def check_intersection(lines, frame):
    print("vericando intersecciones entre " + str(len(lines)))
    print(lines)
    for i in range(0, len(lines)):
        for j in range(0, len(lines)):
            #L1 = line([lines[i][0], lines[i][1]], [lines[i][2], lines[i][3]])
            #L2 = line([lines[j][0], lines[j][1]], [lines[j][2], lines[j][3]])
            #print([lines[i][0], lines[i][1]], [lines[i][2], lines[i][3]], [lines[j][0], lines[j][1]], [lines[j][2], lines[j][3]])

            a1 = (lines[i][0], lines[i][1])
            a2 = (lines[i][2], lines[i][3])
            b1 = (lines[j][0], lines[j][1])
            b2 = (lines[j][2], lines[j][3])

            if i != j:
                # R = intersection(L1, L2)
                R = intersect(a1, a2, b1, b2)
                rows, cols, ch = frame.shape

                if R:
                    print("Intersection detected", a1, a2, b1,b2)
                    overlay = frame.copy()
                    cv2.rectangle(overlay, (0, 0), (cols - 1, rows - 1), (0, 0, 255), -1)
                    opacity = 0.4
                    cv2.addWeighted(overlay, opacity, frame, 1 - opacity, 0, frame)

'''
cap = cv2.VideoCapture('normal2 .avi')
index = 0
frames = []
cv2.namedWindow('win',cv2.WINDOW_NORMAL)
while True:
    ret, frame = cap.read()

    if ret:

        frames.append(frame)

        # solo procesamos yolo una vez
        if index > 1:
            draw_vectors(frames[-2], frames[-1])

        cv2.imshow("win", frame)
        if cv2.waitKey(0) & 0xFF == ord('q'):
            break

        index += 1

    else:
        break

cv2.destroyAllWindows()
'''

''' lukas kanade
cap = cv2.VideoCapture('normal7.mp4')
cv2.namedWindow('frame',cv2.WINDOW_NORMAL)

# params for ShiTomasi corner detection
feature_params = dict( maxCorners = 100,
                       qualityLevel = 0.3,
                       minDistance = 7,
                       blockSize = 7 )

# Parameters for lucas kanade optical flow
lk_params = dict( winSize  = (15,15),
                  maxLevel = 2,
                  criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))

# Create some random colors
color = np.random.randint(0,255,(100,3))

# Take first frame and find corners in it
ret, old_frame = cap.read()
while not ret:
    print()
    ret, old_frame = cap.read()
old_gray = cv2.cvtColor(old_frame, cv2.COLOR_BGR2GRAY)
p0 = cv2.goodFeaturesToTrack(old_gray, mask = None, **feature_params)

# Create a mask image for drawing purposes
mask = np.zeros_like(old_frame)

while(1):
    ret,frame = cap.read()
    if ret:
        frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

        # calculate optical flow
        p1, st, err = cv2.calcOpticalFlowPyrLK(old_gray, frame_gray, p0, None, **lk_params)

        # Select good points
        good_new = p1[st==1]
        good_old = p0[st==1]

        lines = []
        # draw the tracks
        for i, (new, old) in enumerate(zip(good_new,good_old)):
            a, b = new.ravel()
            c, d = old.ravel()
            lines.append((a, b, c, d))
            mask = cv2.line(mask, (a, b), (c, d), color[i].tolist(), 2)
            frame = cv2.circle(frame, (a, b), 5, color[i].tolist(), -1)
            cv2.line(frame, (a, b), (c, d), (0, 0, 255), 1)
            # frame = cv2.circle(frame, (c, d), 5, (255,0 ,0), -1)

        check_intersection(lines, frame)

        img = cv2.add(frame,mask)

        cv2.imshow('frame',img)
        if cv2.waitKey() & 0xFF == ord('q'):
            break
    else:
        break

    # Now update the previous frame and previous points
    old_gray = frame_gray.copy()
    p0 = good_new.reshape(-1,1,2)

cv2.destroyAllWindows()
cap.release()
'''

reason: import numpy,from scipy


text:
## plot Zipfs Law

from nltk.corpus import gutenberg
from nltk import FreqDist
from scipy.stats import lognorm
import matplotlib
import matplotlib.pyplot as plt

all_words = gutenberg.words()

some_words = gutenberg.words('austen-persuasion.txt')

fd = FreqDist(all_words)

## initialize 2 empty lists to hold our ranks and fequencies
ranks = []
freqs = []

for rank, freq in enumerate([freq for word,freq in fd.most_common()]):
    ranks.append(rank+1)
    freqs.append(freq)

# ## Plot Ranks vs frequency
plt.loglog(ranks,freqs)
plt.xlabel('frequency(f)', fontsize=14, fontweight='bold')
plt.ylabel('rank(r)', fontsize=14, fontweight='bold')
plt.grid(True)
plt.show()

## note: make a ranked word list like this:
lst = list([word for word,freq in fd.most_common()])

max_stars_repo_name: jclosure/donkus

max_issues_repo_name: jclosure/donkus

max_forks_repo_name: jclosure/donkus

reason: from scipy


text:
"""
Created on Jan. 04, 2021
@author: Heng-Sheng (Hanson) Chang
"""

import numpy as np
from numba import njit

from elastica._linalg import _batch_matvec, _batch_cross
from elastica._calculus import quadrature_kernel, difference_kernel
from elastica.external_forces import inplace_addition, NoForces

@njit(cache=True)
def _material_to_lab(director_collection, vectors):
    blocksize = vectors.shape[1]
    lab_frame_vectors = np.zeros((3, blocksize))
    for n in range(blocksize):
        for i in range(3):
            for j in range(3):
                lab_frame_vectors[i, n] += (
                    director_collection[j, i, n] * vectors[j, n]
                )
    return lab_frame_vectors

@njit(cache=True)
def _lab_to_material(director_collection, vectors):
    return _batch_matvec(director_collection, vectors)

@njit(cache=True)
def _internal_to_external_load(
    director_collection, kappa, tangents,
    rest_lengths, rest_voronoi_lengths,
    dilatation, voronoi_dilatation,
    internal_forces, internal_couples,
    external_forces, external_couples
    ):

    external_forces[:, :] = difference_kernel(
        _material_to_lab(director_collection, internal_forces)
        )

    external_couples[:, :] = (
        difference_kernel(internal_couples) +
        quadrature_kernel(
            _batch_cross(kappa, internal_couples) * rest_voronoi_lengths
        ) +
        _batch_cross(
                _lab_to_material(director_collection, tangents * dilatation),
                internal_forces
            )  * rest_lengths
    )
    return

@njit(cache=True)
def _force_induced_couple(internal_forces, distance, internal_couples):
    internal_couples[:, :] = quadrature_kernel(
        _batch_cross(distance, internal_forces)
    )[:, 1:-1]
    return

class ContinuousActuation(object):
    """
    This class and the classes ingerited from this should contain the parameters
    and functions that describe and calculate the forces/couples generated from
    the actuation.
    """

    def __init__(self, n_elements: int):
        self.internal_forces = np.zeros((3, n_elements))       # material frame
        self.external_forces = np.zeros((3, n_elements+1))     # global frame
        self.internal_couples = np.zeros((3, n_elements-1))    # material frame
        self.external_couples = np.zeros((3, n_elements))      # material frame

    def __call__(self, system):
        raise NotImplementedError

class ApplyActuation(NoForces):
    def __init__(self, actuation, step_skip: int, callback_params: dict):
        self.current_step = 0
        self.actuation = actuation
        self.every = step_skip
        self.callback_params = callback_params

    def apply_torques(self, system, time: np.float = 0.0):

        self.actuation(system)
        inplace_addition(system.external_forces, self.actuation.external_forces)
        inplace_addition(system.external_torques, self.actuation.external_couples)

        self.make_callback()

    def make_callback(self,):
        if self.current_step % self.every == 0:
            self.callback_func()
        self.current_step += 1

    def callback_func(self):
        self.callback_params['internal_forces'].append(self.actuation.internal_forces.copy())
        self.callback_params['internal_couples'].append(self.actuation.internal_couples.copy())
        self.callback_params['external_forces'].append(self.actuation.external_forces.copy())
        self.callback_params['external_couples'].append(self.actuation.external_couples.copy())

max_stars_repo_name: armantekinalp/gym-softrobot

max_issues_repo_name: armantekinalp/gym-softrobot

max_forks_repo_name: armantekinalp/gym-softrobot

reason: import numpy,from numba


text:
# -*- coding: utf-8 -*-
"""
Module to plot EISCAT radar scans. Scan detection based on azimuth (primary) or
elevation (if stationary in azimuth)

@author: Christer van der Meeren

Initially tested at the ESR 2014-11-26

"""

from __future__ import print_function
import os
import shutil
from os.path import isfile, join
import re
import logging
import pickle
import datetime as dt
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import patheffects as pe
from matplotlib.ticker import LogLocator, MultipleLocator
from eiscat_toolkit import loc2gg, load_param_single_simple
from helper_functions import find_coord_alt, find_coord_alt_2, bm_drawCmappedPoly, cbar_right
from mpl_toolkits.basemap import Basemap
from time import sleep

plt.ioff()  # performance impact is negligible, but when plotting many plots, this allows to plot and save without showing figure
np.seterr(invalid='ignore', divide='ignore')
plt.rcParams['font.size'] = 12

#==============================================================================
# Some defaults
#==============================================================================

# behaviour
RT = True  # realtime plotting. continuously checks data folder for new files
RT_replotAfterScan = True  # after a realtime scan/plot is finished, close and plot everything at once and save this plot instead. This corrects the flat-projection plots and doesn't take long.
onlyDoScanNo = None  # the only scan number to plot. Takes precedence over which scan to start at (which is given in input when running the code)

# plot-related
radarLoc = [78.153, 16.029, 0.438]  # latitude, longitude, alt [km] of radar (and center of map)
mapWidth = 1.8e6  # map width (and height) in map projection coordinates
figSize = 60  # DPI of figure (makes EVERYTHING larger/smaller)

# debug
debugRT = False  # forces update of plot after each file is read

#logging.basicConfig(level=logging.INFO)  # uncomment to turn on info messages (not debug)
#logging.basicConfig(level=logging.DEBUG)  # uncomment to turn on debug messages

#==============================================================================
# End switchboard
#==============================================================================

# get most recent 32m data directory
baseDataDir = '/analysis/results'  # will look for 32m data folders in this directory
if os.path.isdir(baseDataDir):
    subdirs_32m = [join(baseDataDir, d) for d in os.listdir(baseDataDir) if os.path.isdir(os.path.join(baseDataDir, d)) and '@32m' in d]
    dataFolder = join(baseDataDir, max(subdirs_32m, key=os.path.getmtime))
else:
    dataFolder = None

class ScanDetectionError(Exception):

    def __init__(self, value):
        self.value = value

    def __str__(self):
        return repr(self.value)

class Scan(object):

    def __init__(self, init_Time, init_par1D, init_par2D, init_err2D, scDir, scNo):

        # plot-related attributes
        self.fig = None
        self.axes = []
        self.map = None
        self.plotAlts = []  # contains the altitude lines to be plotted
        self.plottedAlts = {}
        self.plottedAltsFlat = {}
        self.altLats = {}
        self.altLons = {}
        self.altLatsFlat = {}
        self.altLonsFlat = {}
        self.cbar = [None]*4
        self.cbarAxis = [None]*4
        self.plottedBeams = []
        self.mainTitle = None
        self.byline = None
        self.elScanDirectionPlotted = False
        self.scDirElev = None
        self.towardAwayPlotted = False

        # data
        self.Time = init_Time
        self.par1D = init_par1D
        self.par2D = init_par2D
        self.err2D = init_err2D

        # scan direction and number
        self.scDir = scDir
        self.scNo = scNo

        # is scan finished? (Used for determining method for rotating data in flat-projection plots)
        self.finished = False

        self.midTimes = self.Time[0, :] + (self.Time[1, :] - self.Time[0, :]) / 2
        self.scanStart = self.Time[0, 0]
        self.scanEnd = self.Time[-1, -1]

        # convenience shortcuts
        self.tStart = self.Time[0, :]
        self.tEnd = self.Time[1, :]
        self.Az = self.par1D[:, 0]
        self.El = self.par1D[:, 1]
        self.Pt = self.par1D[:, 2]
        self.Tsys = self.par1D[:, 3]
        try:
            self.Oppd_Php = self.par1D[:, 4]
        except IndexError:
            self.Oppd_Php = None
        self.Ran = self.par2D[:, :, 0]
        self.Alt = self.par2D[:, :, 1]
        self.Ne = self.par2D[:, :, 2]
        self.Te = self.par2D[:, :, 3]
        self.Ti = self.par2D[:, :, 4]
        self.Vi = self.par2D[:, :, 5]
        self.Coll = self.par2D[:, :, 6]
        self.Comp = self.par2D[:, :, 7]
        self.Res = self.par2D[:, :, 8]
        self.errNe = self.err2D[:, :, 0]
        self.errTe = self.err2D[:, :, 1]
        self.errTi = self.err2D[:, :, 2]
        self.errVi = self.err2D[:, :, 3]
        self.errColl = self.err2D[:, :, 4]

    def add_data(self, new_Time, new_par1D, new_par2D, new_err2D):

        # correct dimensions if number of range gates have changed
        self.par2D, new_par2D, self.err2D, new_err2D = fix_dimensions(self.par2D, new_par2D, self.err2D, new_err2D)

        self.Time = np.append(self.Time, new_Time, axis=1)
        self.par1D = np.append(self.par1D, new_par1D, axis=0)
        self.par2D = np.append(self.par2D, new_par2D, axis=1)
        self.err2D = np.append(self.err2D, new_err2D, axis=1)

        self.midTimes = self.Time[0, :] + (self.Time[1, :] - self.Time[0, :]) / 2
        self.scanEnd = self.Time[-1, -1]

        self.tStart = self.Time[0, :]
        self.tEnd = self.Time[1, :]
        self.Az = self.par1D[:, 0]
        self.El = self.par1D[:, 1]
        self.Pt = self.par1D[:, 2]
        self.Tsys = self.par1D[:, 3]
        try:
            self.Oppd_Php = self.par1D[:, 4]
        except IndexError:
            self.Oppd_Php = None
        self.Ran = self.par2D[:, :, 0]
        self.Alt = self.par2D[:, :, 1]
        self.Ne = self.par2D[:, :, 2]
        self.Te = self.par2D[:, :, 3]
        self.Ti = self.par2D[:, :, 4]
        self.Vi = self.par2D[:, :, 5]
        self.Coll = self.par2D[:, :, 6]
        self.Comp = self.par2D[:, :, 7]
        self.Res = self.par2D[:, :, 8]
        self.errNe = self.err2D[:, :, 0]
        self.errTe = self.err2D[:, :, 1]
        self.errTi = self.err2D[:, :, 2]
        self.errVi = self.err2D[:, :, 3]
        self.errColl = self.err2D[:, :, 4]

    def _find_az_speed(self, ind1, ind2):
        """Returns the azimuth speed in deg/sec between the beams corresponding
        to the two indices (scan direction from ind1 to ind2)"""
        azimDist = self.Az[ind2]-self.Az[ind1]
        timeDist = self.tStart[ind2]-self.tStart[ind1]
        return azimDist/timeDist.total_seconds()

    def _find_el_speed(self, ind1, ind2):
        """Returns the azimuth speed in deg/sec between the beams corresponding
        to the two indices (scan direction from ind1 to ind2)"""
        elevDist = self.El[ind2]-self.El[ind1]
        timeDist = self.tStart[ind2]-self.tStart[ind1]
        return elevDist/timeDist.total_seconds()

    def _data_skip_before(self, ind):
        """Returns True if there is a data skip before the beam corresponding
        to ind (as determined by the difference between end time of ind-1 and
        the start time of ind)"""
        if ind == 0 or ind == len(self.Az):
            # data skip before first beam and after last beam
            return True
        else:
            # currently the limit is set at 2 seconds between end[i-1] and start[i]
            return abs(self.tStart[ind]-self.tEnd[ind-1]).total_seconds() >= 2

    def _data_skip_after(self, ind):
        """Same as _data_skip_before, but checks for data skip AFTER scan"""
        return self._data_skip_before(ind+1)

    def _vertex_array(self, alts, radarLoc, doAll=False):
        """Returns an array of vertices to be used with :function:`~cvdm.bm_drawCmappedPoly`"""

        # skip already plotted beams
        allBeams = np.array(range(0, len(self.Az)))
        plottedBeams = np.array(self.plottedBeams)
        if doAll:
            self.newBeamNos = allBeams
            self.altLats = {}
            self.altLons = {}
        else:
            self.newBeamNos = np.setdiff1d(allBeams, plottedBeams)

        # initiate data structures
        vertexLats = np.ones((self.Ran.shape[0] * len(self.newBeamNos), 4))
        vertexLons = np.ones((self.Ran.shape[0] * len(self.newBeamNos), 4))
        vertexLats_flat = np.ones((self.Ran.shape[0] * len(self.newBeamNos), 4))
        vertexLons_flat = np.ones((self.Ran.shape[0] * len(self.newBeamNos), 4))
        vertexX_elev = np.ones((self.Ran.shape[0] * len(self.newBeamNos), 4))
        vertexY_elev = np.ones((self.Ran.shape[0] * len(self.newBeamNos), 4))

        p = 0  # polygon counter

        for beamNo in self.newBeamNos:

            # find beam width (az1 and az2)
            if not self._data_skip_before(beamNo) and not self._data_skip_after(beamNo):  # not including first/last beams
                az1 = (self.Az[beamNo] + self.Az[beamNo-1]) / 2
                az2 = (self.Az[beamNo] + self.Az[beamNo+1]) / 2
            elif self._data_skip_before(beamNo) and not self._data_skip_after(beamNo):  # including first beam
                az2 = (self.Az[beamNo] + self.Az[beamNo+1]) / 2
                az1 = self.Az[beamNo] - (az2 - self.Az[beamNo])
            elif not self._data_skip_before(beamNo) and self._data_skip_after(beamNo):  # including last beam
                az1 = (self.Az[beamNo] + self.Az[beamNo-1]) / 2
                az2 = self.Az[beamNo] + (self.Az[beamNo] - az1)
            else:  # data skip both before and after
                azSpeed = self._find_az_speed(beamNo-1, beamNo)
                beamWidth = abs(azSpeed * (self.tEnd[beamNo] - self.tStart[beamNo]).total_seconds())
                az1 = self.Az[beamNo] - beamWidth/2
                az2 = self.Az[beamNo] + beamWidth/2

            # do the same for elevation
            if not self._data_skip_before(beamNo) and not self._data_skip_after(beamNo):  # not including first/last beams
                el1 = (self.El[beamNo] + self.El[beamNo-1]) / 2
                el2 = (self.El[beamNo] + self.El[beamNo+1]) / 2
            elif self._data_skip_before(beamNo) and not self._data_skip_after(beamNo):  # including first beam
                el2 = (self.El[beamNo] + self.El[beamNo+1]) / 2
                el1 = self.El[beamNo] - (el2 - self.El[beamNo])
            elif not self._data_skip_before(beamNo) and self._data_skip_after(beamNo):  # including last beam
                el1 = (self.El[beamNo] + self.El[beamNo-1]) / 2
                el2 = self.El[beamNo] + (self.El[beamNo] - el1)
            else:  # data skip both before and after
                elSpeed = self._find_el_speed(beamNo-1, beamNo)
                beamWidth = abs(elSpeed * (self.tEnd[beamNo] - self.tStart[beamNo]).total_seconds())
                el1 = self.El[beamNo] - beamWidth/2
                el2 = self.El[beamNo] + beamWidth/2

            # get coordinates for altitude lines
            # TODO: Find better solution than below for avoiding "stair lines"?
            if beamNo == 0:
                azList = [az1, az2]
            else:
                azList = [az2]
            for az in azList:
                for alt in alts:
                    lat, lon, h = find_coord_alt(radarLoc, self.El[beamNo], az, alt)
                    azFlat = az if self.El[beamNo] <= 90 else az + 180
                    latFlat, lonFlat, hFlat = find_coord_alt_2(radarLoc, self.El[beamNo], 30, azFlat, alt)
                    try:  # append coordinates to this altitude
                        self.altLats[alt] = np.append(self.altLats[alt], lat)
                        self.altLons[alt] = np.append(self.altLons[alt], lon)
                        self.altLatsFlat[alt] = np.append(self.altLatsFlat[alt], latFlat)
                        self.altLonsFlat[alt] = np.append(self.altLonsFlat[alt], lonFlat)
                    except KeyError:  # first values for this altitude
                        self.altLats[alt] = np.array(lat)
                        self.altLons[alt] = np.array(lon)
                        self.altLatsFlat[alt] = np.array(latFlat)
                        self.altLonsFlat[alt] = np.array(lonFlat)

            # loop over all range gates within this beam
            beamRans = self.Ran[:, beamNo]  # ranges for this beam
            for ranNo, ran in enumerate(beamRans):

                if ranNo == 0:  # first range gate
                    ran2 = (ran + beamRans[ranNo+1]) / 2
                    ran1 = ran - (ran2 - ran)
                elif ranNo == len(beamRans)-1:  # last range gate
                    ran1 = (ran + beamRans[ranNo-1]) / 2
                    ran2 = ran + (ran - ran1)
                else:  # any beam not first or last
                    ran1 = (ran + beamRans[ranNo-1]) / 2
                    ran2 = (ran + beamRans[ranNo+1]) / 2

                # compute lats and lons of vertices and add to structure
                if not np.isnan(ran1) and not np.isnan(ran2):
                    v1 = loc2gg(radarLoc, [self.El[beamNo], az1, ran1])
                    v2 = loc2gg(radarLoc, [self.El[beamNo], az1, ran2])
                    v3 = loc2gg(radarLoc, [self.El[beamNo], az2, ran2])
                    v4 = loc2gg(radarLoc, [self.El[beamNo], az2, ran1])
                    vertexLats[p, :] = [v1[0], v2[0], v3[0], v4[0]]
                    vertexLons[p, :] = [v1[1], v2[1], v3[1], v4[1]]

                    azFlat1 = az1 if self.El[beamNo] <= 90 else az1 + 180
                    azFlat2 = az2 if self.El[beamNo] <= 90 else az2 + 180
                    v1_flat = loc2gg(radarLoc, [30, azFlat1, ran1])
                    v2_flat = loc2gg(radarLoc, [30, azFlat1, ran2])
                    v3_flat = loc2gg(radarLoc, [30, azFlat2, ran2])
                    v4_flat = loc2gg(radarLoc, [30, azFlat2, ran1])
                    vertexLats_flat[p, :] = [v1_flat[0], v2_flat[0], v3_flat[0], v4_flat[0]]
                    vertexLons_flat[p, :] = [v1_flat[1], v2_flat[1], v3_flat[1], v4_flat[1]]

                    v1_elev = (np.cos(np.deg2rad(el1))*ran1, np.sin(np.deg2rad(el1))*ran1)
                    v2_elev = (np.cos(np.deg2rad(el1))*ran2, np.sin(np.deg2rad(el1))*ran2)
                    v3_elev = (np.cos(np.deg2rad(el2))*ran2, np.sin(np.deg2rad(el2))*ran2)
                    v4_elev = (np.cos(np.deg2rad(el2))*ran1, np.sin(np.deg2rad(el2))*ran1)
                    vertexX_elev[p, :] = [v1_elev[0], v2_elev[0], v3_elev[0], v4_elev[0]]
                    vertexY_elev[p, :] = [v1_elev[1], v2_elev[1], v3_elev[1], v4_elev[1]]
                else:
                    vertexLats[p, :] = [np.nan]*4
                    vertexLons[p, :] = [np.nan]*4
                    vertexLats_flat[p, :] = [np.nan]*4
                    vertexLons_flat[p, :] = [np.nan]*4
                    vertexX_elev[p, :] = [np.nan]*4
                    vertexY_elev[p, :] = [np.nan]*4

                p = p + 1

            self.plottedBeams.append(beamNo)

        return vertexLats, vertexLons, vertexLats_flat, vertexLons_flat, vertexX_elev, vertexY_elev

    def plot_alts(self, plotFlat=False, rotFlat=None, mapObj=None, ax=None):
        """Plots altitude lines for all altitude in the list `alts`"""

        mapObj = self.map if mapObj is None else mapObj

        for alt in self.altLats:

            # transform to map coordinates
            if hasattr(mapObj, 'coords'):
                x, y = mapObj(self.altLons[alt], self.altLats[alt], coords='geo')
                xFlat, yFlat = mapObj(self.altLonsFlat[alt], self.altLatsFlat[alt], coords='geo')
            else:
                x, y = mapObj(self.altLons[alt], self.altLats[alt])
                xFlat, yFlat = mapObj(self.altLonsFlat[alt], self.altLatsFlat[alt])

            # perform rotation for flat-projection plots if given
            if rotFlat is not None:
                angle = np.deg2rad(rotFlat[0])
                rotMatrix = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]])
                # shift rotation point to origin
                if hasattr(mapObj, 'coords'):
                    x0, y0 = mapObj(rotFlat[2], rotFlat[1], coords='geo')
                else:
                    x0, y0 = mapObj(rotFlat[2], rotFlat[1])
                xFlat = xFlat - x0
                yFlat = yFlat - y0

                # rotate
                for i, xy in enumerate(zip(xFlat, yFlat)):
                    xFlat[i] = np.dot(xy, rotMatrix)[0]
                    yFlat[i] = np.dot(xy, rotMatrix)[1]

                # shift back
                xFlat = xFlat + x0
                yFlat = yFlat + y0

            # text alignment for main plots
            xMid = (mapObj.urcrnrx - mapObj.llcrnrx) / 2
            yMid = (mapObj.urcrnry - mapObj.llcrnry) / 2

            if (x[0] >= xMid and y[0] <= yMid and self.scDir == 'cw') or (x[0] <= xMid and y[0] >= yMid and self.scDir == 'ccw'):
                # lower right quadrant and clockwise scan, or upper left quadrant and ccw
                hor = 'left'
                vert = 'bottom'
#                text_alignment = (0, 0)
            elif (x[0] <= xMid and y[0] <= yMid and self.scDir == 'cw') or (x[0] >= xMid and y[0] >= yMid and self.scDir == 'ccw'):
                # lower left quadrant and clockwise scan, or upper right quadrant and ccw scan
                hor = 'left'
                vert = 'top'
#                text_alignment = (0, 1)
            elif (x[0] <= xMid and y[0] >= yMid and self.scDir == 'cw') or (x[0] >= xMid and y[0] <= yMid and self.scDir == 'ccw'):
                # upper left quadrant and clockwise scan, or lower right quadrant and ccw scan
                hor = 'right'
                vert = 'top'
#                text_alignment = (1, 1)
            else:  # upper right quadrant
                hor = 'right'
                vert = 'bottom'
#                text_alignment = (1, 0)

            # text alignment for flat-projection plots
            if self.scDir == 'cw':
                horFlat = 'center'
                vertFlat = 'top'
#                text_alignment_Flat = (0.5, 1)
                angleFlat = -30
                sFlat = str(alt)
            else:
                horFlat = 'center'
                vertFlat = 'top'
#                text_alignment_Flat = (0.5, 1)
                angleFlat = 30
                sFlat = str(alt)

            # initiate dict item for this alt
            try:
                self.plottedAlts[alt]
            except KeyError:
                self.plottedAlts[alt] = []
                self.plottedAltsFlat[alt] = []

            if ax is None:  # plotting by self.plot()

                for _ax in self.axes[0:4]:
                    self.plottedAlts[alt].append(mapObj.plot(x, y, linewidth=2.5, color='w', ax=_ax)[0])
                    self.plottedAlts[alt].append(mapObj.plot(x, y, linewidth=1.5, color='k', ax=_ax)[0])
                    self.plottedAlts[alt].append(_ax.text(x[0], y[0], s=' ' + str(alt) + '$\,$km ', ha=hor, va=vert, path_effects=[pe.withStroke(linewidth=3, foreground='w')]))

                for _ax in self.axes[4:8]:
                    if self.scDir not in ['elev incr', 'elev decr']:  # don't plot in flat-projection for elevation scans
                        self.plottedAltsFlat[alt].append(mapObj.plot(xFlat, yFlat, linewidth=2.5, color='w', ax=_ax)[0])
                        self.plottedAltsFlat[alt].append(mapObj.plot(xFlat, yFlat, linewidth=1.5, color='k', ax=_ax)[0])
                        self.plottedAltsFlat[alt].append(_ax.text(xFlat[0], yFlat[0], s=' ' + str(alt) + '$\,$km ', ha=hor, va=vert, path_effects=[pe.withStroke(linewidth=3, foreground='w')]))

            else:  # plotting by self.plot_overlay() or something else
                if plotFlat:
                    self.plottedAltsFlat[alt].append(mapObj.plot(xFlat, yFlat, linewidth=2.5, color='w', ax=ax)[0])
                    self.plottedAltsFlat[alt].append(mapObj.plot(xFlat, yFlat, linewidth=1.5, color='k', ax=ax)[0])
                    self.plottedAltsFlat[alt].append(ax.text(xFlat[0], yFlat[0], s=sFlat, horizontalalignment=horFlat, verticalalignment=vertFlat, rotation=angleFlat, path_effects=[pe.withStroke(linewidth=3, foreground='w')]))
                else:
                    self.plottedAlts[alt].append(mapObj.plot(x, y, linewidth=2.5, color='w', ax=ax)[0])
                    self.plottedAlts[alt].append(mapObj.plot(x, y, linewidth=1.5, color='k', ax=ax)[0])
                    self.plottedAlts[alt].append(ax.text(x[0], y[0], s=str(alt) + '$\,$km', horizontalalignment=hor, verticalalignment=vert, path_effects=[pe.withStroke(linewidth=3, foreground='w')]))

    def update_alts(self, rotFlat=None):

        # lines in large subplots
        for alt, lines in self.plottedAlts.items():
            for line in lines:
                try:
                    line.set_data(*self.map(self.altLons[alt], self.altLats[alt]))
                except AttributeError:  # Probably a Text object
                    pass

        # lines in small subplots
        for alt, lines in self.plottedAltsFlat.items():

            xFlat, yFlat = self.map(self.altLonsFlat[alt], self.altLatsFlat[alt])

            # perform rotation for flat-projection plots if given
            if rotFlat is not None:
                angle = np.deg2rad(rotFlat[0])
                rotMatrix = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]])
                # shift rotation point to origin
                x0, y0 = self.map(rotFlat[2], rotFlat[1])
                xFlat = xFlat - x0
                yFlat = yFlat - y0

                # rotate
                for i, xy in enumerate(zip(xFlat, yFlat)):
                    xFlat[i] = np.dot(xy, rotMatrix)[0]
                    yFlat[i] = np.dot(xy, rotMatrix)[1]

                # shift back
                xFlat = xFlat + x0
                yFlat = yFlat + y0

            for line in lines:
                try:
                    line.set_data(xFlat, yFlat)
                except AttributeError:  # Probably a Text object
                    pass

    def plot(self, drawMag=False, removeLargeErrs=False):

        # for some reason, plotting only two beams doesn't work...
        if len(self.Az) < 3:
            return

        # initiate figure/axes
        if self.fig is None:

            self.fig = plt.figure(figsize=(22, 14), dpi=figSize)

            gs_main = mpl.gridspec.GridSpec(3, 1, height_ratios=[2, 2, 0.8])
            gs_map = mpl.gridspec.GridSpecFromSubplotSpec(1, 4, subplot_spec=gs_main[0])
            gs_flat = mpl.gridspec.GridSpecFromSubplotSpec(1, 4, subplot_spec=gs_main[1])
            gs_elev = mpl.gridspec.GridSpecFromSubplotSpec(1, 4, subplot_spec=gs_main[2])

            self.axes.append(plt.subplot(gs_map[0], aspect='equal', adjustable='box-forced'))
            self.axes.append(plt.subplot(gs_map[1], aspect='equal', adjustable='box-forced', sharex=self.axes[0], sharey=self.axes[0]))
            self.axes.append(plt.subplot(gs_map[2], aspect='equal', adjustable='box-forced', sharex=self.axes[0], sharey=self.axes[0]))
            self.axes.append(plt.subplot(gs_map[3], aspect='equal', adjustable='box-forced', sharex=self.axes[0], sharey=self.axes[0]))
            self.axes.append(plt.subplot(gs_flat[0], aspect='equal', adjustable='box-forced'))
            self.axes.append(plt.subplot(gs_flat[1], aspect='equal', adjustable='box-forced', sharex=self.axes[4], sharey=self.axes[4]))
            self.axes.append(plt.subplot(gs_flat[2], aspect='equal', adjustable='box-forced', sharex=self.axes[4], sharey=self.axes[4]))
            self.axes.append(plt.subplot(gs_flat[3], aspect='equal', adjustable='box-forced', sharex=self.axes[4], sharey=self.axes[4]))
            self.axes.append(plt.subplot(gs_elev[0], aspect='equal', adjustable='box-forced'))
            self.axes.append(plt.subplot(gs_elev[1], aspect='equal', adjustable='box-forced', sharex=self.axes[8], sharey=self.axes[8]))
            self.axes.append(plt.subplot(gs_elev[2], aspect='equal', adjustable='box-forced', sharex=self.axes[8], sharey=self.axes[8]))
            self.axes.append(plt.subplot(gs_elev[3], aspect='equal', adjustable='box-forced', sharex=self.axes[8], sharey=self.axes[8]))

            gs_main.update(left=0.035, right=0.96, bottom=0.06, top=0.95, wspace=0.1, hspace=0)

            # initiate maps
            self.map = Basemap(width=mapWidth,
                               height=mapWidth,
                               projection='aeqd',
                               lat_0=radarLoc[0],
                               lon_0=radarLoc[1],
                               resolution='l')

            # draw coastlines
            for ax in self.axes[0:4]:
                self.map.drawcoastlines(linewidth=0.5, color="k", ax=ax)

            # draw magnetic meridian/parallel through ESR
            if drawMag:
                with open('magLatLon.pckl', 'r') as f:
                    magLatLon = pickle.load(f)
                for ax in self.axes[0:4]:
                    self.map.plot(magLatLon[1], magLatLon[0], latlon=True, ax=ax, color='k', linestyle='dotted', linewidth=3, path_effects=[pe.withStroke(foreground='w', linewidth=5)])
                    self.map.plot(magLatLon[3], magLatLon[2], latlon=True, ax=ax, color='k', linestyle='dotted', linewidth=3, path_effects=[pe.withStroke(foreground='w', linewidth=5)])

            # draw inivible coastlines in flat plots
            for ax in self.axes[4:8]:
                self.map.drawcoastlines(linewidth=0, color="w", zorder=-100, ax=ax)

            # small plot "titles"
            for ax in self.axes[4:8]:
                ax.set_title(u'Flattened to 30\u00B0 elev\n' +
                             u'(useful e.g. for for mixed az/el scans)')

        # get coordinates of things to plot
        vertexLats, vertexLons, vertexLats_flat, vertexLons_flat, vertexX_elev, vertexY_elev = self._vertex_array(self.plotAlts, radarLoc)

        # plotting parameters
        #          data   clims          ticks                           colormap       logarithmic colormap
#        toPlot = [('Ne',  (0, 1e12),     MultipleLocator(base=2e11),     'jet',         False),
        toPlot = [('Ne',  (1e10, 1e12),  LogLocator(subs=range(1, 10)),  'nipy_spectral_pinktop',         True),
                  ('Vi',  (-500, 500),   [-500, -250, 0, 250, 500],      'RdBu_r',  False),
                  ('Te',  (0, 4000),     [0, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000],    'nipy_spectral_pinktop',         False),
                  ('Ti',  (0, 3000),     [0, 500, 1000, 1500, 2000, 2500, 3000],          'nipy_spectral_pinktop',         False)]

        for i, e in enumerate(toPlot):

            data = np.ma.masked_invalid(eval("self." + e[0])[:, self.newBeamNos])  # mask nans (and infs)
            if removeLargeErrs:
                data = np.ma.masked_where(eval("self.err" + e[0])[:, self.newBeamNos] > abs(eval("self." + e[0])[:, self.newBeamNos]), data)  # mask where error > value
            data = data.flatten(order='F')

            # main plots
            pCol = bm_drawCmappedPoly(self.map,
                                      vertexLons,
                                      vertexLats,
                                      data,
                                      clims=e[1],
                                      cmap=e[3],
                                      logMap=e[4],
                                      linewidths=0)
            self.axes[i].add_collection(pCol)
            if None in self.cbar:
                self.cbar[i], self.cbarAxis[i] = cbar_right(self.axes[i], leftOffset=0.003, width=0.01, relHeight=1, mappable=pCol, ticks=e[2])

            # flat-projection plots
            if self.scDir in ['cw', 'ccw']:
                pCol = bm_drawCmappedPoly(self.map,
                                          vertexLons_flat,
                                          vertexLats_flat,
                                          data,
                                          clims=e[1],
                                          cmap=e[3],
                                          logMap=e[4],
                                          linewidths=0)
                self.axes[i+4].add_collection(pCol)

            # elev plots
            self.axes[i+8].set_xlim((-800, 800))
            self.axes[i+8].set_ylim((100, 700))
            self.axes[i+8].set_yticks([200, 400, 600])
            self.axes[i+8].set_xlabel('Flat ground distance [km]')
            self.axes[i+8].set_ylabel('Altitude [km]')
            self.axes[i+8].grid('on')
            self.axes[i+8].set_title('Elevation-only scan')

            if self.scDir in ['elev incr', 'elev decr']:
                verts = map(zip, vertexX_elev, vertexY_elev)
                if e[4]:
                    norm = mpl.colors.LogNorm(vmin=e[1][0], vmax=e[1][1])
                else:
                    norm = mpl.colors.Normalize(vmin=e[1][0], vmax=e[1][1])
                pCol = mpl.collections.PolyCollection(verts, array=data, cmap=e[3], norm=norm, linewidths=0)
                self.axes[i+8].add_collection(pCol)

                # show cardinal directions on elevation scans
                if not self.elScanDirectionPlotted:
                    azDir = self.Az[0]
                    if self.scDir == 'elev decr':
                        azDir += 180

                    # function to find out if azDir is within 22.5 degrees of some direction
                    az_near = lambda x: np.cos(np.deg2rad(azDir-x)) >= np.cos(np.deg2rad(22.5))

                    # dictionary of direction string for positive/negative elevation at a given azimuth
                    dirDict = {0: ['N', 'S'],
                               45: ['NE', 'SW'],
                               90: ['E', 'W'],
                               135: ['SE', 'NW'],
                               180: ['S', 'N'],
                               225: ['SW', 'NE'],
                               270: ['W', 'E'],
                               315: ['NW', 'SE']}

                    # find the correct direction string
                    for direction in dirDict:
                        if az_near(direction):
                            dirs = dirDict[direction]
                            self.scDirElev = '-'.join(dirs)
                            break

                    left, right = dirs if self.scDir == 'elev decr' else dirs[::-1]
                    self.axes[i+8].text(0, -0.25, left, ha='left', fontweight='bold', transform=self.axes[i+8].transAxes)
                    self.axes[i+8].text(1, -0.25, right, ha='right', fontweight='bold', transform=self.axes[i+8].transAxes)

        self.elScanDirectionPlotted = True

        # toward/away on velocity cbar
        if not self.towardAwayPlotted:
            self.cbarAxis[1].text(0, 1.03, 'away', ha='left', va='bottom', transform=self.cbarAxis[1].transAxes)
            self.cbarAxis[1].text(0, -0.03, 'toward', ha='left', va='top', transform=self.cbarAxis[1].transAxes)
            self.towardAwayPlotted = True

        # big plot titles
        self.axes[0].set_title(u'Electron density [m$\mathregular{^{−3}}$]')
        self.axes[1].set_title('Ion velocity [m/s]')
        self.axes[2].set_title('Electron temperature [K]')
        self.axes[3].set_title('Ion temperature [K]')

        for ax in [self.axes[0], self.axes[4], self.axes[8]]:
            ax.text(0, 1, 'Ne', transform=ax.transAxes, ha='left', va='bottom', fontweight='bold')
        for ax in [self.axes[1], self.axes[5], self.axes[9]]:
            ax.text(0, 1, 'Vi', transform=ax.transAxes, ha='left', va='bottom', fontweight='bold')
        for ax in [self.axes[2], self.axes[6], self.axes[10]]:
            ax.text(0, 1, 'Te', transform=ax.transAxes, ha='left', va='bottom', fontweight='bold')
        for ax in [self.axes[3], self.axes[7], self.axes[11]]:
            ax.text(0, 1, 'Ti', transform=ax.transAxes, ha='left', va='bottom', fontweight='bold')

        # plot altitude levels
        if self.plottedAlts:
            self.update_alts()
        else:
            self.plot_alts()

        # set main title of plot
        mainTitleStr = u'Scan #{0} ({1})   {2}\u2212{3} UT'.format(self.scNo, self.scDirElev or self.scDir, self.tStart[0].strftime('%d %b %Y %H:%M:%S'), self.tEnd[-1].strftime('%H:%M:%S'))
        if self.mainTitle is None:
            self.mainTitle = self.fig.text(0.5, 0.98, mainTitleStr, horizontalalignment='center', verticalalignment='top', fontsize='xx-large')
        else:
            self.mainTitle.set_text(mainTitleStr)

        if self.byline is None:
            self.byline = self.fig.text(0.98, 0.98, 'Plotting software by Christer van der Meeren (cmeeren@gmail.com)\nCode available at https://github.com/cmeeren/eiscatscanplot', ha='right', va='top')

        # draw/update plot if realtime or only plotting a single scan
        if ((RT or debugRT) and not self.finished) or onlyDoScanNo is not None:
            plt.pause(0.001)

    def plot_overlay(self, mapObj, ax=None, data='Ne', altLines=None, radarLoc=None, removeLargeErrs=False, clims=None, cmap=None, logMap=None, flatProjection=False, defScanSpeedPerIP=1.0*3.2, **kwargs):

        altLines = altLines or []
        radarLoc = radarLoc or [78.153, 16.029, 0.438]

        # input handling: data
        if data is 'Ne':
            clims = [1e10, 1e12] if clims is None else clims
            cmap = 'jet' if cmap is None else cmap
            logMap = True if logMap is None else logMap
        elif data is 'Vi':
            clims = [-500, 500] if clims is None else clims
            cmap = 'esrBlueRed' if cmap is None else cmap
            logMap = False if logMap is None else logMap
        elif data is 'Te':
            clims = [0, 4000] if clims is None else clims
            cmap = 'jet' if cmap is None else cmap
            logMap = False if logMap is None else logMap
        elif data is 'Ti':
            clims = [0, 3000] if clims is None else clims
            cmap = 'jet' if cmap is None else cmap
            logMap = False if logMap is None else logMap
        else:
            print('Invalid data specification:', data)
            return

        # input handling: ax
        ax = plt.gca() if ax is None else ax

        # get (ccw) rotation angle for flat-projection plots (ideally, zenith should be vertical)
        ind_maxEl = np.argmin(np.abs(self.El-90))
        # ind_maxEl corresponds to the last data point of the upleg or the first data point of the downleg
        if np.abs(self.El[ind_maxEl]-90) < 30:
            if np.abs(self.El[ind_maxEl]-90) < defScanSpeedPerIP or self._data_skip_before(ind_maxEl) or self._data_skip_after(ind_maxEl):
                # zenith dump is close enough to 90 deg (or the data skips), plot should be centered at 90 degrees
                # use congruence of scan to find azimuth where el=90
                # TODO: What if ind_maxEl corresponds to one of first two indices?
                az1 = self.Az[ind_maxEl-2]
                az2 = self.Az[ind_maxEl-1]
                el1 = self.El[ind_maxEl-2]
                el2 = self.El[ind_maxEl-1]
                smallPlotRotAngle = az1 + ((90 - el1) * (az2 - az1)) / (el2 - el1)
            else:
                # plot should be centered not at 90 degrees, but at the highest elevation
                smallPlotRotAngle = self.Az[ind_maxEl]
#            if self.El[ind_maxEl-1] < self.El[ind_maxEl]:  # upleg
#                smallPlotRotAngle = self.Az[ind_maxEl]+(self.El[ind_maxEl]-90)
#            else:
#                smallPlotRotAngle = self.Az[ind_maxEl-1]+(self.El[ind_maxEl-1]-90)
        else:
            # max elevation of scan not close to 90, make symmetric instead
            smallPlotRotAngle = (self.Az[0] + self.Az[-1]) / 2

        if flatProjection and hasattr(mapObj, 'coords') and mapObj.coords != 'geo':
            print('WARNING: Flat-projected ESR plot will only work on geographic axes')

        # get coordinates of patch vertices
        vertexLats, vertexLons, vertexLats_flat, vertexLons_flat, vertexX_elev, vertexY_elev = self._vertex_array(altLines, radarLoc, doAll=True)

        if flatProjection:
            lats = vertexLats_flat
            lons = vertexLons_flat
            rot = [-smallPlotRotAngle, radarLoc[0], radarLoc[1]]
        else:
            lats = vertexLats
            lons = vertexLons
            rot = None

        # get data to plot
        data = np.ma.masked_invalid(eval("self." + data)[:, self.newBeamNos])  # mask nans (and infs)

        # mask where error > value, if specified
        if removeLargeErrs:
            data = np.ma.masked_where(eval("self.err" + data)[:, self.newBeamNos] > abs(eval("self." + data)[:, self.newBeamNos]), data)

        # flatten data to be used with custom function below
        data = data.flatten(order='F')

        # plot
        pCol = bm_drawCmappedPoly(mapObj,
                                  lons,
                                  lats,
                                  data,
                                  clims=clims,
                                  cmap=cmap,
                                  logMap=logMap,
                                  rot=rot,
                                  **kwargs)
        ax.add_collection(pCol)

        if flatProjection:
            self.plot_alts(mapObj=mapObj, ax=ax, plotFlat=True, rotFlat=rot)
        else:
            self.plot_alts(mapObj=mapObj, ax=ax)

        return pCol

    def saveFig(self, webAccessFolder='', webAccessFolderExternal=''):

        if self.fig is not None:
            fn = '{0}-{1}_{2}_scan{3:03d}.png'.format(self.tStart[0].strftime('%Y-%m-%d_%H%M'), self.tEnd[-1].strftime('%H%M'), self.scDirElev or self.scDir, self.scNo)
            saveTo = os.path.abspath(os.path.expanduser(savePath))
            if not os.path.isdir(saveTo):
                os.makedirs(saveTo)
            plt.savefig(os.path.join(saveTo, fn), dpi=figSize)
            if RT or debugRT:
                print('Saved file ' + os.path.join(saveTo, fn))
#                logging.info('Saved file ' + os.path.join(saveTo, fn))

            # make web page for external access
            if webAccessFolder:
                make_html(webAccessFolder, webAccessFolderExternal, saveTo, lastFile=fn)

        else:
            logging.warn('Scan #{}: No figure to save (perhaps because number of beams in scan is < 3)'.format(self.scNo))

    def closeFig(self):

        if self.fig is not None:
            logging.info('Closing figure')
            plt.close(self.fig)
        else:
            logging.info('Scan #{}: No figure to close (perhaps because number of beams in scan is < 3)'.format(self.scNo))

        # plot-related attributes
        self.fig = None
        self.axes = []
        self.map = None
        self.plottedAlts = {}
        self.plottedAltsFlat = {}
        self.altLats = {}
        self.altLons = {}
        self.altLatsFlat = {}
        self.altLonsFlat = {}
        self.cbar = [None]*4
        self.cbarAxis = [None]*4
        self.plottedBeams = []
        self.mainTitle = None
        self.byline = None
        self.elScanDirectionPlotted = False
        self.towardAwayPlotted = False

def scan_dir(currentAzim, lastAzim):
    '''Returns the direction of the scan (cw, ccw or stat) given two azimuths.'''
    if currentAzim is None or lastAzim is None:
        return None
    elif currentAzim > lastAzim:  # azimuth increases (not crossing 360/0 degrees)
        return 'cw'
    elif currentAzim < lastAzim:  # azimuth decreases (not crossing 360/0 degrees)
        return 'ccw'
    elif currentAzim == lastAzim:  # stationary
        return 'stat'

def scan_dir_elev(currentElev, lastElev):
    '''Returns the direction of the scan (cw, ccw or stat) given two azimuths.'''
    if currentElev is None or lastElev is None:
        return None
    elif currentElev > lastElev:  # azimuth increases (not crossing 360/0 degrees)
        return 'elev incr'
    elif currentElev < lastElev:  # azimuth decreases (not crossing 360/0 degrees)
        return 'elev decr'
    elif currentElev == lastElev:  # stationary
        return 'stat'

def fix_dimensions(par2D, new_par2D, err2D, new_err2D):
    """ correct dimensions if number of range gates have changed"""

    if par2D[:, 0, 0].shape[0] < new_par2D[:, 0, 0].shape[0]:
        par2D = np.append(par2D, np.ones((new_par2D.shape[0]-par2D.shape[0], par2D.shape[1], par2D.shape[2]))*np.nan, axis=0)
        err2D = np.append(err2D, np.ones((new_err2D.shape[0]-err2D.shape[0], err2D.shape[1], err2D.shape[2]))*np.nan, axis=0)
    if par2D[:, 0, 0].shape[0] > new_par2D[:, 0, 0].shape[0]:
        new_par2D = np.append(new_par2D, np.ones((err2D.shape[0]-new_par2D.shape[0], new_par2D.shape[1], new_par2D.shape[2]))*np.nan, axis=0)
        new_err2D = np.append(new_err2D, np.ones((err2D.shape[0]-new_err2D.shape[0], new_err2D.shape[1], new_err2D.shape[2]))*np.nan, axis=0)

    return par2D, new_par2D, err2D, new_err2D

def scan_parse(dataFolder, savePath,
               doPlot=False, onlyDoScanNo=None, startAt=None, removeLargeErrs=False, RT=False, RT_replotAfterScan=True,
               scanWidth=120, defScanSpeedPerIP=0.62*3, alts=None, drawMag=False, radarLoc=None, mapWidth=1.8e6, figSize=72,
               debugRT=False, webAccessFolder='', webAccessFolderExternal=''):
    '''docstring'''

    alts = [250, 500] if alts is None else alts
    radarLoc = radarLoc or [78.153, 16.029, 0.438]
    startAt = '1' if onlyDoScanNo else startAt

    IPerrorShown = False

    oldfiles = set([])  # already loaded data files
    currentScDir = None   # current scan direction
    lastScDir = None  # last scan direction
    currentAzim = None  # current azimuth
    lastAzim = None  # last azimuth
    currentScDirElev = None   # current scan direction in elevation
    lastScDirElev = None  # last scan direction in elevation
    currentElev = None  # current elevation
    lastElev = None  # last elevation

    # not really needed, but included to make the code editor stop complaining about using uninitiated variables
    Time_prev = None
    par1D_prev = None
    par2D_prev = None
    err2D_prev = None

    scanNoThis = 0

    # list to save scans to if not plotting
    allScans = []

    # find out which scan no / time we should start at
    if startAt is not None and ':' in startAt:
        startHour, startMinute = map(int, startAt.split(':'))
        startTime = dt.time(startHour, startMinute)
        startAtScanNo = 1
    elif startAt is not None:
        startAtScanNo = int(startAt)
        startTime = dt.time(0, 0)
    else:
        startAtScanNo = 1
        startTime = dt.time(0, 0)

    try:
        while True:
            allfiles = set([fn for fn in os.listdir(dataFolder) if (isfile(join(dataFolder, fn)) and re.match('\d*\.mat', fn))])
            newfiles = allfiles - oldfiles

            if len(allfiles) < 2:
                if RT:
                    print('Less than 2 files in folder, waiting 3 sec...')
                    sleep(3)
                else:
                    print('Less than 2 files in folder, exiting')
                    return
            elif len(newfiles) > 0:
                if RT:
                    print('\n{} new files discovered, analyzing...'.format(len(newfiles)))
                else:
                    print('Analyzing {} files...'.format(len(newfiles)))

                newfiles_iterator = sorted(list(newfiles))

                # try to do a progress bar if not realtime
                if not RT and not debugRT:
                    try:
                        import frogress
                        newfiles_iterator = frogress.bar(newfiles_iterator)
                    except:
                        pass

                ###########################################
                # LOOP OVER ALL FILES CURRENTLY IN FOLDER #
                ###########################################

                for fn in newfiles_iterator:

                    # load params
                    Time, par2D, par1D, rpar2D, err2D = load_param_single_simple(join(dataFolder, fn), trueAzEl=True)

                    # avoid crashing if a new beam somehow starts before the last beam
                    if Time_prev is not None and Time[0, 0] < Time_prev[0, 0]:
                        if not IPerrorShown:
                            raw_input('ERROR: Start timestamp of current integration is before start timestamp previous integration. This may happen if the GUISDAP analysis has been re-started using another integration period or start time. Things may not work properly from now. It is recommended that you delete the analysed datafiles and the plotted figures and start anew.\n\nPress Enter to continue >> ')
                            IPerrorShown = True
                        continue

                    # get current azim and scan direction
                    currentAzim = par1D[0, 0]
                    currentScDir = scan_dir(currentAzim, lastAzim)
                    currentElev = par1D[0, 1]
                    currentScDirElev = scan_dir_elev(currentElev, lastElev)
                    logging.debug('currentAzim = {} | lastAzim = {} | currentScDir = {} | lastScDir = {} | currentElev = {} | lastElev = {} | currentScDirElev = {} | lastScDirElev = {} | scanNoThis = {}'.format(currentAzim, lastAzim, currentScDir, lastScDir, currentElev, lastElev, currentScDirElev, lastScDirElev, scanNoThis))

                    if Time[0, 0].time() < startTime:
                        startAtScanNo = scanNoThis + 1

                    # logic to determine if this scan should be plotted
                    skipThisScan = (onlyDoScanNo is not None and scanNoThis != onlyDoScanNo) or (startAtScanNo is not None and scanNoThis < startAtScanNo)
                    skipNextScan = (onlyDoScanNo is not None and scanNoThis+1 != onlyDoScanNo) or (startAtScanNo is not None and scanNoThis+1 < startAtScanNo)
                    doOnlyThisScan = onlyDoScanNo is not None and scanNoThis == onlyDoScanNo

                    # CONTROL logics to determine where we are in relation to scans
                    noScanYet = currentScDir is None and lastScDir is None and currentScDirElev is None and lastScDirElev is None
                    staticPeriod = currentScDir == 'stat' and lastScDir in ['stat', None] and currentScDirElev == 'stat' and lastScDirElev in ['stat', None]
                    startOfVeryFirstScan = (currentScDir in ['cw', 'ccw'] and lastScDir is None) or (currentScDirElev in ['elev incr', 'elev decr'] and lastScDirElev is None)
                    endOfStaticPeriod = lastScDir == 'stat' and lastScDirElev == 'stat' and (currentScDir in ['cw', 'ccw'] or currentScDirElev in ['elev incr', 'elev decr'])
                    sameScanAsBefore = (lastScDir == currentScDir and lastScDir in ['cw', 'ccw']) or (lastScDir == 'stat' and currentScDir == 'stat' and lastScDirElev == currentScDirElev and lastScDirElev in ['elev incr', 'elev decr'])
                    endOfScan_newScan = (currentScDir in ['cw', 'ccw'] or currentScDirElev in ['elev incr', 'elev decr']) and (lastScDir in ['cw', 'ccw'] or lastScDirElev in ['elev incr', 'elev decr']) and (currentScDir != lastScDir or (currentScDirElev != lastScDirElev and currentScDir == 'stat'))
                    endOfScan_static = currentScDir == 'stat' and currentScDirElev == 'stat' and ((lastScDir in ['cw', 'ccw']) ^ (lastScDirElev in ['elev incr', 'elev decr']))

#                    # old control logics commented out below, will only detect azimuthal scans
#                    noScanYet = currentScDir is None and lastScDir is None
#                    staticPeriod = currentScDir == 'stat' and lastScDir in ['stat', None]
#                    startOfVeryFirstScan = currentScDir in ['cw', 'ccw'] and lastScDir is None
#                    endOfStaticPeriod = currentScDir in ['cw', 'ccw'] and lastScDir == 'stat'
#                    sameScanAsBefore = lastScDir == currentScDir and lastScDir in ['cw', 'ccw']
#                    endOfScan_newScan = lastScDir != currentScDir and lastScDir in ['cw', 'ccw'] and currentScDir in ['cw', 'ccw']
#                    endOfScan_static = lastScDir != currentScDir and lastScDir in ['cw', 'ccw'] and currentScDir == 'stat'

                    if (endOfScan_newScan or startOfVeryFirstScan or endOfStaticPeriod) and Time[-1, -1].time() < startTime:
                        skipNextScan = True

                    # check that only one of the above is True
                    controlList = [noScanYet, staticPeriod, startOfVeryFirstScan, endOfStaticPeriod, sameScanAsBefore, endOfScan_newScan, endOfScan_static]
                    if sum(controlList) != 1:
                        raise ScanDetectionError('none or >1 of control logics are True: {}'.format(controlList))

                    # remove altitude lines if elevation-only scan
                    if startOfVeryFirstScan or endOfStaticPeriod or endOfScan_newScan:
                        scDir = currentScDir if currentScDir != 'stat' else currentScDirElev
                        if scDir in ['elev incr', 'elev decr']:
                            useAlts = [400]
                        else:
                            useAlts = alts

                    #################################################
                    # DETERMINE WHAT GETS GONE WITH THE LOADED DATA #
                    #################################################

                    if noScanYet:
                        logging.info('{}: Skipping data file (very first integration period)'.format(fn))
                        pass

                    if staticPeriod:
                        logging.info('{}: Skipping data file (radar stationary)'.format(fn))
                        pass

                    elif startOfVeryFirstScan or endOfStaticPeriod:
                        logging.info('{}: New scan detected (very first scan or end of stationary period)'.format(fn))
                        scanNoThis += 1
                        if not skipNextScan:
                            if RT or debugRT:
                                print('\nInitializing scan #{} (scan start: {})'.format(scanNoThis, Time[0, 0]))
                            par2D_prev, par2D, err2D_prev, err2D = fix_dimensions(par2D_prev, par2D, err2D_prev, err2D)
                            thisScan = Scan(np.append(Time_prev, Time, axis=1), np.append(par1D_prev, par1D, axis=0), np.append(par2D_prev, par2D, axis=1), np.append(err2D_prev, err2D, axis=1), scDir, scanNoThis)
                            thisScan.plotAlts = useAlts

                    elif sameScanAsBefore:
                        logging.info('{}: Adding data to same scan'.format(fn))
                        if not skipThisScan:
                            thisScan.add_data(Time, par1D, par2D, err2D)
                            if debugRT and doPlot:
                                thisScan.plot()

                    elif endOfScan_newScan or endOfScan_static:
                        if endOfScan_newScan:
                            logging.info('{}: End of scan (new scan detected)'.format(fn))
                        else:
                            logging.info('{}: End of scan (entering static period)'.format(fn))
                        if not skipThisScan:
                            logging.info('   Plotting scan')
                            thisScan.finished = True
                            if (RT or debugRT) and RT_replotAfterScan:
                                print('Plotting scan' if thisScan.fig is None else 'Closing realtime scan and replotting')
                                thisScan.closeFig()
                            if doPlot:
                                thisScan.plot()
                                thisScan.saveFig(webAccessFolder, webAccessFolderExternal)
                                if not doOnlyThisScan:
                                    thisScan.closeFig()
                                else:  # return scan in addition to plotting if only plotting this scan
                                    return thisScan
                            else:  # no plotting, append scan to list instead
                                if not doOnlyThisScan:  # append to list if parsing many scans
                                    allScans.append(thisScan)
                                else:  # return only this scan if doing one scan
                                    return thisScan
                        if endOfScan_newScan:
                            scanNoThis += 1
                            if not skipNextScan:
                                if RT or debugRT:
                                    print('\nInitializing scan #{} (scan start: {})'.format(scanNoThis, Time[0, 0]))
                                par2D_prev, par2D, err2D_prev, err2D = fix_dimensions(par2D_prev, par2D, err2D_prev, err2D)
                                thisScan = Scan(np.append(Time_prev, Time, axis=1), np.append(par1D_prev, par1D, axis=0), np.append(par2D_prev, par2D, axis=1), np.append(err2D_prev, err2D, axis=1), scDir, scanNoThis)
                                thisScan.plotAlts = useAlts

                    # save current data for next loop
                    lastAzim = currentAzim
                    lastScDir = currentScDir
                    lastElev = currentElev
                    lastScDirElev = currentScDirElev
                    Time_prev = Time
                    par1D_prev = par1D
                    par2D_prev = par2D
                    err2D_prev = err2D

                ###########################################################
                # ALL THE PREVIOUSLY DISCOVERED DATA FILES ARE NOW LOADED #
                ###########################################################

                # update list of loaded files
                oldfiles = allfiles

                # if realtime, plot what we currently have
                if (RT or doOnlyThisScan) and not skipThisScan:
                    logging.info('End of list with new files. Plotting current data...')
                    thisScan.plot()

                # finished processing previously loaded files, look for more (if RT) or quit
                if RT or debugRT:
                    print('Loading done! Looking for more files...')
                else:
                    print('All files loaded, finalizing last scan...')
                    if not skipThisScan:
                        thisScan.finished = True
                        if doPlot:
                            thisScan.plot()
                            thisScan.saveFig(webAccessFolder, webAccessFolderExternal)
                            thisScan.closeFig()
                        else:
                            if doOnlyThisScan:
                                return thisScan
                            else:
                                allScans.append(thisScan)
                                return allScans
                    break
            else:  # no new files found
                print('No new files, waiting 3 sec... (Ctrl-C to quit and finalize current scan)')
                sleep(3)
    except KeyboardInterrupt:  # user has pressed Ctrl-C, finalize current scan and quit
        print('Aborting, finalizing current scan...'),
        if (RT or debugRT) and RT_replotAfterScan:
            print('Plotting scan' if thisScan.fig is None else 'Closing realtime scan and replotting')
            thisScan.finished = True
            thisScan.closeFig()
        if doPlot:
            thisScan.plot()
            thisScan.saveFig(webAccessFolder, webAccessFolderExternal)
            thisScan.closeFig()
        else:
            if doOnlyThisScan:
                return thisScan
            else:
                allScans.append(thisScan)
        print('Done!')
        return allScans

def make_html(webAccessFolder, webAccessFolderExternal, imageFolder, lastFile):
    '''htmlInternal and htmlExternal point to the same folder, e.g.
          internal: /www_kstdev/display/
          external: http://158.39.70.130/~kstdev/display/
    '''

    webAccessFolder = os.path.abspath(os.path.expanduser(webAccessFolder))
    imageFolder = os.path.abspath(os.path.expanduser(imageFolder))

    # local and remote folders where plots are found
    plotFoldersIn = os.path.join(webAccessFolder, 'eiscatscanplot')
    plotFoldersIn_external = os.path.join(webAccessFolderExternal, 'eiscatscanplot')
    plotsIn = os.path.join(plotFoldersIn, os.path.basename(imageFolder))

    # create plot directory
    if not os.path.isdir(plotsIn) and os.path.isdir(imageFolder):
        shutil.copytree(imageFolder, plotsIn)
    elif not os.path.isdir(plotsIn):
        os.makedirs(plotsIn)

    # copy last file
    shutil.copyfile(os.path.join(imageFolder, lastFile), os.path.join(plotsIn, lastFile))

    # read html template code from source file
    with open('rt_src.html', 'r') as f:
        s1 = f.read()
    with open('rt_src.html', 'r') as f:
        s2 = f.read()

    # make list of files
    files = [fn for fn in os.listdir(imageFolder) if '.png' in fn]
    files.sort()

    # copy latest file to permanent name
    shutil.copyfile(os.path.join(imageFolder, files[-1]), os.path.join(webAccessFolder, 'ESRlatest.png'))

    files_external = [os.path.join(plotFoldersIn_external, os.path.basename(imageFolder), fn) for fn in files]

    # insert file list, latest image and link to image folder
    s1 = s1.replace('[!filenames]', 'filenames = ' + str(files_external))
    s1 = s1.replace('[!latestImg]', files_external[-1])
    s1 = s1.replace('[!latestImgPermalink]', os.path.join(webAccessFolderExternal, 'ESRlatest.png'))
    s1 = s1.replace('[!picFolder]', plotFoldersIn_external)

    # do the same for the plot folder html file
    s2 = s2.replace('[!filenames]', 'filenames = ' + str(files))
    s2 = s2.replace('[!latestImg]', files[-1])
    s2 = s2.replace('[!latestImgPermalink]', os.path.join(webAccessFolderExternal, 'ESRlatest.png'))
    s2 = s2.replace('[!picFolder]', plotFoldersIn_external)

    # write html page
    with open(os.path.join(webAccessFolder, 'scans.html'), 'w') as f:
        f.write(s1)

    # write html page
    with open(os.path.join(plotsIn, '!_SCAN_BROWSER.html'), 'w') as f:
        f.write(s2)

if __name__ == "__main__":

    # datafolder input
    if dataFolder is None:
        dataFolderPrompt = '\nData folder could not be auto-detected, please specify: >> '
    else:
        dataFolderPrompt = '\nData folder [default: {}] >> '.format(dataFolder)
    while True:
        dataFolderOverride = raw_input(dataFolderPrompt)
        if not dataFolderOverride and dataFolder is not None:
            break  # use default
        elif not dataFolderOverride and dataFolder is None:
            print('Please specify a data folder.')
        elif os.path.isdir(dataFolderOverride):
            dataFolder = dataFolderOverride
            break
        else:
            print('Folder {} doesn\'t exist, please try again.'.format(dataFolderOverride))

    savePath = '~/users/eiscatscanplot/plotted/esr_scans_{}'.format(os.path.basename(os.path.normpath(dataFolder)))  # save plotted figures to this path.

    webAccessFolder = '/www_kstdev/display/'
    webAccessFolderExternal = 'http://158.39.70.130/~kstdev/display/'
#    webAccessFolderExternal = 'file:///C:/www_kstdev/display/'
    webAccessFolderMsg = '{} --> {}'.format(os.path.join(webAccessFolder, 'scans.html'), os.path.join(webAccessFolderExternal, 'scans.html')) if webAccessFolder else 'disabled'

    # guess which scan to start at
    startAt = '1'
    startAtMsg = ''
    # guess from files in savePath
    if os.path.isdir(os.path.abspath(os.path.expanduser(savePath))):
        files = [fn for fn in os.listdir(os.path.abspath(os.path.expanduser(savePath))) if '.png' in fn]
        if not len(files) == 0:
            try:
                scanNos = [fn[-7:-4] for fn in files]
                startAt = str(int(max(scanNos)))  # will overwrite the last plot, which might be incomplete
                startAtMsg = ' (auto-detected)'
            except:
                startAtMsg = ' (unable to parse files in plot folder)'

    scanSpeedDegPerSec = 0.63  # scan speed per second
    IPsec = 6.4  # integration period in seconds
    removeLargeErrs = False   # remove data where error > |value|
    alts = [250, 500] if '_tau0_' in dataFolder else [200, 400]  # altitude lines to plot [km]. Set to empty list [] to disable
    drawMag = True

    # additional settings
    while True:
        additionalSettings = raw_input('\nPress Enter to plot or type a number to change the settings below. The following additional settings may be edited.\n\n' +
                                       '   1. Save figures to: {}\n'.format(savePath) +
                                       '   2. Web access: {}\n'.format(webAccessFolderMsg) +
                                       '   3. Start at (scan no. x or time HH:MM): {}\n'.format(startAt) +
                                       '   4. Altitude lines: {}\n'.format(' '.join(map(str, alts))) +
                                       '   5. Draw magnetic parallel/meridian: {}\n'.format(drawMag) +
                                       '   6. Remove data where error > |value|: {}\n'.format(removeLargeErrs) +
                                       '   7. Realtime plotting: {}\n'.format(RT) +
                                       '\nPlease select a number or press Enter to start plotting >> ')
        if not additionalSettings:
            break
        elif additionalSettings == '1':
            # plot save path
            savePathOverride = raw_input('\nSave plots in folder [default: {}] >> '.format(savePath))
            savePath = savePathOverride or savePath
            savePath = os.path.abspath(os.path.expanduser(savePath))
        elif additionalSettings == '2':
            # 'remote web access to scans
            webAccessFolderOverride = raw_input('\nLocal folder in which to put the html file [default: {}, single space to disable] >> '.format(webAccessFolder))
            if webAccessFolderOverride == ' ':
                webAccessFolder = ''
            else:
                latestImagePath = webAccessFolderOverride
                webAccessFolderExternalOverride = raw_input('\nRemote (web) address to the same folder [default: {}] >> '.format(webAccessFolderExternal))
                webAccessFolderExternal = webAccessFolderExternalOverride or webAccessFolderExternal
            webAccessFolderMsg = '{} --> {}'.format(os.path.join(webAccessFolder, 'scans.html'), os.path.join(webAccessFolderExternal, 'scans.html')) if webAccessFolder else 'disabled'
        elif additionalSettings == '3':
            # start scan no. or start time
            startAtOverride = raw_input('\nEnter scan number or HH:MM from which to start [default: scan no. {}{}] >> '.format(startAt, startAtMsg))
            if startAtOverride == '0':
                startAt = '1'
            elif startAtOverride:
                startAt = startAtOverride
        elif additionalSettings == '4':
            # altitude lines
            altsOverride = raw_input('\nDraw lines at which altitudes? Space-separated list of numbers, blank for default, single space to disable [default: {}] >> '.format(' '.join(map(str, alts))))
            if altsOverride:
                alts = map(int, altsOverride.split())
        elif additionalSettings == '5':
            # switch draw magnetic meridian/parallel through ESR
            drawMag = not drawMag
            print('Draw magnetic parallel/meridian through ESR turned {}'.format({True: 'ON', False: 'OFF'}[drawMag]))
        elif additionalSettings == '6':
            # switch error filter
            removeLargeErrs = not removeLargeErrs
            print('Error filter turned {}'.format({True: 'ON', False: 'OFF'}[removeLargeErrs]))
        elif additionalSettings == '7':
            # switch realtime plotting
            RT = not RT
            print('Realtime turned {}'.format({True: 'ON', False: 'OFF'}[RT]))
        else:
            print('Invalid choice')

    print('Starting scan parsing')

    scan_parse(dataFolder=dataFolder, savePath=savePath, doPlot=True, RT=RT,
               onlyDoScanNo=onlyDoScanNo, startAt=startAt, removeLargeErrs=removeLargeErrs, RT_replotAfterScan=RT_replotAfterScan,
               alts=alts, drawMag=drawMag, radarLoc=radarLoc, mapWidth=mapWidth, figSize=figSize,
               debugRT=debugRT, webAccessFolder=webAccessFolder, webAccessFolderExternal=webAccessFolderExternal)

    raw_input('\nPlotting finished, figures saved to {}. Press Enter to close >> '.format(savePath))

max_stars_repo_name: cmeeren/eiscatscanplot

max_issues_repo_name: cmeeren/eiscatscanplot

max_forks_repo_name: cmeeren/eiscatscanplot


text:
# -*- coding: utf-8 -*-
"""
Created on Wed Jun 22 14:06:51 2016

@author: damodara

In this file,
#%%
gradForm1   -- Computes the gradient of the loss function based on formulation 1
     Input:
         Data -- data to be used in the loss function
         W   -- precomputed pseudo random fourier expansion (upto L-1 expansion)
         w   -- pseudo random fourier expansion to be computed (L expansion)
         sigma -- sigma parameter to compute the RBF kernel
         K1  -- precomputed RBF kernel with sigma
    Output:
         L -- is the computed gradient vector
#
gradForm2  -- Computes the gradient of the loss function based on formulation 2
    The description is same as the gradForm1 (see above)
#%%
LossfunctionForm1  -- Computes the loss function based on formulation 1
LossfunctionForm2  -- Computes the loss function based on formulation 2
    Input:
        W  -- Pseudo Random Fourier Coefficient (or Random Fourier Coefficient) matrix
    Output:
        L  -- the computed value of the loss function
ApproxKernelLossfunction1 -- computes the approximated kernel based on formulation 1
ApproxKernelLossfunction2 -- computes the approximated kernel based on formulation 1
"""

#%% This function computes gradient for the loss function based formulation 1
    # gradient of L= 1/2N (k(x1,x2)-cos(W(x1-x2)))
def gradForm1(Data,W,w,sigma, K1=None):
    import numpy as np
    from sklearn.metrics import pairwise

    Nexp=np.shape(W)[1]
    Ndata,Nfeat=np.shape(Data)
    if K1 is None:
        K1=pairwise.rbf_kernel(Data,gamma=sigma**2/2)

    from kernelCentralization import kernelCentralization
    #KC=kernelCentralization(K1)

    nn = np.shape(W)[1]
    C=np.zeros((Nfeat,Ndata))
    B=np.zeros((Nfeat,Ndata))
    for i in range(Ndata):
        X=Data-np.tile(Data[i,],(Ndata,1))
        K=K1[i,:]
        K=K[np.newaxis,]
 #       K=KC[i,:]
#        K=K[np.newaxis,]
        tmp=np.sin(np.dot(X,w))
        B[:,i]=np.dot((X.T), np.diag(np.dot(tmp,K)))
        tmp2=np.dot(tmp,np.ones((1,Nexp)))
        tmp3=(np.cos(np.dot(X,W)))*tmp2
        C[:,i]=np.squeeze(np.dot((X.T),np.sum(tmp3,axis=1)))

    L=(1/Ndata**2)*((np.sum(B,axis=1)) - (1/Nexp)*np.sum(C,axis=1))

    L = L[:,np.newaxis]
    return L

#%% This function computes the value of the Loss function of formulation 1
def LossfunctionForm1(Data,W,gamma,K=None):
    import numpy as np
    from sklearn.metrics import pairwise
    Nexp=np.shape(W)[1]
    Ndata,Nfeat=np.shape(Data)
    Np, Nexp=np.shape(W)
#    if (Np != Nfeat):
#        W = np.delete(W, Np-1,0)
#        Data = np.concatenate((Data, np.ones((Ndata,1), dtype= np.float64)), axis =1)
    W = np.delete(W, Np-1,0)
    # if kernel is not provided
    if K is None:
        K=pairwise.rbf_kernel(Data,gamma=gamma)
    # explicit feature computation
    Dummy = np.concatenate((np.cos(np.dot(Data,W)), np.sin(np.dot(Data,W))), axis=1)
    AK = (1.0/Nexp)*np.dot(Dummy, Dummy.T)
    L = np.linalg.norm(K-AK)/np.linalg.norm(K)
    L = np.sum((K-AK)**2)/(Ndata**2)
    # normal implementation as loss function
#    AK=np.zeros((Ndata,Ndata))
#    for i in range(Ndata):
#        X=Data-np.tile(Data[i,],(Ndata,1))
#        B[i]=np.sum(np.cos(np.dot(X,W)))
#        AK[i,:] = np.sum(np.cos(np.dot(X,W)), axis=1)
#
#    #L=(0.5/Ndata**2)*(A-(1/Nexp)*np.sum(B))**2
#    L = np.linalg.norm(K-(1/Nexp)*AK)/np.linalg.norm(K)
    return L

#%%  This function computes for the gradient of the formulation 2
def gradForm2(Data,W,w,sigma, K1=None):
    # gradient of L= 1/2N (k(x1,x2)-cos(Wx1)cos(Wx2))
    import numpy as np
    from sklearn.metrics import pairwise

    Nexp = np.shape(W)[1]
    Ndata,Nfeat = np.shape(Data)

    if K1 is None:
        K1 = pairwise.rbf_kernel(Data,gamma=sigma**2/2)

    C=np.zeros((Nfeat,Ndata))
    B=np.zeros((Nfeat,Ndata))
    for i in range(Ndata):
        X = np.tile(Data[i,], (Ndata,1))
        K = K1[i,:].copy()
        K = K[:,np.newaxis]
        #
        c1 = np.outer(np.cos(np.dot(Data, w))*K, np.sin(np.dot(Data[i,], w)))
        c2 = np.outer(np.sin(np.dot(Data, w))*K, np.cos(np.dot(Data[i,], w)))
        C[:,i] = np.squeeze((np.dot((X.T), c1)+np.dot((Data.T), c2)))
        #

        AK = np.dot(np.cos(np.dot(Data[i,], W)), np.cos(np.dot(Data,W)).T)
        AK = AK[:,np.newaxis]
        #
        b1 = np.outer(np.cos(np.dot(Data, w))*AK, np.sin(np.dot(Data[i,], w)))
        b2 = np.outer(np.sin(np.dot(Data, w))*AK, np.cos(np.dot(Data[i,], w)))
        B[:,i] = np.squeeze((np.dot((X.T), b1)+np.dot((Data.T), b2)))
        #

    L = (1.0/Ndata**2)*(np.sum(C,axis=1) - (2.0/Nexp)*np.sum(B, axis=1))
    L = L[:,np.newaxis]
    return L

#%%  This function computes loss function value of formulation 2
def LossfunctionForm2(NData,W,gamma,K=None):
    import numpy as np
    from sklearn.metrics import pairwise

    Data = NData.copy()
    Np, Nexp=np.shape(W)
    Ndata,Nfeat=np.shape(Data)

    if (Np != Nfeat):
        Data = np.concatenate((Data, np.ones((Ndata,1), dtype= np.float64)), axis =1)
    # if kernel is not provided
    if K is None:
        K=pairwise.rbf_kernel(Data,gamma=gamma)
    Phi = np.cos(np.dot(Data,W))
    AK= (2.0/Nexp)*np.dot(Phi, Phi.T)
    #L = np.linalg.norm(K-AK)/np.linalg.norm(K)
    L = np.sum((K-AK)**2)/(Ndata**2)

    ## normal implementation as the loss function
    #A=np.sum(K)
    #B=np.zeros((Ndata,1))
#    for i in range(Ndata):
#        X=Data-np.tile(Data[i,],(Ndata,1))
#        B[i]=np.sum(np.cos(np.dot(X,W)))

    #L=(0.5/Ndata**2)*(A-(1/Nexp)*np.sum(B))**2

    return L

#%%  ######################################################################
def ApproxKernelLossfunction1(Data,W,sigma,K=None):
    import numpy as np
    from sklearn.metrics import pairwise
    Nexp=np.shape(W)[1]
    Ndata,Nfeat=np.shape(Data)
    # if kernel is not provided
    if K is None:
        K=pairwise.rbf_kernel(Data,gamma=sigma**2/2.0)
    #

    Dummy = np.concatenate((np.cos(np.dot(Data,W)), np.sin(np.dot(Data,W))), axis=1)
    AK = (2.0/Nexp)*np.dot(Dummy, Dummy.T)
    L = np.linalg.norm(K-AK)/np.linalg.norm(K)
    # normal implementation as loss function
#    AK=np.zeros((Ndata,Ndata))
#    for i in range(Ndata):
#        X=Data-np.tile(Data[i,],(Ndata,1))
#        #B[i]=np.sum(np.cos(np.dot(X,W)))
#        AK[i,:] = np.sum(np.cos(np.dot(X,W)), axis=1)
#
#    #L=(0.5/Ndata**2)*(A-(1/Nexp)*np.sum(B))**2
#    AK = (1/Nexp)*AK
#    L = np.linalg.norm(K-AK)/np.linalg.norm(K)
    return L, AK
#%%
def ApproxKernelLossfunction2(NData,W,sigma,K=None):
    import numpy as np
    from sklearn.metrics import pairwise
    Data = NData.copy()
    Nexp=np.shape(W)[1]
    Ndata,Nfeat=np.shape(Data)
    # if kernel is not provided
    if K is None:
        K=pairwise.rbf_kernel(Data,gamma=sigma**2/2.0)
    Phi =  np.cos(np.dot(Data,W))
    AK= (2.0/Nexp)*np.dot(Phi, Phi.T)
    L = np.linalg.norm(K-AK)/np.linalg.norm(K)
    return L, AK


text:
import naturalize.crossover.core as c
import naturalize.crossover.strategies as st
from naturalize.solutionClass import Individual
import numpy as np
import pytest

# from naturalize.crossover.strategies import crossGene, basicCrossover
from naturalize.solutionClass import Individual
import numpy as np

np.random.seed(25)

# df.defaultFitness(individual, env)
genea = np.empty(5)
geneb = np.empty(5)

individualA = Individual([genea])
individualB = Individual([geneb])
# genea = np.array([1,1,1,1,1,0])
# geneb = np.array([1,2,2,2,2,1])

basicCrossover = c.getCrossover()

def test_crossGene():

    out1, out2 = st.crossGeneSingleCut(genea, geneb)

    check1 = np.all(genea[:4] == out1[:4]) and (geneb[4] == out2[4])
    check2 = np.all(geneb[:4] == out2[:4]) and (genea[4] == out1[4])

    assert np.all(check1 == check2)

def test_default_crossGene():

    Indout1, Indout2 = basicCrossover(individualA, individualB)

    genea = individualA.genotype[0]
    geneb = individualB.genotype[0]

    out1 = Indout1.genotype[0]
    out2 = Indout2.genotype[0]

    check1 = np.all(genea[:2] == out1[:2]) and (geneb[2:] == out2[2:])
    check2 = np.all(geneb[:2] == out2[:2]) and (genea[2:] == out1[2:])

    assert np.all(check1 == check2)

# test_crossGene()
# test_default_crossGene()

max_stars_repo_name: cslotboom/Naturalize

max_issues_repo_name: cslotboom/Naturalize

max_forks_repo_name: cslotboom/Naturalize