text:
import numpy as np
from netCDF4 import Dataset
import qc

class Climatology:
    """
    Class for dealing with climatologies, reading, extracting values etc.
    Automatically detects if this is a single field, pentad or daily climatology
    """

    def __init__(self, infield):
        """
        Read in the climatology for variable var from infile

        :param infield: numpy array containing the climatology
        :type infield: numpy array
        """
        self.field = infield
        self.n = self.field.shape[0]
        assert self.n in [1, 73, 365], 'weird shaped field'
        self.res = 180. / self.field.shape[1]

    @classmethod
    def from_filename(cls, infile, var):
        """
        Read in the climatology for variable var from infile

        :param infile: filename of a netcdf file
        :param var: the variable name to be extracted from the netcdf file
        :type infile: string
        :type var: string
        """
        if infile is not None:
            climatology = Dataset(infile)
            field = climatology.variables[var][:]

            lat_synonyms = ['lat', 'lats', 'latitude', 'latitudes']
            found_lat = False
            for lat in lat_synonyms:
                if lat in climatology.variables:
                    latitudes = climatology.variables[lat][:]
                    found_lat = True
            assert found_lat, 'no readable latitude information in NetCDF file: ' + infile

            lon_synonyms = ['lon', 'lons', 'long', 'longs', 'longitude', 'longitudes']
            found_lon = False
            for lon in lon_synonyms:
                if lon in climatology.variables:
                    longitudes = climatology.variables[lon][:]
                    found_lon = True
            assert found_lon, 'no readable longitude information in NetCDF file: ' + infile

            climatology.close()

            # transpose the fields if the second axis is longitude
            if (field.shape[0] == 1 and
                    field.shape[1] == 360 and
                    field.shape[2] == 180):
                field = field.transpose(0, 2, 1)

            # added an exception here as the OSTIA background variances are specified from S to N so first
            # element in latitudes is negative
            if latitudes[0] < 0:
                field = np.flip(field, 1)

            # if the longitudes start near zero then roll the array along its longitude axis
            if longitudes[0] > 0.0 and longitudes[0] < 1.0:
                lon_len = field.shape[2]
                field = np.roll(field, lon_len / 2, axis=2)

            if field.ndim == 4:
                field = field[:, 0, :, :]

        else:
            field = np.ma.array(np.zeros((1, 180 * 20, 360 * 20)), mask=True)

        return cls(field)

    def get_tindex(self, month, day):
        """
        Get the time index of the input month and day

        :param month: month for which the time index is required
        :param day: day for which the time index is required
        :type month: integer
        :type day: integer
        :return: time index for specified month and day.
        :rtype: integer
        """

        tindex = None

        if self.n == 1:
            tindex = 0
        if self.n == 73:
            tindex = qc.which_pentad(month, day) - 1
        if self.n == 365:
            tindex = qc.day_in_year(month, day) - 1

        return tindex

    def get_value_ostia(self, lat, lon):
        """
        :param lat: latitude of location to extract value from in degrees of arc
        :param lon: longitude of location to extract value from in degrees of arc
        :return: SST at that location or None

        :type lat: float
        :type lon: float
        :rtype: float
        """
        yindex = qc.mds_lat_to_yindex(lat, res=0.05)
        xindex = qc.mds_lon_to_xindex(lon, res=0.05)
        tindex = 0

        result = self.field[tindex, yindex, xindex]

        if type(result) is np.float64 or type(result) is np.float32:
            pass
        else:
            if result.mask:
                result = None
            else:
                result = result.data[0]

        return result

    def get_value_mds_style(self, lat, lon, month, day):
        """
        Get the value from the climatology at the give position and time using the MDS
        method for deciding which grid cell borderline cases fall into

        :param lat: latitude of location to extract value from in degrees
        :param lon: longitude of location to extract value from in degrees
        :param month: month for which the value is required
        :param day: day for which the value is required
        :type lat: float
        :type lon: float
        :type month: integer
        :type day: integer
        :return: climatology value at specified location and time.
        :rtype: float
        """
        if month < 1 or month > 12:
            return None
        ml = qc.get_month_lengths(2004)
        if day < 1 or day > ml[month - 1]:
            return None

        yindex = qc.mds_lat_to_yindex(lat)
        xindex = qc.mds_lon_to_xindex(lon)
        tindex = self.get_tindex(month, day)

        result = self.field[tindex, yindex, xindex]

        if type(result) is np.float64 or type(result) is np.float32:
            pass
        else:
            if result.mask:
                result = None
            else:
                result = result.data[0]

        return result

    def get_value(self, lat, lon, month, day):
        """
        Get the value from the climatology at the give position and time

        :param lat: latitude of location to extract value from in degrees
        :param lon: longitude of location to extract value from in degrees
        :param month: month for which the value is required
        :param day: day for which the value is required
        :type lat: float
        :type lon: float
        :type month: integer
        :type day: integer
        :return: climatology value at specified location and time.
        :rtype: float
        """
        if month < 1 or month > 12:
            return None
        ml = qc.get_month_lengths(2004)
        if day < 1 or day > ml[month - 1]:
            return None

        yindex = qc.lat_to_yindex(lat, self.res)
        xindex = qc.lon_to_xindex(lon, self.res)
        tindex = self.get_tindex(month, day)

        result = self.field[tindex, yindex, xindex]

        if type(result) is np.float64 or type(result) is np.float32:
            pass
        else:
            if result.mask:
                result = None
            else:
                result = result.data[0]

        return result

    def get_interpolated_value(self, lat, lon, mo, dy):
        """
        Get the value from the climatology interpolated to the precise location of
        the observation in time and space

        :param lat: latitude of location to extract value from in degrees
        :param lon: longitude of location to extract value from in degrees
        :param mo: month for which the value is required
        :param dy: day for which the value is required
        :type lat: float
        :type lon: float
        :type mo: integer
        :type dy: integer
        :return: climatology value at specified location and time.
        :rtype: float
        """
        # check that the lat lon point falls in a grid cell with a value or on
        # the border of one
        if lat + 0.001 < 90:
            pert1 = self.get_value(lat + 0.001, lon + 0.001, mo, dy)
            pert2 = self.get_value(lat + 0.001, lon - 0.001, mo, dy)
        else:
            pert1 = None
            pert2 = None
        if lat - 0.001 > -90:
            pert3 = self.get_value(lat - 0.001, lon + 0.001, mo, dy)
            pert4 = self.get_value(lat - 0.001, lon - 0.001, mo, dy)
        else:
            pert3 = None
            pert4 = None

        if (pert1 is None and pert2 is None and
                pert3 is None and pert4 is None):
            return None

        x1, x2, y1, y2 = qc.get_four_surrounding_points(lat, lon, 1)

        try:
            q11 = self.get_value(y1, x1, mo, dy)
        except:
            q11 = None
        if q11 is not None:
            q11 = float(q11)

        try:
            q22 = self.get_value(y2, x2, mo, dy)
        except:
            q22 = None
        if q22 is not None:
            q22 = float(q22)

        try:
            q12 = self.get_value(y2, x1, mo, dy)
        except:
            q12 = None
        if q12 is not None:
            q12 = float(q12)

        try:
            q21 = self.get_value(y1, x2, mo, dy)
        except:
            q21 = None
        if q21 is not None:
            q21 = float(q21)

        q11, q12, q21, q22 = qc.fill_missing_vals(q11, q12, q21, q22)

        x1, x2, y1, y2 = qc.get_four_surrounding_points(lat, lon, 0)

        return qc.bilinear_interp(x1, x2, y1, y2,
                                  lon, lat,
                                  q11, q12, q21, q22)

max_stars_repo_name: ET-NCMP/MarineQC

max_issues_repo_name: ET-NCMP/MarineQC

max_forks_repo_name: ET-NCMP/MarineQC


text:
"""
Objective: Dashboard of the NSE Nifty & NSE Next Index stocks and its distribution.

Plots: Weightages, Returns (1yr, 6m, 3m, 1m) and Free Float Market Capitalization.
"""

import datetime
from datetime import date, timedelta
import matplotlib.pyplot as plt
import mplfinance as mpf
from mplfinance.original_flavor import candlestick_ohlc
import matplotlib.dates as mdates
import nsepy
from nsepy import get_history
import numpy as np
import pandas as pd
from sklearn import preprocessing

#Top Level Settings:
plt.style.use('seaborn-darkgrid')

#================================================================================
#Top 10 nifty stock weightage
#================================================================================

t = datetime.datetime.today()

day = t.strftime('%d')
month = t.strftime('%m')
year = t.strftime('%y')

#Dynamic URL
top10_url = 'https://www1.nseindia.com/content/indices/top10nifty50_{}{}{}.csv'.format(day, month, year)

df_weightage = pd.DataFrame()

try:
    df_weightage = pd.read_csv(top10_url)
    if not df_weightage.empty:
        print('Updating cache file.')
        df_weightage.to_csv('nifty_weightage.csv')
except:
    print ("Error reading file online: Reading cached file.")
    try:
        df_weightage = pd.read_csv('nifty_weightage.csv', index_col = 0)
    except:
        print ("Error: File Not Found: Did not find 'nifty_weightage.csv' file.")

else:
    if df_weightage.empty == True:
        print ("Error: File Not Found: Did not find 'nifty_weightage.csv' file.")
        print ("Error: Data Read: Check if the date is valid business day.")

#================================================================================
#Monthly Index Returns
#================================================================================

#Static URL
monthly_index_returns = 'https://www1.nseindia.com/content/indices/mir.csv'

df_monthly_index_returns = pd.read_csv(monthly_index_returns, skiprows = 2)
df_monthly_index_returns.rename(columns={'Unnamed: 0': 'indices'}, inplace=True)

#================================================================================
#Market Capitalization and Weightages for Nifty and Next nifty
#================================================================================

#Static URL
next_url = 'https://www1.nseindia.com/content/indices/niftynext50_mcwb.csv'
nifty_url = 'https://www1.nseindia.com/content/indices/nifty50_mcwb.csv'

df_nifty_mcap = pd.read_csv(nifty_url, skiprows = 2, index_col = 0)
df_next_mcap = pd.read_csv(next_url, skiprows = 2, index_col = 0)

df_columns_names_weightage = {
                                'Security Symbol' : 'symbol',
                                'Security Name' : 'security',
                                'Industry' : 'industry',
                                'Equity Capital (In Rs.)' : 'equity',
                                'Free Float Market Capitalisation (Rs. Crores)' : 'FFMC',
                                'Weightage (%)' : 'weightage',
                                'Volatility (%)' : 'volatility',
                                'Monthly Return' : 'monthly_return',
                                'Avg. Impact Cost (%)' : 'impact_cost'
}

df_nifty_mcap.rename(columns = df_columns_names_weightage, inplace = True)
df_nifty_mcap = df_nifty_mcap[0:49]
df_nifty_mcap['industry'] = df_nifty_mcap['industry'].str[:15]
df_nifty_mcap.sort_values('FFMC', ascending = False, inplace = True)

df_next_mcap.rename(columns = df_columns_names_weightage, inplace = True)
df_next_mcap = df_next_mcap[0:49]
df_next_mcap['industry'] = df_next_mcap['industry'].str[:15]
df_next_mcap.sort_values('FFMC', ascending = False, inplace = True)

#================================================================================
# Plotting Charts
#================================================================================

#Plot Settings
fig = plt.figure(figsize=(15,20))

#Index Returns
ax1 = plt.subplot2grid((5, 4), (0, 0), rowspan = 1, colspan = 3)

#Top 10 weightage
ax2 = plt.subplot2grid((5, 4), (0, 3), rowspan = 1, colspan = 1)

#Monthly returns
ax6 = plt.subplot2grid((5, 4), (1, 3), rowspan = 2, colspan = 1)
ax3 = plt.subplot2grid((5, 4), (1, 0), rowspan = 2, colspan = 1, sharey=ax6)
ax4 = plt.subplot2grid((5, 4), (1, 1), rowspan = 2, colspan = 1, sharey=ax6)
ax5 = plt.subplot2grid((5, 4), (1, 2), rowspan = 2, colspan = 1, sharey=ax6)

#FreeFloat & MarketCaps & Weightages
ax7 = plt.subplot2grid((5, 4), (3, 0), rowspan = 1, colspan = 2)
ax8 = plt.subplot2grid((5, 4), (4, 0), rowspan = 1, colspan = 2)
ax9 = plt.subplot2grid((5, 4), (3, 2), rowspan = 1, colspan = 2)
ax10 = plt.subplot2grid((5, 4), (4, 2), rowspan = 1, colspan = 2)

#================================================================================
# Charts & Settings
#================================================================================

#Plotting Monthly Index Returns
colored = ['tab:green' if x else 'tab:red' for x in (df_monthly_index_returns['1 year'] >= 0)]
ax1.barh(df_monthly_index_returns['indices'], df_monthly_index_returns['1 year'], color = colored)
plt.xticks(rotation=90)

#Basic bar plot of Nifty Top 10 weights
k_temp = df_weightage['WEIGHTAGE(%)'].sum().round(2)
ax2.barh(df_weightage['SYMBOL'], df_weightage['WEIGHTAGE(%)'], color = 'tab:blue')
#ax2 = df_weightage.plot.bar('SYMBOL', 'WEIGHTAGE(%)')
ax2.annotate('TotalWeight={}'.format(k_temp), (5, 9))
#ax2.xticks(rotation=90)
#ax2.set_xticklabels(df_weightage['SYMBOL'], rotation=90, ha='center')

#Plotting Monthly Index Returns - 1 year, 6 months, 3 months, 1 month
#
colored_ax3 = ['tab:green' if x else 'tab:red' for x in (df_monthly_index_returns['1 year'] >= 0)]
ax3.barh(df_monthly_index_returns['indices'], df_monthly_index_returns['1 year'], color = colored_ax3)
ax3.axes.get_yaxis().set_visible(True)
ax3.set_title('1-year', loc='center', color = 'tab:gray')
ax3.grid(which='minor', axis='both', color = 'tab:gray', linestyle='--', linewidth=1 )

colored_ax4 = ['tab:green' if x else 'tab:red' for x in (df_monthly_index_returns['6 month'] >= 0)]
ax4.barh(df_monthly_index_returns['indices'], df_monthly_index_returns['6 month'], color = colored_ax4)
ax4.axes.get_yaxis().set_visible(False)
ax4.set_title('6-month', loc='center', color = 'tab:gray')
ax4.grid(which='minor', axis='both', color = 'tab:gray', linestyle='--', linewidth=1 )

colored_ax5 = ['tab:green' if x else 'tab:red' for x in (df_monthly_index_returns['3 month'] >= 0)]
ax5.barh(df_monthly_index_returns['indices'], df_monthly_index_returns['3 month'], color = colored_ax5)
ax5.axes.get_yaxis().set_visible(False)
ax5.set_title('3-month', loc='center', color = 'tab:gray')
ax5.grid(which='minor', axis='both', color = 'tab:gray', linestyle='--', linewidth=1 )

colored_ax6 = ['tab:green' if x else 'tab:red' for x in (df_monthly_index_returns['1 month'] >= 0)]
ax6.barh(df_monthly_index_returns['indices'], df_monthly_index_returns['1 month'], color = colored_ax6)
ax6.axes.get_yaxis().set_visible(False)
ax6.set_title('1-month', loc='center', color = 'tab:gray')
ax6.grid(which='minor', axis='both', color = 'tab:gray', linestyle='--', linewidth=1 )

#================================================================================
#First Column - Plotting
#ax7 = df_nifty_mcap.plot(x='symbol', y= 'FFMC', kind = 'bar', ax=ax7, label = 'Free Float Market Cap (Rs Crs)', secondary_y=True)
ax7.bar(df_nifty_mcap['symbol'], df_nifty_mcap['FFMC'], label = 'Free Float Market Cap (Rs Crs)')
#ax7.set_xlabel("Security")
#ax7.set_ylabel("Market Capitalization in (Rs Crs)")
ax7.set_xticklabels(df_nifty_mcap['symbol'].str.upper().str[:8], rotation=90, ha='center')
ax7.yaxis.tick_right()
ax7.legend(['Free Float Market Cap (Rs Crs)'])

xmin, xmax = ax7.get_xlim()
mean_FFMC_nifty = round(df_nifty_mcap["FFMC"].mean(), 2)
mean_FFMC_displace_nifty = mean_FFMC_nifty * 1.05

ax7.axhline(y=df_nifty_mcap["FFMC"].mean(), xmin=0, xmax=1, color = 'black', alpha = 0.5)
ax7.xaxis.grid(False, which='major', color = 'whitesmoke', linestyle='--', linewidth=0.02)
ax7.yaxis.grid(True, which='major', color = 'gray', linestyle='--', linewidth=0.25)
ax7.annotate(mean_FFMC_nifty, (xmax-10, mean_FFMC_displace_nifty) )

#FirstColumn = Industry Weightage chart
c = df_nifty_mcap.pivot_table(index = ['security', 'industry'], values = ['FFMC','weightage']).sort_values('weightage', ascending = False)
Weight = c.groupby('industry')[['weightage', 'FFMC']].sum().sort_values('FFMC', ascending = False)
ax8 = Weight.plot(kind = 'bar', secondary_y= 'weightage', grid = True, ax=ax8)
ax8.xaxis.grid(True, which='major', color = 'gray', linestyle='--', linewidth=0.25)
ax8.yaxis.grid(True, which='major', color = 'gray', linestyle='--', linewidth=0.25)
ax8.get_xaxis().get_label().set_visible(False)

#================================================================================

#Second Column - Plotting
#ax9 = df_next_mcap.plot(x='symbol', y= 'FFMC', kind = 'bar', ax=ax9, label = 'Free Float Market Cap (Rs Crs)', secondary_y=True)
ax9.bar(df_next_mcap['symbol'], df_next_mcap['FFMC'], label = 'Free Float Market Cap (Rs Crs)')
#ax9.set_xlabel("Security")
#ax9.set_ylabel("Market Capitalization in (Rs Crs)")
ax9.set_xticklabels(df_next_mcap['symbol'].str.upper().str[:8], rotation=90, ha='center')
ax9.yaxis.tick_right()
ax9.legend(['Free Float Market Cap (Rs Crs)'])

xmin, xmax = ax9.get_xlim()
mean_FFMC_next = round(df_next_mcap["FFMC"].mean(), 2)
mean_FFMC_displace_next = mean_FFMC_next * 1.05

ax9.axhline(y=df_next_mcap["FFMC"].mean(), xmin=0, xmax=1, color = 'black', alpha = 0.5)
ax9.xaxis.grid(False, which='major', color = 'whitesmoke', linestyle='--', linewidth=0.02)
ax9.yaxis.grid(True, which='major', color = 'gray', linestyle='--', linewidth=0.25)
ax9.annotate(mean_FFMC_next, (xmax-10, mean_FFMC_displace_next) )

#SecondColumn = Industry Weightage chart
c = df_next_mcap.pivot_table(index = ['security', 'industry'], values = ['FFMC','weightage']).sort_values('weightage', ascending = False)
Weight = c.groupby('industry')[['weightage', 'FFMC']].sum().sort_values('FFMC', ascending = False)
ax10 = Weight.plot(kind = 'bar', secondary_y= 'weightage', grid = True, ax=ax10)
ax10.xaxis.grid(True, which='major', color = 'gray', linestyle='--', linewidth=0.25)
ax10.yaxis.grid(True, which='major', color = 'gray', linestyle='--', linewidth=0.25)
ax10.get_xaxis().get_label().set_visible(False)

#================================================================================

plt.subplots_adjust(wspace=0.1, hspace=0.15)
fig.tight_layout()

ext: py

reason: import numpy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Author: Benjamin Vial
# License: MIT

"""
Band diagram of 2D photonic crystal
===================================

Calculation of the band diagram of a two-dimensional photonic crystal.
"""

# sphinx_gallery_thumbnail_number = -1

import matplotlib.pyplot as plt
import numpy as np

import gyptis as gy
from gyptis import pi

##############################################################################
# Reference results are taken from  :cite:p:`Joannopoulos2008` (Chapter 5 Fig. 2).
#
# The structure is a square lattice of dielectric
# columns, with radius r and dielectric constant :math:`\varepsilon`.
# The material is invariant along the z direction  and periodic along
# :math:`x` and :math:`y` with lattice constant :math:`a`.
# We will define the geometry using the class :class:`~gyptis.Lattice`,
# defining the two vectors of periodicity:

a = 1  # unit cell size
vectors = (a, 0), (0, a)  # vectors defining the unit cell
R = 0.2 * a  # inclusion radius

lattice = gy.Lattice(dim=2, vectors=vectors)

##############################################################################
# Next, we add a cylinder and compute the boolean fragments

circ = lattice.add_circle(a / 2, a / 2, 0, R)
circ, cell = lattice.fragment(circ, lattice.cell)

##############################################################################
# One needs to define physical domains associated with the basic geometrical
# entities:

lattice.add_physical(cell, "background")
_ = lattice.add_physical(circ, "inclusion")

##############################################################################
# Set minimum mesh size

lattice.set_size("background", a / 10)
lattice.set_size("inclusion", R / 10)

##############################################################################
# Finally, we can build the geometry, which will also construct the mesh.

lattice.build()

##############################################################################
# Material parameters are defined with a python dictionary:

epsilon = dict(background=1, inclusion=8.9)
mu = dict(background=1, inclusion=1)

##############################################################################
# We can now instanciate the simulation class :class:`~gyptis.PhotonicCrystal`.
# We will compute eigenpairs at the :math:`X` point of the Brillouin zone, *i.e.*
# the propagation vector is :math:`\mathbf k = (\pi/a,0)`.

phc = gy.PhotonicCrystal(
    lattice,
    epsilon,
    mu,
    propagation_vector=(pi / a, 0),
    polarization="TE",
    degree=2,
)

##############################################################################
# To calculate the eigenvalues and eigenvectors, we call the
# :meth:`~gyptis.PhotonicCrystal.eigensolve` method.

solution = phc.eigensolve(n_eig=6, wavevector_target=0.5)

##############################################################################
# The results can be accessed through the `phc.solution` attribute
# (a dictionary).

ev_norma = phc.solution["eigenvalues"].real * a / (2 * pi)
print("Normalized eigenfrequencies")
print("---------------------------")
print(ev_norma)

##############################################################################
# Lets plot the field map of the modes.

eig_vects = phc.solution["eigenvectors"]
for mode, eval in zip(eig_vects, ev_norma):
    gy.plot(mode.real, cmap="gyptis")
    plt.title(fr"$\omega a/2\pi c = {eval:0.3f}$")
    H = phc.formulation.get_dual(mode, 1)
    gy.dolfin.plot(H.real, cmap="Greys")
    lattice.plot_subdomains()
    plt.axis("off")

##############################################################################
# We define here the wavevector path:

Nb = 21
K = np.linspace(0, pi / a, Nb)
bands = np.zeros((3 * Nb - 3, 2))
bands[:Nb, 0] = K
bands[Nb : 2 * Nb, 0] = K[-1]
bands[Nb : 2 * Nb - 1, 1] = K[1:]
bands[2 * Nb - 1 : 3 * Nb - 3, 0] = bands[2 * Nb - 1 : 3 * Nb - 3, 1] = np.flipud(K)[
    1:-1
]

##############################################################################
# Calculate the band diagram:

BD = {}
for polarization in ["TE", "TM"]:
    ev_band = []
    for kx, ky in bands:
        phc = gy.PhotonicCrystal(
            lattice,
            epsilon,
            mu,
            propagation_vector=(kx, ky),
            polarization=polarization,
            degree=1,
        )
        phc.eigensolve(n_eig=6, wavevector_target=0.1)
        ev_norma = phc.solution["eigenvalues"].real * a / (2 * pi)
        ev_band.append(ev_norma)
    # append first value since this is the same point
    ev_band.append(ev_band[0])
    BD[polarization] = ev_band

##############################################################################
# Plot the bands:

bands_plot = np.zeros(3 * Nb - 2)
bands_plot[:Nb] = K
bands_plot[Nb : 2 * Nb - 1] = K[-1] + K[1:]
bands_plot[2 * Nb - 1 : 3 * Nb - 2] = 2 * K[-1] + 2 ** 0.5 * K[1:]

plt.figure(figsize=(3.2, 2.5))
plotTM = plt.plot(bands_plot, BD["TM"], c="#4199b0")
plotTE = plt.plot(bands_plot, BD["TE"], c="#cf5268")
plt.annotate("TM modes", (1, 0.05), c="#4199b0")
plt.annotate("TE modes", (0.33, 0.33), c="#cf5268")
plt.ylim(0, 0.8)
plt.xlim(0, bands_plot[-1])
plt.xticks(
    [0, K[-1], 2 * K[-1], bands_plot[-1]], ["$\Gamma$", "$X$", "$M$", "$\Gamma$"]
)
plt.axvline(K[-1], c="k", lw=0.3)
plt.axvline(2 * K[-1], c="k", lw=0.3)
plt.ylabel(r"Frequency $\omega a/2\pi c$")
plt.tight_layout()

max_stars_repo_name: benvial/gyptis

max_issues_repo_name: benvial/gyptis

max_forks_repo_name: benvial/gyptis


text:
import numpy as np

# select dataset from: krogan2006_core, krogan2006_extended, collins2007
dataset = 'krogan2006_core'

P = np.load('p_' + dataset + '.npy')
print('finished loading matrix P for dataset ' + dataset)

n = len(P)

# recommended k value for
# krogran2006_core (validation against SGD): 1350
# krogran2006_core (validation against MIPS): 950
# krogran2006_extended (validation against SGD, MIPS): 1900
# collins2007 (validation against SGD): 800
# collins2007 (validation against MIPS): 450
k = 950

J = []
for j in range(k):
    col_norms = np.sum(np.multiply(P, P), axis = 0)
    max_ind = np.where(col_norms == np.max(col_norms))[0][0]
    J.append(max_ind)
    v = np.reshape(P[:, max_ind], (n, 1))
    mat1 = np.dot(np.transpose(v), P)
    mat2 = np.dot(v, mat1)/(col_norms[max_ind])
    P -= mat2
    if j % 10 == 0:
        print('found ' + str(j) + ' almost pure nodes out of ' + str(k))

# save pure nodes
with open('J_' + dataset + '.txt', 'w') as fp:
    for item in J:
        fp.write("%s\n" % item)
print('finished computing and saving pure nodes for dataset ' + dataset)


text:
# Copyright (c) 2017 Chao Ji

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
# ==============================================================================
import numpy as np

from autodiff.core import base_node

class _BinaryOp(base_node.Node):
  """Base class for all binary arithmetic operations. Supports broadcasting.
  """
  def __init__(self, x, y, graph=None):
    """Constructor.

    Args:
      x: a Node instance, the first operand.
      y: a Node instance, the second operand.
      graph: a Graph instance.
    """
    shape = self._compute_shape(x._shape, y._shape)
    super(_BinaryOp, self).__init__(shape, graph)
    self._arguments['x'] = x
    self._arguments['y'] = y

  def _compute_shape(self, x_shape, y_shape):
    """Computs the shape of the resulting Node.

    Args:
      x_shape: a TensorShape instance, shape of the first operand.
      y_shape: a TensorShape instance, shape of the second operand.

    Returns:
      a tuple of integers (or None), shape of the resulting Node.
    """
    if x_shape._ndims == 0:
      return y_shape._raw_shape
    elif y_shape._ndims == 0:
      return x_shape._raw_shape

    shape = []
    min_ndims = np.minimum(x_shape._ndims, y_shape._ndims)
    for i in range(min_ndims):
      dx = x_shape._raw_shape[-1 - i]
      dy = y_shape._raw_shape[-1 - i]
      if dx is not None and dy is not None:
        if dx == dy:
          shape.append(dx)
        elif dx == 1:
          shape.append(dy)
        elif dy == 1:
          shape.append(dx)
        else:
          raise ValueError('operands x(%s) and y(%s) have incompatible shapes for '
              'broadcasting.' % (x_shape, y_shape))
      elif dx is not None:
        shape.append(None if dx == 1 else dx)
      elif dy is not None:
        shape.append(None if dy == 1 else dy)
      else:
        shape.append(None)

    if min_ndims == x_shape._ndims:
      return tuple(y_shape._raw_shape[:-min_ndims]) + tuple(shape[::-1])
    else: # min_ndims == y_shape._ndims
      return tuple(x_shape._raw_shape[:-min_ndims]) + tuple(shape[::-1])

  def _compute_reduce_dims(
      self, dynamic_x_shape, dynamic_y_shape, dynamic_shape):
    """Compute list of dimensions to be reduces.

    Args:
      dynamic_x_shape: a tuple of ints, dynamic shape of the first operand.
      dynamic_y_shape: a tuple of ints, dynamic shape of the second operand.
      dynamic_shape: a tuple of ints, dynamic shape of the resulting Node.
    """
    if len(dynamic_shape) > len(dynamic_x_shape):
      x_shape_expand = tuple(np.ones(len(dynamic_shape) - len(dynamic_x_shape),
          dtype=int)) + dynamic_x_shape
    else:
      x_shape_expand = dynamic_x_shape
    if len(dynamic_shape) > len(dynamic_y_shape):
      y_shape_expand = tuple(np.ones(len(dynamic_shape) - len(dynamic_y_shape),
          dtype=int)) + dynamic_y_shape
    else:
      y_shape_expand = dynamic_y_shape

    reduce_dims_x = tuple(np.where(np.array(dynamic_shape) !=
        np.array(x_shape_expand))[0])
    reduce_dims_y = tuple(np.where(np.array(dynamic_shape) !=
        np.array(y_shape_expand))[0])
    return reduce_dims_x, reduce_dims_y

class Add(_BinaryOp):
  """Add operation. Supports broadcasting."""
  def __init__(self, x, y, graph=None):
    """Constructor.

    Args:
      x: a Node instance, the first operand.
      y: a Node instance, the second operand.
      graph: a RunTime instance.
    """
    super(Add, self).__init__(x, y, graph)

  def _forward(self, feed_dict):
    """Compute the forward pass value of the node.

    Args:
      feed_dict: a dict mapping from a `Node` instance to a numpy array.

    Returns:
      the forward pass value of the node.
    """
    return (self._arguments['x'].forward(feed_dict) +
        self._arguments['y'].forward(feed_dict))

  def _backward(self, feed_dict):
    """Retrieve the gradient value of the current node. Then compute and
    backprop the gradient w.r.t. the argument nodes of the current node.

    Args:
      feed_dict: a dict mapping from a `Node` instance to a numpy array.

    Returns:
      grad_dict: a dict mapping from a `Node` instance to a numpy array, holding
        the gradient w.r.t. `self`'s arguments that pass through `self`.
    """
    grad_val = self._graph.get_runtime()._bwval[self.name]

    dynamic_x_shape = self._arguments['x'].forward(feed_dict).shape
    dynamic_y_shape = self._arguments['y'].forward(feed_dict).shape
    dynamic_shape = self.forward(feed_dict).shape

    reduce_dims_x, reduce_dims_y = self._compute_reduce_dims(
        dynamic_x_shape, dynamic_y_shape, dynamic_shape)

    dx_val = grad_val.sum(axis=reduce_dims_x).reshape(dynamic_x_shape)
    dy_val = grad_val.sum(axis=reduce_dims_y).reshape(dynamic_y_shape)

    if self._arguments['x'] is self._arguments['y']:
      grad_dict = {self._arguments['x']: np.ones(dynamic_x_shape,
          dtype='float32') * 2}
    else:
      grad_dict = {self._arguments['x']: dx_val,
                   self._arguments['y']: dy_val}
    return grad_dict

class Multiply(_BinaryOp):
  """Multiply operation. Supports broadcasting."""
  def __init__(self, x, y, graph=None):
    """Constructor.

    Args:
      x: a Node instance, the first operand.
      y: a Node instance, the second operand.
      graph: a RunTime instance.
    """
    super(Multiply, self).__init__(x, y, graph)

  def _forward(self, feed_dict):
    """Compute the forward pass value of the node.

    Args:
      feed_dict: a dict mapping from a `Node` instance to a numpy array.

    Returns:
      the forward pass value of the node.
    """
    return (self._arguments['x'].forward(feed_dict) *
        self._arguments['y'].forward(feed_dict))

  def _backward(self, feed_dict):
    """Retrieve the gradient value of the current node. Then compute and
    backprop the gradient w.r.t. the argument nodes of the current node.

    Args:
      feed_dict: a dict mapping from a `Node` instance to a numpy array.

    Returns:
      grad_dict: a dict mapping from a `Node` instance to a numpy array, holding
        the gradient w.r.t. `self`'s arguments that pass through `self`.
    """
    grad_val = self._graph.get_runtime()._bwval[self.name]

    x_val = self._arguments['x'].forward(feed_dict)
    y_val = self._arguments['y'].forward(feed_dict)
    dynamic_x_shape = x_val.shape
    dynamic_y_shape = y_val.shape
    dynamic_shape = self.forward(feed_dict).shape

    reduce_dims_x, reduce_dims_y = self._compute_reduce_dims(
        dynamic_x_shape, dynamic_y_shape, dynamic_shape)

    dx_val = (grad_val * y_val).sum(axis=reduce_dims_x).reshape(dynamic_x_shape)
    dy_val = (grad_val * x_val).sum(axis=reduce_dims_y).reshape(dynamic_y_shape)

    if self._arguments['x'] is self._arguments['y']:
      grad_dict = {self._arguments['x']: dx_val + dy_val}
    else:
      grad_dict = {self._arguments['x']: dx_val,
                   self._arguments['y']: dy_val}

    return grad_dict

class Subtract(_BinaryOp):
  """Subtract operation. Supports broadcasting."""
  def __init__(self, x, y, graph=None):
    """Constructor.

    Args:
      x: a Node instance, the first operand.
      y: a Node instance, the second operand.
      graph: a RunTime instance.
    """
    super(Subtract, self).__init__(x, y, graph)

  def _forward(self, feed_dict):
    """Compute the forward pass value of the node.

    Args:
      feed_dict: a dict mapping from a `Node` instance to a numpy array.

    Returns:
      the forward pass value of the node.
    """
    return (self._arguments['x'].forward(feed_dict) -
        self._arguments['y'].forward(feed_dict))

  def _backward(self, feed_dict):
    """Retrieve the gradient value of the current node. Then compute and
    backprop the gradient w.r.t. the argument nodes of the current node.

    Args:
      feed_dict: a dict mapping from a `Node` instance to a numpy array.

    Returns:
      grad_dict: a dict mapping from a `Node` instance to a numpy array, holding
        the gradient w.r.t. `self`'s arguments that pass through `self`.
    """
    grad_val = self._graph.get_runtime()._bwval[self.name]

    dynamic_x_shape = self._arguments['x'].forward(feed_dict).shape
    dynamic_y_shape = self._arguments['y'].forward(feed_dict).shape
    dynamic_shape = self.forward(feed_dict).shape

    reduce_dims_x, reduce_dims_y = self._compute_reduce_dims(
        dynamic_x_shape, dynamic_y_shape, dynamic_shape)

    dx_val = grad_val.sum(axis=reduce_dims_x).reshape(dynamic_x_shape)
    dy_val = grad_val.sum(axis=reduce_dims_y).reshape(dynamic_y_shape)

    if self._arguments['x'] is self._arguments['y']:
      grad_dict = {self._arguments['x']: np.zeros(dynamic_x_shape,
          dtype='float32')}
    else:
      grad_dict = {self._arguments['x']: dx_val,
                   self._arguments['y']: -dy_val}
    return grad_dict

class Exp(base_node.Node):
  """Computes the elementwise Exp(x) function."""
  def __init__(self, x, graph=None):
    """Constructor.

    Args:
      x: a Node instance, the input tensor.
      graph: a RunTime instance.
    """
    super(Exp, self).__init__(x._shape, graph)
    self._arguments['x'] = x

  def _forward(self, feed_dict):
    """Compute the forward pass value of the node.

    Args:
      feed_dict: a dict mapping from a `Node` instance to a numpy array.

    Returns:
      the forward pass value of the node.
    """
    x_val = self._arguments['x'].forward(feed_dict)
    return np.exp(x_val)

  def _backward(self, feed_dict):
    """Retrieve the gradient value of the current node. Then compute and
    backprop the gradient w.r.t. the argument nodes of the current node.

    Args:
      feed_dict: a dict mapping from a `Node` instance to a numpy array.

    Returns:
      grad_dict: a dict mapping from a `Node` instance to a numpy array, holding
        the gradient w.r.t. `self`'s arguments that pass through `self`.
    """
    val = self.forward(feed_dict)
    grad_val = self._graph.get_runtime()._bwval[self.name]
    grad_val = grad_val * val
    grad_dict = {self._arguments['x']: grad_val}
    return grad_dict

max_stars_repo_name: chao-ji/np-auto-diff


text:
from __future__ import print_function
from future.utils import string_types
import os
import re
import glob, copy
from collections import OrderedDict as odict
from functools import partial
import warnings
import numpy as np
import netCDF4 as nc
import dimarray as da
#from geo.data.index import to_slice, _slice3D

from dimarray.tools import format_doc, isscalar
from dimarray.dataset import Dataset, concatenate_ds, stack_ds
from dimarray.core import DimArray, Axis, Axes
from dimarray.config import get_option
from dimarray.core.bases import AbstractDimArray, AbstractDataset, AbstractAxis, GetSetDelAttrMixin, AbstractAxes
# from dimarray.core.metadata import _repr_metadata
from dimarray.prettyprinting import repr_axis, repr_axes, repr_dimarray, repr_dataset, repr_attrs

from .conventions import encode_cf_datetime, decode_cf_datetime

__all__ = ['read_nc','summary_nc', 'write_nc']

#
# Global variables
#
FORMAT = get_option('io.nc.format') # for the doc

#
# Helper functions
#
def maybe_encode_values(values, format=None):
    """ strings are given "object" type in Axis object
    ==> assume all objects are actually strings
    NOTE: this will fail for other object-typed axes such as tuples
    """
    # if dtype is np.dtype('O'):
    values = np.asarray(values)
    dtype = values.dtype
    cf_attrs = {}

    if dtype.kind in ('S','O'):
        encoded = np.asarray(values, dtype=object)
        dtype = str
    elif dtype.kind == 'M': # 'datetime64'
        # values = np.asarray(np.datetime_as_string(values), dtype=object)
        # dtype = str
        encoded, units, calendar = encode_cf_datetime(values)
        dtype = encoded.dtype
        cf_attrs['calendar'] = calendar
        cf_attrs['units'] = units
    elif dtype.kind == 'm': # 'timedelta64'
    # elif np.issubdtype(dtype, np.timedelta64):
        encoded, units = encode_cf_timedelta(values)
        cf_attrs['units'] = units
        dtype = encoded.dtype
    else:
        encoded = values

    # if the format is NETCDF3, uses int32 instead of int64
    if format == "NETCDF3_CLASSIC" and dtype is np.dtype("int64"):
        warnings.warn("convert int64 into int32 for writing to NETCDF3 format", RuntimeWarning)
        dtype = "int32"
    elif format == "NETCDF3_64BIT" and dtype is np.dtype("int64"):
        # it is strange that a warnings also occurs in the _64BIT version !
        warnings.warn("convert int64 into int32 for writing to NETCDF3 format", RuntimeWarning)
        dtype = "int32"
    return encoded, dtype, cf_attrs

# quick and dirty conversion from string
TIME_UNITS = ['years','months','days','hours','minutes', 'seconds']

def hastimeunits(ncvar):
    regexpr = "{} +since".format("|".join(TIME_UNITS))
    return hasattr(ncvar, 'units') and re.match(regexpr, ncvar.units)

def istimevariable(ncvar):
    """Return True if a netCDF4 variable represents time
    """
    # return hastimeunits(ncvar) or \
    return (
        ncvar.size > 0 and isinstance(ncvar[0], string_types) \
         and re.match('\d\d(\d\d)?-\d\d-\d\d',ncvar[0]))

class NCTimeVariableWrapper(object):
    """Subclass of netCDF4 variable that performs conversions
    to datetime64 object when indexing. Should only be used for
    object or string types
    """
    def __init__(self, ncvar):
        self.ncvar = ncvar
    def __getitem__(self, idx):
        arr = self.ncvar[idx]
        if arr.dtype.kind in ('S','O'):
            try:
                arr = np.asarray(arr, dtype='datetime64')
            except Exception as error:
                warnings.warn("Failed to convert string time to datetime64\n"+str(error), RuntimeWarning)
        else:
            try:
                arr = decode_cf_datetime(arr, getattr(self.ncvar, 'units',None), getattr(self.ncvar, 'calendar', None))
            except Exception as error:
                warnings.warn("Failed to convert {} to datetime64\n".format(arr.dtype)+str(error), RuntimeWarning)
        return arr
    def __getattr__(self, att):
        return getattr(self.ncvar, att)

#
# Define an open_nc function return an on-disk Dataset object, more similar to
# netCDF4's Dataset.
#

#
# Specific to NetCDF I/O
#
class NetCDFOnDisk(object):
    " DimArray wrapper for netCDF4 Variables and Datasets "
    @property
    def nc(self):
        return self._ds
    @property
    def _obj(self):
        return NotImplementedError("need to be subclassed")
    @property
    def attrs(self):
        return AttrsOnDisk(self._obj)
    def close(self):
        return self._ds.close()
    def __exit__(self, exception_type, exception_value, tracebook):
        self.close()
    def __enter__(self):
        return self

class DatasetOnDisk(GetSetDelAttrMixin, NetCDFOnDisk, AbstractDataset):
    """ Dataset on Disk

    .. versionadded :: 0.2

    See open_nc for examples of use.
    """
    @format_doc(format=FORMAT)
    def __init__(self, f, mode="r", clobber=True, diskless=False, persist=False, format=FORMAT, **kwargs):
        """
        Parameters
        ----------
        f : file name or netCDF4.Dataset instance
        mode : str, optional
            access mode. This include "r" (read) (the default), "w" (write)
            and "a" (append)
            See netCDF4.Dataset for full information.
        clobber : `bool`, optional
            if True (default), opening a file with mode='w'
            will clobber an existing file with the same name.  if False, an
            exception will be raised if a file with the same name already exists.
        diskless : see netCDF4.Dataset help
        persist : see netCDF4.Dataset help
        format : `str`
            netCDF file format. Default is '{format}' (only accounted for during file creation)

        **kwargs : key-word arguments, used as default argument when creating a variable

        Notes
        -----
        See the netCDF4-python module documentation for more information about the use
        of keyword arguments to DatasetOnDisk
        """
        if isinstance(f, nc.Dataset):
            self._ds = f
        else:
            try:
                self._ds = nc.Dataset(f, mode=mode, clobber=clobber, diskless=diskless, persist=persist, format=format)
            except UserWarning as error:
                print(error)
            except Exception as error: # indicate file name when issuing error
                raise IOError("{}\n=> failed to open {} (mode={}, clobber={})".format(str(error), f, mode, clobber)) # easier to handle
        # assert self._kwargs is not None
        self._kwargs = kwargs
    @property
    def _obj(self):
        return self._ds

    def read(self, names=None, indices=None, axis=0, indexing=None, tol=None, keepdims=False,
             verbose=False, # back-compatibility
             ):
        """ Read values from disk

        Parameters
        ----------
        names : list of variables to read, optional
        indices : int or list or slice (single-dimensional indices)
                   or a tuple of those (multi-dimensional)
                   or `dict` of { axis name : axis indices }
            Indices refer to Dataset axes. Any item that does not possess
            one of the dimensions will not be indexed along that dimension.
            For example, scalar items will be left unchanged whatever indices
            are provided.
        axis : None or int or str, optional
            if specified and indices is a slice, scalar or an array, assumes
            indexing is along this axis.
        indexing : {'label', 'position'}, optional
            Indexing mode.
            - "label": indexing on axis labels (default)
            - "position": use numpy-like position index
            Default value can be changed in dimarray.rcParams['indexing.by']
        tol : float, optional
            tolerance when looking for numerical values, e.g. to use nearest
            neighbor search, default `None`.
        keepdims : bool, optional
            keep singleton dimensions (default False)

        Returns
        -------
        Dataset

        See Also
        --------
        open_nc : examples of use
        DimArrayOnDisk.read, DatasetOnDisk.write, DimArray.take
        """
        if verbose:
            # print "Read ",self.filename
            pass
        # automatically read all variables to load (except for the dimensions)
        if names is None:
            dims = self.dims
            names = self.keys()
        elif isinstance(names, string_types):
            return self[names].read(indices=indices, axis=axis, indexing=indexing, tol=tol, keepdims=keepdims)
        else:
            dims = []
        # else:
        #     raise TypeError("Expected list or str for 'names=', got {}".format(names))

        tuple_indices = self._get_indices(indices, axis=axis, tol=tol, keepdims=keepdims, indexing=indexing)
        dict_indices = {dim:tuple_indices[i] for i, dim in enumerate(self.dims)}

        data = Dataset()

        # first load dimensions
        for dim in dims:
            data.axes.append(self.axes[dim][dict_indices[dim]])

        # then normal variables
        for nm in names:
            data[nm] = self[nm].read(indices={dim:dict_indices[dim] for dim in self[nm].dims}, indexing='position')
        data.attrs.update(self.attrs) # dataset's metadata

        # reorder the axes in the dataset to match input
        data.axes = Axes(data.axes[dim] for dim in self.dims if dim in data.dims)

        return data

    # so that loc, iloc, sel, isel, nloc work:
    _getitem = read

    def __getitem__(self, name):
        if name in self.dims:
            return DimArrayOnDisk(self._ds, name)
            # raise KeyError("Use 'axes' property to access dimension variables: {}".format(name))
            # return AxisOnDisk(self._ds, name)

        if name not in self.keys():
            raise KeyError("Variable not found in the dataset: {}.\nExisting variables: {}".format(name, self.keys()))
        return DimArrayOnDisk(self._ds, name)

    # def write(self, name, dataset, zlib=False, **kwargs):
    def write(self, name, dima, **kwargs):
        """ Write a variable to a netCDF4 dataset.

        Parameters
        ----------
        name : variable name
        dima : DimArray instance
        zlib : `bool`
            Enable zlib compression if True. Default is False (no compression).
        complevel : `int`
            integer between 1 and 9 describing the level of compression desired. Ignored if zlib=False.
        **kwargs : key-word arguments
            Any additional keyword arguments accepted by `netCDF4.Dataset.createVariable`

        Notes
        -----
        See the netCDF4-python module documentation for more information about the use
        of keyword arguments to write_nc.

        See also
        --------
        DimArrayOnDisk.write

        Examples
        --------
        >>> import numpy as np
        >>> import dimarray as da

        Create a DimArray (in memory), with metadata

        >>> dima = da.DimArray([[1,2,3],[4,5,6]], axes=[('time',[2000,2045.5]),('scenario',['a','b','c'])])
        >>> dima.units = 'myunits' # metadata
        >>> dima.axes['time'].units = 'metadata-dim-in-memory'

        Write it to disk, and add some additional metadata

        >>> ds = da.open_nc('/tmp/test.nc', mode='w')
        >>> ds['myvar'] = dima
        >>> ds['myvar'].bla = 'bla'
        >>> ds['myvar'].axes['time'].yo = 'metadata-dim-on-disk'
        >>> ds.axes['scenario'].ya = 'metadata-var-on-disk'
        >>> ds.yi = 'metadata-dataset-on-disk'
        >>> ds.close()

        Check the result with ncdump utility from a terminal (need to be installed for the test below to work)

        :> ncdump -h /tmp/test.nc
        netcdf test {
        dimensions:
            time = 2 ;
            scenario = 3 ;
        variables:
            double time(time) ;
                time:units = "metadata-dim-in-memory" ;
                time:yo = "metadata-dim-on-disk" ;
            string scenario(scenario) ;
                scenario:ya = "metadata-var-on-disk" ;
            int64 myvar(time, scenario) ;
                myvar:units = "myunits" ;
                myvar:bla = "bla" ;

        // global attributes:
                :yi = "metadata-dataset-on-disk" ;
        }
        """
        # if not isinstance(obj, da.DimArray):
        #     raise TypeError("Can only write Dataset, use `ds[name] = dima` to write a DimArray")
        _, nctype, cf_attrs = maybe_encode_values(dima, format=self._ds.file_format)

        name = name or getattr(self, "name", None)
        if not name:
            raise ValueError("Need to provide variable name")

        kw = self._kwargs.copy()
        kw.update(kwargs)

        if name not in self._ds.variables.keys():

            if np.isscalar(dima):
                dima = da.DimArray(dima)
            if not isinstance(dima, DimArray):
                raise TypeError("Expected DimArray, got {}".format(type(dima)))
            if dima.ndim > 0:
                # create Dimension and associated variable
                for ax in dima.axes:
                    if ax.name not in self.dims:
                        self.axes.append(ax, **kw)

            # add _FillValue or missing_value attribute to fill_value
            if 'fill_value' not in kw and hasattr(dima, '_FillValue'):
                kw['fill_value'] = dima._FillValue
            elif 'fill_value' not in kw and hasattr(dima, 'missing_value'):
                kw['fill_value'] = dima.missing_value

            # if the variable is already present (e.g. as present as variable
            # AND dimension in dimarray Dataset), just skip it
            if  name in self._ds.variables.keys():
                pass
                # warnings.warn("dimension variable present both as variable and Axis in DimArra", RuntimeWarning)
            else:
                self._ds.createVariable(name, nctype, dima.dims, **kw)

        dimaondisk = DimArrayOnDisk(self._ds, name)
        dimaondisk[()] = dima
        dimaondisk.attrs.update(cf_attrs) # calendar?

    __setitem__ = write

    def __delitem__(self, name):
        " will raise an Exception if on-disk "
        del self._ds.variables[name]

    def keys(self):
        " all variables except for dimensions"
        return [nm for nm in self._ds.variables.keys() if nm not in self.dims]

    def values(self):
        return [self[k] for k in self.keys()]

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

    @property
    def dims(self):
        return tuple(self._ds.dimensions.keys())
    @dims.setter
    def dims(self, newdims):
        self._set_dims(newdims)

    @property
    def axes(self):
        return AxesOnDisk(self._ds, self.dims, **self._kwargs)
    @axes.setter
    def axes(self, newaxes):
        for ax in newaxes:
            if ax.name in self.dims:
                self.axes[ax.name] = ax
            else:
                self.axes.append(ax)

    _repr = repr_dataset

    def __iter__(self):
        for k in self.keys():
            yield k

class NetCDFVariable(NetCDFOnDisk):
    @property
    def _obj(self):
        return self.values
    @property
    def values(self):
        return self._ds.variables[self._name] # simply the variable to be indexed and returns values like netCDF4
    @values.setter
    def values(self, values):
        self[()] = values
    @property
    def __array__(self):
        return self.values[()].__array__ # returns a numpy array

class DimArrayOnDisk(GetSetDelAttrMixin, NetCDFVariable, AbstractDimArray):
    _constructor = DimArray
    _broadcast = False
    def __init__(self, ds, name, _indexing=None):
        self._indexing = _indexing or get_option('indexing.by')
        self._ds = ds
        self._name = name
    @property
    def name(self):
        return self._name
    @name.setter
    def name(self, newname):
        self._ds.renameVariable(self, self._name, newname)
        self._name = newname
    @property
    def dims(self):
        return tuple(self._obj.dimensions)
    @property
    def axes(self):
        return AxesOnDisk(self._ds, self.dims)

    def write(self, indices, values, axis=0, indexing=None, tol=None):
        """ Write numpy array or DimArray to netCDF file (The
        variable must have been created previously)

        Parameters
        ----------
        indices : int or list or slice (single-dimensional indices)
                   or a tuple of those (multi-dimensional)
                   or `dict` of { axis name : axis values }
        values : np.ndarray or DimArray
        axis : None or int or str, optional
            if specified and indices is a slice, scalar or an array, assumes
            indexing is along this axis.
        indexing : {'label', 'position'}, optional
            Indexing mode.
            - "label": indexing on axis labels (default)
            - "position": use numpy-like position index
            Default value can be changed in dimarray.rcParams['indexing.by']
        tol : float, optional
            tolerance when looking for numerical values, e.g. to use nearest
            neighbor search, default `None`.

        See Also
        --------
        DimArray.put, DatasetOnDisk.write, DimArrayOnDisk.read
        """
        # just create the variable and dimensions
        ds = self._ds
        dima = values # internal convention: values is a numpy array
        values, nctype, cf_attrs = maybe_encode_values(values, format=self._ds.file_format)

        assert self._name in ds.variables.keys(), "variable does not exist, should have been created earlier!"

        # add attributes
        if hasattr(dima,'attrs'):
            self.attrs.update(dima.attrs)
            self.attrs.update(cf_attrs) # calendar?

        # special case: index == slice(None) and self.ndim == 0
        # This would fail with numpy, but not with netCDF4
        if type(indices) is slice and indices == slice(None) and self.ndim == 0:
            indices = ()

        indices = self._get_indices(indices,axis=axis, indexing=indexing, tol=tol)

        # Perform additional checks on axes if the Data to assign is a DimArray
        if isinstance(dima, DimArray):
            for i, ax in enumerate(self.axes):
                idx = indices[i]
                axis = dima.axes[ax.name]
                # write unlimited dimensions
                if self._ds.dimensions[ax.name].isunlimited():
                    self.axes[ax.name][idx] = axis
                else:
                    # dimension variable already written, simple check
                    ondisk = self.axes[ax.name][idx if not (np.isscalar(idx) or np.ndim(idx)==0)  else [idx]].values
                    inmemory = axis.values
                    # inmemory, _ = maybe_encode_values(axis.values)
                    if not np.all(ondisk == inmemory):
                        # assert np.all(ondisk == inmemory)
                        warnings.warn("axes values differ in the netCDF file, try using a numpy array instead, or re-index\
                                      the dimarray prior to assigning to netCDF", RuntimeWarning)

        # do all the indexing and assignment via IndexedArray class
        # ==> it will set the values via _setvalues_ortho below
        # super(DimArrayOnDisk)._setitem(self, indices, values, **kwargs)
        AbstractDimArray._setitem(self, indices, values, indexing='position')

    def read(self, indices=None, axis=0, indexing=None, tol=None, keepdims=False):
        """ Read values from disk

        Parameters
        ----------
        indices : int or list or slice (single-dimensional indices)
                   or a tuple of those (multi-dimensional)
                   or `dict` of { axis name : axis values }
        axis : None or int or str, optional
            if specified and indices is a slice, scalar or an array, assumes
            indexing is along this axis.
        indexing : {'label', 'position'}, optional
            Indexing mode.
            - "label": indexing on axis labels (default)
            - "position": use numpy-like position index
            Default value can be changed in dimarray.rcParams['indexing.by']
        tol : float, optional
            tolerance when looking for numerical values, e.g. to use nearest
            neighbor search, default `None`.
        keepdims : bool, optional
            keep singleton dimensions (default False)

        Returns
        -------
        DimArray instance or scalar

        See Also
        --------
        DimArray.take, DatasetOnDisk.read, DimArrayOnDisk.write
        """
        if type(indices) is slice and indices == slice(None) and self.ndim == 0:
            indices = ()
        return AbstractDimArray._getitem(self, indices=indices, axis=axis,
                                         indexing=indexing, tol=tol, keepdims=keepdims)

    __setitem__ = write
    __getitem__ = read

    # so that loc, iloc, sel, isel, nloc work:
    _getitem = read
    _setitem = write

    def _repr(self, **kwargs):
        assert 'lazy' not in kwargs, "lazy parameter cannot be provided, it is always True"
        return repr_dimarray(self, lazy=True, **kwargs)

    def _setvalues_ortho(self, idx_tuple, values, cast=False):
        if cast is True:
            warnings.warn("`cast` parameter is ignored", DeprecationWarning)
        # values = maybe_encode_values(values)
        values, nctype, cf_attrs = maybe_encode_values(values)
        self.values[idx_tuple] = values
        self.attrs.update(cf_attrs)

    def _getvalues_ortho(self, idx_tuple):
        res = self.values[idx_tuple]
        # scalar become arrays with netCDF4# scalar become arrays with netCDF4
        # need convert to ndim=0 numpy array for consistency with axes
        if self.ndim == 0:
            try:
                res[0] + 1 # pb arises only for numerical types
                res = np.array(res[0])
            except:
                res = np.array(res) # str and unicode
        return res

    def _getaxes_ortho(self, idx_tuple):
        " idx: tuple of position indices  of length = ndim (orthogonal indexing)"
        axes = []
        for i, ix in enumerate(idx_tuple):
            ax = self.axes[i][ix]
            if np.ndim(ax) != 0: # do not include scalar axes
                axes.append(ax)
        return axes

class AxisOnDisk(GetSetDelAttrMixin, NetCDFVariable, AbstractAxis):
    def __init__(self, ds, name):
        self._ds = ds
        if not isinstance(name, string_types):
            raise TypeError("only string names allowed")
        self._name = name

    def _isdefined(self):
        return self._name in self._ds.variables
    def _ncvar(self):
        return self._ds.variables[self._name]
    def _ncdim(self):
        return self._ds.dimensions[self._name]
    def _istime(self):
        return self._isdefined() and istimevariable(self._ncvar())

    @property
    def values(self):
        if self._isdefined():
            ncvar = self._ncvar()
            if self._istime():
                values = NCTimeVariableWrapper(ncvar)
            else:
                values = ncvar
        else:
            # default, dummy dimension axis
            msg = "'{}' dimension not found, define integer range".format(self._name)
            warnings.warn(msg, RuntimeWarning)
            values = np.arange(len(self._ncdim()))
        return values

    @property
    def name(self):
        return self._name
    @name.setter
    def name(self, new):
        self._ds.renameDimension(self._name, new)
        if self._name in self._ds.variables.keys():
            self._ds.renameVariable(self._name, new)
        self._name = new

    @property
    def attrs(self):
        if self._istime():
            return AttrsOnDisk(self._obj, exclude=['calendar','units'])
        else:
            return AttrsOnDisk(self._obj)

    def __setitem__(self, indices, ax):
        ds = self._ds
        name = self._name

        assert getattr(ax, 'name', name) == name, "inconsistent axis name"

        values, nctype, cf_attrs = maybe_encode_values(ax, format=self._ds.file_format)

        # assign value to variable (finer-grained control in AxesOnDisk.append)
        if name not in self._ds.variables.keys():
            ds.createVariable(name, nctype, name)

        if np.isscalar(indices):
            assert values.size == 1
            try:
                values = np.asscalar(values)
            except AttributeError:
                values = values[()]

        ds.variables[name][indices] = values

        # add attributes
        attrs = getattr(ax, 'attrs', {})
        try:
            self.attrs.update(attrs)
        except Exception as error:
            warnings.warn(str(error), RuntimeWarning)
        self.attrs.update(cf_attrs)

    def __getitem__(self, indices):
        values = self.values[indices]

        if isinstance(values, np.ma.MaskedArray):
            values = values.filled(0)

        # do not produce an Axis object
        if np.isscalar(values):
            return values
        elif np.ndim(values) == 0:
            return np.asscalar(values)

        ax = Axis(values, self._name)

        # add metadata
        if self._name in self._ds.variables.keys():
            ax.attrs.update(self.attrs)

        return ax

    def __len__(self):
        return len(self._ds.dimensions[self._name])

    @property
    def size(self):
        return len(self)

    def _bounds(self):
        size = len(self) # does not work with unlimited dimensions?
        if self.size == 0:
            first, last = None, None
        elif self._name in self._ds.variables.keys():
            first, last = self._ds.variables[self._name][0], self._ds.variables[self._name][size-1]
        else:
            first, last = 0, size-1
        return first, last

    @property
    def dtype(self):
        if self._name in self._ds.variables.keys():
            return np.dtype(self._ds.variables[self._name].dtype)
        else:
            return np.dtype('i')

class AxesOnDisk(AbstractAxes):
    def __init__(self, ds, dims, **kwargs):
        self._ds = ds
        self._dims = dims
        self._kwargs = kwargs # for variable creation

    @property
    def dims(self):
        return self._dims

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

    def __getitem__(self, dim):
        if not isinstance(dim, string_types):
            dim = self.dims[dim]
        elif dim not in self.dims:
            msg = "{} not found in dimensions (dims={}).".format(dim, self.dims)
            print("""An error ocurred? A new dimension can be created via `axes.append(axis)` syntax,
where axis can be a string (unlimited dimension), an Axis instance,
or a tuple (name, values).  See help on `axes.append` for more information.
Low-level netCDF4 function is also available as `ds.nc.createDimension`
and `ds.nc.createVariable`""")
            raise KeyError(msg)
        return AxisOnDisk(self._ds, dim)

    def __setitem__(self, dim, ax):
        " modify existing axis and possibly create new associated variable "
        if not isinstance(dim, string_types):
            dim = self.dims[dim]

        if dim not in self._ds.dimensions.keys():
            raise ValueError("{} dimension does not exist. Use axes.append() to create a new axis".format(dim))

        values, nctype, cf_attrs = maybe_encode_values(ax, format=self._ds.file_format)

        if dim not in self._ds.variables.keys():
            v = self._ds.createVariable(dim, nctype, dim, **self._kwargs)

        self[dim][:] = ax
        self[dim].attrs.update(cf_attrs)
        # v[:] = values

    def append(self, ax, size=None, **kwargs):
        """ create new dimension

        Parameters
        ----------
        ax : str or Axis or (name, values) tuple
            if str, simply create a dimension without writing the axis values
        size : int or None, optional
            if None, an unlimited dimension is created
            if ax is provided as an array-like or Axis, size is taken
            from the axis values by default (so size does not need to
            be provided)
        **kwargs : passed to createVariable (compression parameters)
        """
        if isinstance(ax, string_types):
            self._ds.createDimension(ax, size)
            self._dims += (ax,)

        else:
            try:
                ax = Axis.as_axis(ax)
            except:
                raise TypeError("can only append Axis instances or (name, values),\
                                or provide str to create unlimited dimension")
            # elif isinstance(ax, Axis):
            size = size or ax.size
            self._ds.createDimension(ax.name, size)
            self._dims += (ax.name,)

            self._kwargs.update(kwargs) # createVariable parameters
            self[ax.name] = ax # assign values and attributes

    def __iter__(self):
        for dim in self.dims:
            yield AxisOnDisk(self._ds, dim)

    _repr = repr_axes

class AttrsOnDisk(object):
    """ represent netCDF Dataset or Variable Attribute
    """
    def __init__(self, obj, exclude=None):
        self.obj = obj
        self.exclude = exclude
    def __setitem__(self, name, value):
        if name == "_FillValue": return
        try:
            self.obj.setncattr(name, value)
        except TypeError as error:
            warnings.warn(str(error), RuntimeWarning)
            if type(value) is bool:
                warnings.warn("convert boolean to integer", RuntimeWarning)
                self.obj.setncattr(name, int(value))

    def __getitem__(self, name):
        return self.obj.getncattr(name)
    def __delitem__(self, name):
        return self.obj.delncattr(name)
    def update(self, attrs):
        for k in attrs.keys():
            self[k] = attrs[k]
    def keys(self):
        keys = getattr(self.obj,'ncattrs', lambda :[])()
        if self.exclude is not None:
            keys = [k for k in keys if k not in self.exclude]
        return keys
    def values(self):
        return [self.obj.getncattr(k) for k in self.obj.ncattrs()]
    def todict(self):
        return odict(zip(self.keys(), self.values()))
    def __iter__(self):
        for k in self.keys():
            yield k
    def __len__(self):
        return len(self.keys())
    def __repr__(self):
        return repr_attrs(self)

###################################################
# Wrappers
###################################################

def open_nc(file_name, *args, **kwargs):
    """ open a netCDF file a la netCDF4, for interactive access to its properties

    Parameters
    ----------
    file_name : netCDF file name
    *args, **kwargs : passed to netCDF4.Dataset

    Returns
    -------
    dimarray.io.nc.DatasetOnDisk class (see help on this class for more info)

    Examples
    --------
    >>> ncfile = da.get_ncfile('greenland_velocity.nc')
    >>> ds = da.open_nc(ncfile)

    Informative Display similar to a in-memory Dataset

    >>> ds    # doctest: +SKIP
    DatasetOnDisk of 6 variables (NETCDF4)
    0 / y1 (113): -3400000.0 to -600000.0
    1 / x1 (61): -800000.0 to 700000.0
    surfvelmag: ('y1', 'x1')
    lat: ('y1', 'x1')
    lon: ('y1', 'x1')
    surfvely: ('y1', 'x1')
    surfvelx: ('y1', 'x1')
    mapping: nan

    Load variables with [:] syntax, like netCDF4 package

    >>> ds['surfvelmag'][:] # load one variable
    dimarray: 6893 non-null elements (0 null)
    0 / y1 (113): -3400000.0 to -600000.0
    1 / x1 (61): -800000.0 to 700000.0
    array(...)

    Indexing is similar to DimArray, and includes the sel, isel methods

    >>> ds['surfvelmag'].ix[:10, -1] # load first 10 y1 values, and last x1 value
    dimarray: 10 non-null elements (0 null)
    0 / y1 (10): -3400000.0 to -3175000.0
    array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], dtype=float32)

    >>> ds['surfvelmag'].sel(x1=700000, y1=-3400000)
    dimarray: 1 non-null elements (0 null)
    array(0., dtype=float32)

    >>> ds['surfvelmag'].isel(x1=-1, y1=0)
    dimarray: 1 non-null elements (0 null)
    array(0., dtype=float32)

    Need to close the Dataset at the end

    >>> ds.close() # close

    Also usable as context manager

    >>> with da.open_nc(ncfile) as ds:
    ...     dataset = ds.read() # load full data set, same as da.read_nc(ncfile)
    >>> dataset
    Dataset of 6 variables
    0 / y1 (113): -3400000.0 to -600000.0
    1 / x1 (61): -800000.0 to 700000.0
    surfvelmag: ('y1', 'x1')
    lat: ('y1', 'x1')
    lon: ('y1', 'x1')
    surfvely: ('y1', 'x1')
    surfvelx: ('y1', 'x1')
    mapping: nan
    """
    return DatasetOnDisk(file_name, *args, **kwargs)

def read_nc(f, names=None, *args, **kwargs):
    """  Wrapper around DatasetOnDisk.read

    Read one or several variables from one or several netCDF file

    Parameters
    ----------
    f : str or netCDF handle
        netCDF file to read from or regular expression

    names : None or list or str, optional
        variable name(s) to read
        default is None

    indices : int or list or slice (single-dimensional indices)
               or a tuple of those (multi-dimensional)
               or `dict` of { axis name : axis indices }
        Indices refer to Dataset axes. Any item that does not possess
        one of the dimensions will not be indexed along that dimension.
        For example, scalar items will be left unchanged whatever indices
        are provided.

    indexing : {'label', 'position'}, optional
        Indexing mode.
        - "label": indexing on axis labels (default)
        - "position": use numpy-like position index
        Default value can be changed in dimarray.rcParams['indexing.by']

    tol : float, optional
        tolerance when looking for numerical values, e.g. to use nearest
        neighbor search, default `None`.

    keepdims : bool, optional
        keep singleton dimensions (default False)

    axis : str, optional
        When reading multiple files, axis along which to join
        the dimarrays or datasets. It the axis already exist,
        the resulting arrays will be concatenated, otherwise
        they will be stacked along a new array (in the sense
        of the numpy functions `concatenate` and `stack`)

    keys : sequence, optional
        When reading multiple files, keys for the join axis. If the
        axis already exists in the dataset, the concatenated dataset/dimarray
        will be re-indexed along the provided key, otherwise the keys
        will be used to create a new axis for stacking. In the latter case,
        keys' length needs to exactly match the number of input files, and
        if not provided, file names will be taken instead. Note you may
        manually rename the axes later, or use the `set_axis` method.

    align : bool, optional
        When reading multiple files, passed to `stack` (new axis) or
        `concatenate` (existing axis) to reindex all arrays onto common axes.
        (in `concatenate` mode, the concatenation axis is *not* re-indexed of course,
        only the secondary axes)
        Default to False.

    **kwargs : optional key-word arguments passed to align, if align is True
        When reading multiple files, passed to `stack` (new axis) or
        This includes: `sort` (False by default) and `join` ('outer' by default)

    Returns
    -------
    obj : DimArray or Dataset
        depending on whether a (single) variable name is passed as argument (names) or not

    See Also
    --------
    DatasetOnDisk.read, stack, concatenate, stack_ds, concatenate_ds, align,
    DimArray.write_nc, Dataset.write_nc

    Examples
    --------
    >>> import os
    >>> from dimarray import read_nc, get_datadir

    Single netCDF file

    >>> ncfile = os.path.join(get_datadir(), 'cmip5.CSIRO-Mk3-6-0.nc')

    >>> data = read_nc(ncfile)  # load full file
    >>> data
    Dataset of 2 variables
    0 / time (451): 1850 to 2300
    1 / scenario (5): 'historical' to 'rcp85'
    tsl: ('time', 'scenario')
    temp: ('time', 'scenario')
    >>> data = read_nc(ncfile,'temp') # only one variable
    >>> data = read_nc(ncfile,'temp', indices={"time":slice(2000,2100), "scenario":"rcp45"})  # load only a chunck of the data
    >>> data = read_nc(ncfile,'temp', indices={"time":1950.3}, tol=0.5)  #  approximate matching, adjust tolerance
    >>> data = read_nc(ncfile,'temp', indices={"time":-1}, indexing='position')  #  integer position indexing

    Multiple files
    Read variable 'temp' across multiple files (representing various climate models)
    In this case the variable is a time series, whose length may
    vary across experiments (thus align=True is passed to reindex axes before stacking)

    >>> direc = get_datadir()
    >>> temp = da.read_nc(direc+'/cmip5.*.nc', 'temp', align=True, axis='model')

    A new 'model' axis is created labeled with file names. It is then
    possible to rename it more appropriately, e.g. keeping only the part
    directly relevant to identify the experiment:

    >>> getmodel = lambda x: os.path.basename(x).split('.')[1] # extract model name from path
    >>> temp.set_axis(getmodel, axis='model') # would return a copy if inplace is not specified
    >>> temp
    dimarray: 9114 non-null elements (6671 null)
    0 / model (7): 'CSIRO-Mk3-6-0' to 'MPI-ESM-MR'
    1 / time (451): 1850 to 2300
    2 / scenario (5): 'historical' to 'rcp85'
    array(...)

    This works on datasets as well:

    >>> ds = da.read_nc(direc+'/cmip5.*.nc', align=True, axis='model')
    >>> ds.set_axis(getmodel, axis='model')
    >>> ds
    Dataset of 2 variables
    0 / model (7): 'CSIRO-Mk3-6-0' to 'MPI-ESM-MR'
    1 / time (451): 1850 to 2300
    2 / scenario (5): 'historical' to 'rcp85'
    tsl: ('model', 'time', 'scenario')
    temp: ('model', 'time', 'scenario')
    """
    # check for regular expression
    if type(f) is str:
        test = glob.glob(f)
        if len(test) == 1:
            f = test[0]
        elif len(test) == 0:
            raise ValueError('File is not present: '+repr(f))
        else:
            f = test
            f.sort()

    # multi-file ?
    if type(f) is list:
        mf = True
    else:
        mf = False

    # Read a single file
    if not mf:
        f, close = _maybe_open_file(f, mode='r')
        obj = DatasetOnDisk(f).read(names, *args, **kwargs)
        if close: f.close()

    # Read multiple files
    else:
        # single variable ==> DimArray (via Dataset)
        if names is not None and isinstance(names, str):
            obj = _read_multinc(f, [names], *args, **kwargs)
            obj = obj[names]

        # single variable ==> DimArray
        else:
            obj = _read_multinc(f, names, *args, **kwargs)

    return obj

def _read_multinc(fnames, names=None, axis=None, keys=None, align=False, sort=False, join='outer', concatenate_only=False, **kwargs):
    """ read multiple netCDF files

    Parameters
    ----------
    fnames : list of file names or file handles to be read
    names : variable names to be read
    axis : str, optional
        dimension along which the files are concatenated
        (created as new dimension if not already existing)
    keys : sequence, optional
        to be passed to stack_nc, if axis is not part of the dataset
    align : `bool`, optional
        if True, reindex axis prior to stacking / concatenating (default to False)
    sort : `bool`, optional
        if True, sort common axes prior to stacking / concatenating (defaut to False)
    join : `str`, optional
        join method in align (default to 'outer')
    concatenate_only : `bool`, optional
        if True, only concatenate along existing axis (and raise error if axis not existing)

    **kwargs : keyword arguments passed to DatasetOnDisk.read  (cannot
    contain 'axis', though, but indices can be passed as a dictionary
    if needed, e.g. {'time':2010})

    Returns
    -------
    dimarray's Dataset instance

    This function reads several files and call stack_ds or concatenate_ds
    depending on whether `axis` is absent from or present in the dataset, respectively

    See `dimarray.read_nc` for more complete documentation.

    See Also
    --------
    read_nc
    """
    variables = None
    dimensions = None

    datasets = []
    for fn in fnames:
        with DatasetOnDisk(fn) as f:
            ds = f.read(names, **kwargs)

        # check that the same variables are present in the file
        if variables is None:
            variables = ds.keys()
        else:
            variables_new = ds.keys()
            assert variables == variables_new, \
                    "netCDF files must contain the same \
                    subset of variables to be concatenated/stacked"

        # check that the same dimensions are present in the file
        if dimensions is None:
            dimensions = ds.dims
        else:
            dimensions_new = ds.dims
            assert dimensions == dimensions_new, \
                    "netCDF files must contain the same \
                    subset of dimensions to be concatenated/stacked"

        datasets.append(ds)

    # check that dimension in dataset if required
    if concatenate_only and axis not in dimensions:
        raise Exception('required axis {} not found, only got {}'.format(axis, dimensions))

    # Join dataset
    if axis in dimensions:
        ds = concatenate_ds(datasets, axis=axis, align=align, sort=sort, join=join)
        if keys is not None:
            ds = ds.reindex_axis(keys, axis=axis)
        # elif sort:
        #     ds = ds.sort_axis(axis=axis)
            #warnings.warn('keys argument will be ignored.', RuntimeWarning)

    else:
        # use file name as keys by default
        if keys is None:
            keys = [os.path.splitext(fn)[0] for fn in fnames]
        ds = stack_ds(datasets, axis=axis, keys=keys, align=align, sort=sort, join=join)

    return ds

#
# display summary information
#

def summary_nc(fname, name=None, metadata=False):
    """ Print summary information about the content of a netCDF file
    Deprecated, see dimarray.open_nc
    """
    warnings.warn("Deprecated. Use dimarray.open_nc", FutureWarning)
    with DatasetOnDisk(fname) as obj:
        if name is not None:
            obj = obj[name]
        print((obj.__repr__(metadata=metadata)))

def _maybe_open_file(f, mode='r', clobber=None, verbose=False, format=None):
    """ open a netCDF4 file

    Parameters
    ----------
    f : file name (str) or netCDF file handle
    mode: changed from original 'r','w','r' & clobber option:

    mode : `str`
        read or write access
        - 'r': read
        - 'w' : write, overwrite if file if present (clobber=True)
        - 'w-': create new file, but raise Exception if file is present (clobber=False)
        - 'a' : append, raise Exception if file is not present
        - 'a+': append if file is present, otherwise create

    format: passed to netCDF4.Dataset, only relevatn when mode = 'w', 'w-', 'a+'
        'NETCDF4', 'NETCDF4_CLASSIC', 'NETCDF3_CLASSIC', 'NETCDF3_64BIT'

    Returns
    -------
    f : netCDF file handle
    close: `bool`, `True` if input f indicated file name
    """
    format = format or get_option('io.nc.format')

    if mode == 'w-':
        mode = 'w'
        if clobber is None: clobber = False

    # mode 'a+' appends if file exists, otherwise create new variable
    elif mode == 'a+' and not isinstance(f, nc.Dataset):
        if os.path.exists(f): mode = 'a'
        else: mode = 'w'
        if clobber is None: clobber=False

    else:
        if clobber is None: clobber=True

    if not isinstance(f, nc.Dataset):
        fname = f

        # make sure the file does not exist if mode == "w"
        if os.path.exists(fname) and clobber and mode == "w":
            os.remove(fname)

        try:
            f = nc.Dataset(fname, mode, clobber=clobber, format=format)
        except UserWarning as msg:
            print(msg)
        except Exception as msg: # raise a weird RuntimeError
            #print "read from",fname
            raise IOError("{} => failed to opend {} in mode {}".format(msg, fname, mode)) # easier to handle

        if verbose:
            if 'r' in mode:
                print("read from",fname)
            elif 'a' in mode:
                print("append to",fname)
            else:
                print("write to",fname)
        close = True
    else:
        close = False # leave it open

    return f, close

max_stars_repo_name: perrette/dimarray

max_issues_repo_name: perrette/dimarray

max_forks_repo_name: perrette/dimarray


text:
from gensim import similarities
import os, logging
import numpy as np
import pandas as pd

class NeighborSampler():
    def __init__(self, chip_name, stage = 0, x_data=None):
        self.index = None
        self.chip_name = chip_name
        self.stage = stage
        self.build(x_data)

    def build(self, x_data=np.array([])):
        corpus = []
        for i in range(x_data.shape[0]):
            gensim_format_vec = []
            for j in range(x_data.shape[1]):
                gensim_format_vec.append((j, x_data[i][j]))
            corpus.append(gensim_format_vec)
        logging.info("#sample to build index: %s" % x_data.shape[0])
        self.get_index(corpus,x_data.shape[1])

    def get_index(self, corpus=None, n_feature=None):
        self.index = similarities.Similarity("%s_%s_neighbor.index.tmp" % (self.chip_name, self.stage), corpus, num_features=n_feature)
        return self.index

    def get_topk(self, vec, k):
        gensim_format_vec = []
        for i in range(len(vec)):
            gensim_format_vec.append((i, vec[i]))
        sims = self.index[gensim_format_vec]

        sim_sort = sorted(list(enumerate(sims)), key=lambda item: item[1])
        top_k = sim_sort[0:k]
        return top_k

def build_sampled_coexpression_matrix(x_data, filename, k = 1000):
    logging.info("Co-expression matrix building process starts.")
    logging.info("data shape: %s %s" % (x_data.shape[0], x_data.shape[1]))
    n_gene = x_data.shape[1]

    neighbor_sampler = NeighborSampler(filename, x_data = x_data)
    with open("matrix.txt", "w") as f:
        for i in range(n_gene):
            list_neighbor = neighbor_sampler.get_topk(x_data[i], k)
            list_neighbor.extend(neighbor_sampler.get_topk(x_data[i], k))
            list_neighbor = list(set(list_neighbor))
            logging.info("sample %s's topk neighbor: %s" % (i, list_neighbor))
            for j, value in list_neighbor:
                pearson_value = pearson(x_data[i], x_data[j])
                f.write("%s\t%s\t%s\n" % (i, j, pearson_value))
            f.flush()
            os.fsync(f.fileno())

def pearson(X, Y):
    return np.corrcoef(X, Y)[0][1]

if __name__ == "__main__":
    x_data = pd.read_csv("F:\gene2\spartan2\live-tutorials\inputData\example.csv", sep=",", header=None)
    build_sampled_coexpression_matrix(x_data, "out")


text:
#!/usr/bin/env python
import numpy as np

def intrinsic_eval_euclid(model, space, testX, testY):
    retval = np.zeros(testX.shape[0])
    denom = float(space.matrix.shape[0])
    testY = testY[:,0]
    pred = model.predict(testX)
    for i in xrange(testX.shape[0]):
        distances = np.sqrt(np.sum(np.square(space.matrix - pred[i]), axis=1))
        idx = testY[i]
        gold = distances[idx]
        rank = np.sum(distances <= gold) / denom
        retval[i] = rank
    return retval.mean()

from sklearn.preprocessing import normalize
def intrinsic_eval(model, space, testX, testY):
    space = space.normalize()
    retval = np.zeros(testX.shape[0])
    denom = float(space.matrix.shape[0])
    testY = testY[:,0]
    pred = model.predict(testX)
    pred = normalize(pred, norm='l2', axis=1)
    #pred += space.matrix[testY]
    pred = normalize(pred, norm='l2', axis=1)
    dots = pred.dot(space.matrix.T)
    goldscores = np.array([dots[i,yid] for i, yid in enumerate(testY)])
    greaterthan = np.mean(dots >= np.array([goldscores]).T, axis=1)
    return np.mean(greaterthan)


text:
import numpy as np

from src.setting import *
from src import gendf
from src.binary_io import NdlBinary

class MF6Like(gendf.MF6Like):
    def __init__(self, target_tape, equiprob_map, mf):
        self._mf = mf
        self.target_tape = target_tape
        self.equiprob_map = equiprob_map
        if mf == 27:
            self.prob_map = equiprob_map
        else:
            self.prob_map = equiprob_map[:,:1]

    def setFromMatrix(self, matrix): # override method
        raise AttributeError("'MF6Like' object has no attribute 'setFromMatrix'")

    def genEquiProbMap(self): #override method
        raise AttributeError("'MF6Like' object has no attribute 'genEquiProbMap'")

    def sampling(self, inc_group):
        if self._mf == 27:
            # sampling exit channel group
            line_pointer, group_pointer = self.target_tape[inc_group]
            if line_pointer < 0:
                return []
            rand = np.random.random()
            while True:
                if rand < self.equiprob_map[line_pointer]:
                    break
                line_pointer += 1
                group_pointer += 1

            return [[self._mf, group_pointer, -2]]

        else:
            # sampling exit channel group
            line_pointer, group_pointer = self.target_tape[inc_group]
            if line_pointer < 0:
                return []
            rand = np.random.random()
            while True:
                if rand < self.equiprob_map[line_pointer,0]:
                    break
                line_pointer += 1
                group_pointer += 1

            # sampling directional cosine
            rand = np.random.randint(0,int(ENV["equiprob_nbin"]))
            rand2 = np.random.random()
            mu = self.equiprob_map[line_pointer,rand+1] * rand2 \
                +self.equiprob_map[line_pointer,rand+2] * (1-rand2)

            if group_pointer > 187:
                raise ValueError

            return [[self._mf, group_pointer, mu]]

class MF16(gendf.MF16):
    def __init__(self, target_tape, prob_map):
        self._mf = 16
        self.target_tape = target_tape
        self.prob_map = prob_map

    def sampling(self, inc_group):
        # sampling the number of photon
        line_start, multiplicity, group_start = self.target_tape[inc_group]
        if line_start < 0:
            return []
        multiplicity = np.frombuffer(self.target_tape[inc_group,1:2].tobytes(),
                                     dtype=np.float32)[0]
        photon = []
        while True:
            rand = np.random.random()
            if rand > multiplicity:
                break
            line_pointer = line_start
            group_pointer = group_start
            # sampling exit channel group
            rand = np.random.random()
            while True:
                if rand < self.prob_map[line_pointer]:
                    break
                line_pointer += 1
                group_pointer += 1

            # sampling directional cosine
            mu = np.random.random() * 2 - 1
            photon += [[self._mf, group_pointer, mu]]
            multiplicity -= 1
        return photon

class Reaction(gendf.Reaction):
    def __init__(self, mt, sampling_rule):
        self.mt = mt
        self.mf = {}
        self._sampling_rule = sampling_rule

    def _add(self, egn, mf, nz, lrflag, data, label):
        raise AttributeError("'Reaction' object has no attribute '_add'")

    def sampling(self, inc_group):
        exit_particles = []
        for i in range(0, len(self._sampling_rule), 2):
            if self._sampling_rule[i+1] < 0:
                continue
            particle_temp = self.mf[self._sampling_rule[i]].sampling(inc_group)
            if len(particle_temp) > 0:
                exit_particles += particle_temp
        return exit_particles

class GPUNDL(gendf.GendfInterface):
    def __init__(self, file_name, verbose=False):
        super().__init__()

        file = NdlBinary(file_name, mode="r")
        self.reactions = {}

        # read atomic mass
        self.za = file.read()[0]
        self.mass = file.read()[0]

        # read total xs
        self.reactions[1] = Reaction(1, None)
        self.reactions[1].mf[3] = gendf.MF3(file.read())

        # read MT table & mt sampling table
        self._mt_target = file.read()
        self._mt_map = file.read()

        xs_prob_cumul = np.zeros_like(self.reactions[1].mf[3].xs)
        # read MT data
        for mt in self._mt_target:
            sampling_rule = file.read()
            self.reactions[mt] = Reaction(mt, sampling_rule)
            xs_prob = self._mt_map[:,np.argmax(self._mt_target == mt)] - xs_prob_cumul
            self.reactions[mt].mf[3] = gendf.MF3(xs_prob * self.reactions[1].mf[3].xs)
            xs_prob_cumul = self._mt_map[:,np.argmax(self._mt_target == mt)]
            tindex = 0
            for i in range(0, len(sampling_rule), 2):
                if sampling_rule[i+1] < tindex:
                    continue
                mf = sampling_rule[i]
                target_tape = file.read()
                prob_map = file.read()
                if mf == 16: # photon production
                    self.reactions[mt].mf[mf] = MF16(target_tape, prob_map)
                else:
                    self.reactions[mt].mf[mf] = MF6Like(target_tape, prob_map, mf)
                tindex += 1

    def getNeutronEnergyGroup(self, file_name):
        self.egn = np.load(file_name)

    def getPhotonEnergyGroup(self, file_name):
        self.egg = np.load(file_name)

    def sampling(self, inc_group):
        # get type of reaction
        rand = np.random.random()
        pointer = 0
        while True:
            if rand < self._mt_map[inc_group,pointer]:
                break
            pointer += 1
        mt = self._mt_target[pointer]
        return self.reactions[mt].sampling(inc_group)


text:
import gym
import numpy as np
from gym.spaces import Discrete, Dict, Box
from compiler_gym.envs.llvm.llvm_env import make_benchmark
import time
import sqlite3

class csr_rl:
    """A :class:
            This task is concerned with determining the LLVM passes and their order to minimize the code size.
            Following the setup of CompilerGym, we compute the code size reduction by measuring the ratio of
            LLVM IR instruction count reduction with respect to the LLVM -Oz code size optimization option.
            This metric is platform-independent and deterministic.

        Source:
            This environment is following the setup of CompilerGym(https://github.com/facebookresearch/CompilerGym).

        Observation:
            Type: Box(0, 9223372036854775807, (56,), int64)
            The Autophase observation space is a 56-dimension integer feature vector summarizing the static LLVM-IR representation.
            It is described in 'AutoPhase: Juggling HLS phase orderings in random forests with deep reinforcement learning'.
            paper link: https://proceedings.mlsys.org/paper/2020/file/4e732ced3463d06de0ca9a15b6153677-Paper.pdf

        Actions:
            Type: Discrete(123)
            Num      Action                         Description
            0        -add-discriminators            Add DWARF path discriminators.
            1        -adce                          Aggressive Dead Code Elimination.
            2        -aggressive-instcombine        Combine pattern based expressions.
            3        -alignment-from-assumptions    Replaces an instruction's opcode with a new one that takes operands of the same type.
            4        -always-inline                 Inliner for always_inline functions.
            ... ...

            120      -strip                         Strip all symbols from a module.
            121      -tailcallelim                  Tail Call Elimination.
            122      -mergereturn                   Unify function exit nodes.

        Reward:
            In all cases, lower code size is better. We compute the code size reduction by measuring the ratio of
            LLVM IR instruction count reduction with respect to the LLVM -Oz code size optimization option.

        Starting State:
            All observations are assigned a uniform random value in [-1..1]

        """

    def __init__(self, env_config):
        """ Defines the reinforcement leaning environment. Initialise with an environment.

                    :param env_config: including  "state_function", "action_function", "reward_function", "observation_space"
                """
        # self.init_from_env_config(env_config)
        self.sql_path = env_config.get("sql_path")
        self.benchmarks = env_config.get("benchmark")
        self.log_path = env_config.get("log_path")
        self.pass_path = env_config.get("pass_path")
        self.deadline = env_config.get("deadline")
        # print(self.benchmarks)
        self.env = gym.make(
            "llvm-autophase-ic-v0",  # selects the compiler to use
            benchmark=make_benchmark(
                [self.benchmarks]
            ),  # selects the program to compile
        )

        self.action_space = Discrete(self.env.action_space.n)

        self.observation_space = Dict(
            {
                "obs": self.env.observation_space,
                "action_mask": Box(low=0, high=1, shape=(self.action_space.n,)),
            }
        )

        self.running_reward = 0
        self.episode_reward = 0

        # record best
        self.view = 0
        self.best_actions = []
        self.time = time.time()

        self.num = 0
        self.seeds = env_config.get("seed")
        self.flag = False
        self.best_command = ""
        self.env.seed(env_config.get("seed"))
        self.score = 0
        # self.env.init()
        print(f"seed:{self.seeds}")

    def reset(self):
        """ reset the RL environment.
                        """
        self.running_reward = 0
        self.episode_reward = 0
        return {
            "obs": self.env.reset(),
            "action_mask": np.ones((self.env.action_space.n,), dtype=int),
        }

    def step(self, action):
        """Take a step.

        :param action: An action, or a sequence of actions. When multiple
                actions are provided the observation and reward are returned after
                running all of the actions.
        :return: A tuple of observation, observation_mask, score, done, and info.
        """

        used_time = time.time() - self.time
        benchmarks = self.benchmarks.split("/")[-1]
        with open(f"{self.pass_path}/{benchmarks}_pass.txt", "a+") as f:
            f.write("=================Step======================\r\n")
            f.write(f"step:{self.step}\r\nactions:{self.env.actions}\r\n")
            f.write(f"command line:{self.env.commandline()}\r\n")
            f.write(
                f"time={used_time:.3},action_len={len(self.env.actions)},episode_reward={self.env.episode_reward:.3%},reward_all={self.view:.3%}\r\n"
            )
        try:
            conn = sqlite3.connect(
                self.sql_path
            )
            c = conn.cursor()
            sql = "INSERT or IGNORE INTO CSR (TIME, BENCHMARK, RESULT, REWARD) \
                            VALUES (?, ?, ?, ?)"
            # print("time",round(time.time() - self.time))
            c.execute(
                sql,
                (
                    time.time() - self.time,
                    benchmarks,
                    self.env.commandline(),
                    self.env.episode_reward,
                ),
            )
            conn.commit()
            conn.close()

        except Exception as e:
            print(e)
        obs, rew, done, info = self.env.step(action)
        # print("action",action)
        self.running_reward += rew
        self.episode_reward = self.episode_reward + rew

        if self.env.episode_reward > self.view:
            self.view = self.env.episode_reward
            self.best_actions = self.env.actions
            self.best_command = self.env.commandline()

        if len(self.env.actions) > self.env.action_space.n:
            used_time = time.time() - self.time

            if used_time > self.deadline and not self.flag:
                self.flag = True
                with open(f"{self.log_path}/{benchmarks}_result_ppo.txt", "a+") as f:
                    # pass_time = time.time() - self.time
                    f.write(
                        f"actions{str(self.best_actions)}\r\ncommand_line:{self.best_command}\r\n"
                    )

                    f.write(
                        f"seeds:{self.seeds}\r\nreward: " + str(self.view) + "\r\n "
                    )
                    f.close()

            done = True
            self.score = self.running_reward if done else 0
            # print(f"time={used_time},quality_finally={self.view:.3%}")
            # print("action_masks",np.ones((self.env.action_space.n,)))

        return (
            {"obs": obs, "action_mask": np.ones((self.env.action_space.n,), dtype=int)},
            self.score,
            done,
            info,
        )

    def set_state(self, state):
        """ Set policy to specific state and action mask.

        :param state: Current reward and environments
        :return: state and action mask
        """
        self.running_reward = state[1]
        self.env = state[0].fork()
        obs = self.env.observation["Autophase"]
        return {
            "obs": obs,
            "action_mask": np.ones((self.env.action_space.n,), dtype=int),
        }

    def get_state(self):
        """Returns actor state.

        :return: current environment and reward
        """
        env_1 = self.env.fork()
        return env_1, self.running_reward

max_stars_repo_name: HuantWang/SUPERSONIC

max_issues_repo_name: HuantWang/SUPERSONIC

max_forks_repo_name: HuantWang/SUPERSONIC