text:
import numpy as np
from time import time as t

from downscale.utils.utils_func import change_dtype_if_required
from downscale.utils.decorators import print_func_executed_decorator, timer_decorator, change_dtype_if_required_decorator
from downscale.utils.context_managers import timer_context

try:
    from numba import jit, prange, float64, float32, int32, int64

    _numba = True
except ModuleNotFoundError:
    _numba = False

try:
    import numexpr as ne

    _numexpr = True
except ModuleNotFoundError:
    _numexpr = False

class Wind_utils:

    def __init__(self):
        pass

    @staticmethod
    def apply_log_profile(z_in=None, z_out=None, wind_in=None, z0=None, library="numexpr",
                          verbose=True, z_in_verbose=None, z_out_verbose=None):

        """
        Apply log profile to a wind time serie.

        Parameters
        ----------
        z_in : Input elevation
        z_out : Output elevation
        wind_in : Input wind
        z0 : Roughness length for momentum
        verbose: boolean
        z_in_verbose: str
        z_out_verbose: str

        Returns
        -------
        wind_out: Output wind, after logarithmic profile
        """

        print(f"Applied log profile: {z_in_verbose} => {z_out_verbose}") if verbose else None

        if _numexpr and library == "numexpr":
            return ne.evaluate("(wind_in * log(z_out / z0)) / log(z_in / z0)")
        else:
            return (wind_in * np.log(z_out / z0)) / np.log(z_in / z0)

    def _exposed_wind_speed_xarray(self, xarray_data, z_out=10):
        dims = xarray_data.Wind.dims
        return xarray_data.assign(exp_Wind=(dims,
                                            self._exposed_wind_speed_num(wind_speed=xarray_data.Wind.values,
                                                                         z_out=z_out,
                                                                         z0=xarray_data.Z0.values,
                                                                         z0_rel=xarray_data.Z0REL.values)))

    def _exposed_wind_speed_num(self, wind_speed=None, z_out=None, z0=None, z0_rel=None, library="numexpr"):
        if _numexpr and library == "numexpr":
            acceleration_factor = ne.evaluate(
                "log((z_out) / z0) * (z0 / (z0_rel+z0))**0.0706 / (log((z_out) / (z0_rel+z0)))")
            acceleration_factor = ne.evaluate("where(acceleration_factor > 0, acceleration_factor, 1)")
            exp_Wind = ne.evaluate("wind_speed * acceleration_factor")

        else:
            acceleration_factor = np.log(z_out / z0) * (z0 / (z0_rel + z0)) ** 0.0706 / (
                np.log(z_out / (z0_rel + z0)))
            acceleration_factor = np.where(acceleration_factor > 0, acceleration_factor, 1)
            exp_Wind = wind_speed * acceleration_factor

        return exp_Wind, acceleration_factor

    @print_func_executed_decorator("expose wind speed", level_begin="\n__", level_end="__")
    @timer_decorator("expose wind speed", unit="second", level=". . . . ")
    def exposed_wind_speed(self, wind_speed=None, z_out=None, z0=None, z0_rel=None,
                           library="numexpr", xarray_data=None, level_verbose="____", verbose=True):
        """
        Expose wind speed (deoperate subgrid parameterization from NWP)

        Parameters
        ----------
        wind_speed : ndarray
            Wind speed from NWP
        z_out : ndarray
            Elevation at which desinfluencing is performed
        z0 : ndarray
            Roughness length for momentum
        z0_rel : ndarray
            Roughness length for momentum associated with mean topography
        apply_directly_to_nwp : boolean, optional
            If True, updates the nwp directly by adding a new variable (Default: False)
        library: str
        verbose: boolean

        Returns
        -------
        exp_Wind: Unexposed wind
        acceleration_factor: Acceleration related to unexposition (usually >= 1)
        """

        if library == "xarray":
            return self._exposed_wind_speed_xarray(xarray_data, z_out=z_out)

        else:
            return self._exposed_wind_speed_num(wind_speed=wind_speed, z_out=z_out, z0=z0,
                                                z0_rel=z0_rel, library=library)

    @staticmethod
    def _scale_wind_for_ideal_case(wind_speed, wind_dir, input_speed, input_dir, verbose=True):
        """
        Specify specific wind speed and direction for ideal cases.

        Parameters
        ----------
        wind_speed : array
        wind_dir : array
        input_speed : constant
        input_dir : constant

        Returns
        -------
        wind_speed : array
            Constant wind speed.
        wind_direction : array
            Constant wind direction.
        """
        wind_speed = np.full(wind_speed.shape, input_speed)
        wind_dir = np.full(wind_dir.shape, input_dir)
        print("__Wind speed scaled for ideal cases") if verbose else None
        return wind_speed, wind_dir

    @print_func_executed_decorator("wind speed ratio (acceleration)", level_begin="\n__", level_end="__")
    @timer_decorator("wind speed ratio (acceleration)", unit="second", level=". . . . ")
    def wind_speed_ratio(self, num=None, den=None, library="numexpr", verbose=True):
        if _numexpr and library == "numexpr":
            print("____Library: Numexpr") if verbose else None
            a1 = ne.evaluate("where(den > 0, num / den, 1)")
        else:
            print("____Library: Numpy") if verbose else None
            a1 = np.where(den > 0, num / den, 1)
        return a1

    def _3D_wind_speed(self, U=None, V=None, W=None, out=None, library="numexpr", verbose_level="____", verbose=True):
        """
        UVW = np.sqrt(U**2 + V**2 + W**2)

        Parameters
        ----------
        U : dnarray
            Horizontal wind speed component U
        V : string, optional
            Horizontal wind speed component V
        W : string, optional
            Vertical wind speed component W
        out : ndarray, optional
            If specified, the output of the caluclation is directly written in out,
            which is best for memory consumption (Default: None)

        Returns
        -------
        UVW : ndarray
            Wind speed
        """

        if verbose:
            t0 = t()

        if out is None:
            if _numexpr and library == "numexpr":
                print("____Library: Numexpr") if verbose else None
                wind_speed = ne.evaluate("sqrt(U**2 + V**2 + W**2)")
            else:
                print("____Library: Numpy") if verbose else None
                wind_speed = np.sqrt(U ** 2 + V ** 2 + W ** 2)

            if verbose:
                print(f"{verbose_level}Wind speed computed")
                print(f". . . . Time to calculate wind speed: {np.round(t() - t0)} seconds")

            return wind_speed

        else:
            if _numexpr and library == "numexpr":
                ne.evaluate("sqrt(U**2 + V**2 + W**2)", out=out)
            else:
                np.sqrt(U ** 2 + V ** 2 + W ** 2, out=out)

            if verbose:
                print(f"{verbose_level}Wind speed computed")
                print(f". . . . Time to calculate wind speed: {np.round(t()-t0)}")

    def _2D_wind_speed(self, U=None, V=None, out=None, library="numexpr", verbose_level="____", verbose=True):
        """
        UV = np.sqrt(U**2 + V**2)

        Parameters
        ----------
        U : dnarray
            Horizontal wind speed component U
        V : string, optional
            Horizontal wind speed component V
        out : ndarray, optional
            If specified, the output of the caluclation is directly written in out,
            which is best for memory consumption (Default: None)

        Returns
        -------
        UV : ndarray
            Wind speed
        """

        if verbose:
            t0 = t()

        if out is None:
            if _numexpr and library == "numexpr":
                return ne.evaluate("sqrt(U**2 + V**2)")
            else:
                return np.sqrt(U ** 2 + V ** 2)

        else:
            if _numexpr and library == "numexpr":
                ne.evaluate("sqrt(U**2 + V**2)", out=out)
            else:
                np.sqrt(U ** 2 + V ** 2, out=out)

            if verbose:
                print(f"{verbose_level}Wind speed computed")
                print(f". . . . Time to calculate wind speed: {np.round(t()-t0)}")

    def _compute_wind_speed_num(self, U=None, V=None, W=None, verbose=True):
        if W is None:
            wind_speed = self._2D_wind_speed(U=U, V=V, verbose=verbose)
        else:
            wind_speed = self._3D_wind_speed(U=U, V=V, W=W, verbose=verbose)
        return wind_speed

    @change_dtype_if_required_decorator(np.float32)
    def _compute_wind_speed_num_out(self, U=None, V=None, W=None, out=None, verbose=True):
        if W is None:
            self._2D_wind_speed(U=U, V=V, out=out, verbose=verbose)
        else:
            self._3D_wind_speed(U=U, V=V, W=W, out=out, verbose=verbose)

    def compute_wind_speed_xarray(self, xarray_data, u_name="U", v_name="V", verbose=True):
        dims = xarray_data[u_name].dims
        return xarray_data.assign(Wind=(dims, self._compute_wind_speed_num(xarray_data[u_name].values,
                                                                           xarray_data[v_name].values,
                                                                      verbose=verbose)))

    def compute_wind_speed(self, U=None, V=None, W=None,
                           library='num', out=None,
                           xarray_data=None, u_name="U", v_name="V",
                           verbose=True, time_level=". . . . ", name_to_print="compute wind speed", unit="second"):
        """
        Calculates wind speed from wind speed components.

        First detects the number of wind component then calculates wind speed.
        The calculation can be performed on numexpr, numpy or xarray dataset

        Parameters
        ----------
        U : dnarray
            Horizontal wind speed component U
        V : string, optional
            Horizontal wind speed component V
        W : string, optional
            Vertical wind speed component V
        out : ndarray, optional
            If specified, the output of the calculation is directly written in out,
            which is best for memory consumption (Default: None)
        library : str, optional
            Select the librairie to use for calculation. If 'num' first test numexpr, if not available select numpy.
            If 'xarray', a xarray dataframe needs to be specified with the corresponding names for wind components.
            (Default: 'num')

        Returns
        -------
        UV : ndarray
            Wind speed
        """
        # Numexpr or numpy
        with timer_context(name_to_print, level=time_level, unit=unit, verbose=verbose):
            if library == 'num':
                if out is None:
                    return self._compute_wind_speed_num(U=U, V=V, W=W, verbose=verbose)
                else:
                    self._compute_wind_speed_num_out(U=U, V=V, W=W, out=out, verbose=verbose)

            if library == 'xarray':
                return self.compute_wind_speed_xarray(xarray_data, u_name=u_name, v_name=v_name, verbose=verbose)

    def compute_speed_and_direction_xarray(self, xarray_data, u_name="U", v_name="V", verbose=True):

        assert u_name in xarray_data
        assert v_name in xarray_data

        xarray_data = self.compute_wind_speed(library="xarray", xarray_data=xarray_data,
                                              u_name=u_name, v_name=v_name, verbose=verbose)
        xarray_data = self.direction_from_u_and_v(library="xarray", xarray_data=xarray_data,
                                                  u_name=u_name, v_name=v_name, verbose=verbose)
        print("____Wind and Wind_DIR calculated on xarray") if verbose else None

        return xarray_data

    @print_func_executed_decorator("wind speed scaling", level_begin="\n__", level_end="__")
    @timer_decorator("wind speed scaling", unit="second", level=". . . . ")
    def wind_speed_scaling(self, scaling_wind, prediction, type_scaling="linear", library="numexpr", verbose=True):
        """
        Linear: scaling_wind * prediction / 3

        Decrease acceleration and deceleration:
        Arctan_30_1: 30*np.arctan(scaling_wind/30) * prediction / 3

        Decrease acceleration only:
        Arctan_30_2: 30*np.arctan((scaling_wind* prediction / 3)/30)

        Parameters
        ----------
        scaling_wind : ndarray
            Scaling wind (ex: NWP wind)
        prediction : ndarray
            CNN ouptut
        linear : boolean, optional
            Linear scaling (Default: True)

        Returns
        -------
        prediction : ndarray
            Scaled wind
        """
        scaling_wind = scaling_wind.astype(np.float32)
        prediction = prediction.astype(np.float32)
        thirty = np.float32(30)
        three = np.float32(3)

        if type_scaling == "linear":
            if _numexpr and library == "numexpr":
                print("____Library: Numexpr")
                prediction = ne.evaluate("scaling_wind * prediction / 3")
            else:
                print("____Library: Numpy")
                prediction = scaling_wind * prediction / 3

        if type_scaling == "Arctan_30_1":
            if _numexpr and library == "numexpr":
                prediction = ne.evaluate("30*arctan(scaling_wind/30) * prediction / 3")
            else:
                prediction = 30*np.arctan(scaling_wind/30) * prediction / 3

        if type_scaling == "Arctan_30_2":
            if _numexpr and library == "numexpr":
                print("____Library: Numexpr")
                prediction = ne.evaluate("thirty*arctan((scaling_wind * prediction / three)/thirty)")
            else:
                print("____Library: Numpy")
                prediction = 30*np.arctan((scaling_wind * prediction / 3)/30)

        if type_scaling == "Arctan_10_2":
            if _numexpr and library == "numexpr":
                prediction = ne.evaluate("10*arctan((scaling_wind * prediction / 3)/10)")
            else:
                prediction = 10*np.arctan((scaling_wind * prediction / 3)/10)

        if type_scaling == "Arctan_20_2":
            if _numexpr and library == "numexpr":
                prediction = ne.evaluate("30*arctan((scaling_wind * prediction / 3)/30)")
            else:
                prediction = 20*np.arctan((scaling_wind * prediction / 3)/20)

        if type_scaling == "Arctan_38_2_2":
            constant = np.float32(38.2)
            scaling_wind = change_dtype_if_required(scaling_wind, np.float32)
            prediction = change_dtype_if_required(prediction, np.float32)
            if _numexpr and library == "numexpr":
                prediction = ne.evaluate("constant*arctan((scaling_wind * prediction / three)/constant)")
            else:
                prediction = constant*np.arctan((scaling_wind * prediction / three)/constant)

        if type_scaling == "Arctan_40_2":
            if _numexpr and library == "numexpr":
                prediction = ne.evaluate("30*arctan((scaling_wind * prediction / 3)/30)")
            else:
                prediction = 40*np.arctan((scaling_wind * prediction / 3)/40)

        if type_scaling == "Arctan_50_2":
            if _numexpr and library == "numexpr":
                prediction = ne.evaluate("30*arctan((scaling_wind * prediction / 3)/30)")
            else:
                prediction = 50*np.arctan((scaling_wind * prediction / 3)/50)

        prediction = change_dtype_if_required(prediction, np.float32)
        print(f"____Type {type_scaling}") if verbose else None
        return prediction

    @change_dtype_if_required_decorator(np.float32)
    @print_func_executed_decorator("computing angular deviation", level_begin="____", level_end="____")
    @timer_decorator("computing angular deviation", unit="second", level=". . . . ")
    def angular_deviation(self, U, V, unit_output="radian", library="numexpr", verbose=True):
        """
        Angular deviation from incoming flow.

        The incoming flow is supposed from the West so that V=0. If deviated, V != 0.
        The angular deviation is then np.arctan(V / U)

        Parameters
        ----------
        U : dnarray
            Horizontal wind speed component U
        V : string, optional
            Horizontal wind speed component V

        Returns
        -------
        alpha : ndarray
            Angular deviation [rad]
        """
        if _numexpr and library == "numexpr":
            print("____Library: Numexpr") if verbose else None
            alpha = ne.evaluate("where(U == 0, where(V == 0, 0, V/abs(V) * 3.14159 / 2), arctan(V / U))")
        else:
            print("____Library: Numpy") if verbose else None
            alpha = np.where(U == 0,
                             np.where(V == 0, 0, np.sign(V) * np.pi / 2),
                             np.arctan(V / U))

        if unit_output == "degree":
            alpha = self._rad_to_deg_num(alpha)

        return alpha

    @change_dtype_if_required_decorator(np.float32)
    @print_func_executed_decorator("computing wind direction from angular deviation", level_begin="____", level_end="____")
    @timer_decorator("computing wind direction from angular deviation", unit="second", level=". . . . ")
    def direction_from_alpha(self, wind_dir, alpha, unit_direction ="degree", unit_alpha="radian", unit_output="radian",
                             library="numexpr", verbose=True):
        """
        wind_dir - alpha

        Wind direction modified by angular deviation due to wind/topography interaction.
        Warning: this function might return a new wind direction in a rotated coordinates.

        Parameters
        ----------
        wind_dir : dnarray
            Initial NWP wind direction
        alpha : dnarray
            Angular deviation
        input_dir_in_degree : boolean, optionnal
            If True, converts the input wind direction from radans to degrees (Default: True)

        Returns
        -------
        alpha : ndarray
            Modified wind direction
        """
        if unit_direction == "degree":
            wind_dir = self._deg_to_rad_num(wind_dir, library="numexpr")

        if unit_alpha == "degree":
            alpha = self._deg_to_rad_num(alpha, library="numexpr")

        if _numexpr and library == "numexpr":
            print("____Library: Numexpr") if verbose else None
            UV_DIR = ne.evaluate("wind_dir - alpha")
        else:
            print("____Library: Numpy") if verbose else None
            UV_DIR = wind_dir - alpha

        if unit_output == "degree":
            UV_DIR = self._rad_to_deg_num(UV_DIR)

        return UV_DIR

    @change_dtype_if_required_decorator(np.float32)
    def _u_zonal_component(self, UV=None, UV_DIR=None, library="numexpr", unit_direction="degree", verbose=True):

        assert unit_direction == "radian"

        if _numexpr and library == "numexpr":
            print(f"____Library: {library}") if verbose else None
            return ne.evaluate("-sin(UV_DIR) * UV")
        else:
            print("____Library: numpy") if verbose else None
            return -np.sin(UV_DIR) * UV

    @change_dtype_if_required_decorator(np.float32)
    def _v_meridional_component(self, UV=None, UV_DIR=None, library="numexpr", unit_direction="degree", verbose=True):

        assert unit_direction == "radian"

        if _numexpr and library == "numexpr":
            print(f"____Library: {library}") if verbose else None
            return ne.evaluate("-cos(UV_DIR) * UV")
        else:
            print("____Library: numpy") if verbose else None
            return -np.cos(UV_DIR) * UV

    @print_func_executed_decorator("degree to radians", level_begin="________", level_end="________")
    @timer_decorator("converting degrees to radians", unit="second", level=". . . . ")
    @change_dtype_if_required_decorator(np.float32)
    def _deg_to_rad_num(self, UV_DIR, library="numexpr"):
        if _numexpr and library == "numexpr":
            return ne.evaluate("3.14159 * UV_DIR/180")
        else:
            return np.deg2rad(UV_DIR)

    @print_func_executed_decorator("radians to degree", level_begin="________", level_end="________")
    @timer_decorator("converting radians to degrees", unit="second", level=". . . . ")
    @change_dtype_if_required_decorator(np.float32)
    def _rad_to_deg_num(self, UV_DIR, library="numexpr"):
        if _numexpr and library == "numexpr":
            return ne.evaluate("180 * UV_DIR/3.14159")
        else:
            return np.rad2deg(UV_DIR)

    def _horizontal_wind_component_num(self, UV, UV_DIR, library="numexpr", unit_direction="degree", verbose=True):
        if unit_direction == "degree":
            UV_DIR = self._deg_to_rad_num(UV_DIR, library=library)
            unit_direction = "radian"
        U = self._u_zonal_component(UV=UV, UV_DIR=UV_DIR, library=library, unit_direction=unit_direction,
                                    verbose=verbose)
        V = self._v_meridional_component(UV=UV, UV_DIR=UV_DIR, library=library, unit_direction=unit_direction,
                                         verbose=verbose)
        return U, V

    def _horizontal_wind_component_xarray(self, xarray_data, wind_name="UV", wind_dir_name="UV_DIR",
                                          unit_direction="degree", verbose=True):
        if unit_direction == "degree":
            xarray_data = xarray_data.assign(theta=lambda x: (np.pi / 180) * (x[wind_dir_name] % 360))
        xarray_data = xarray_data.assign(U=lambda x: -x[wind_name] * np.sin(x["theta"]))
        xarray_data = xarray_data.assign(V=lambda x: -x[wind_name] * np.cos(x["theta"]))
        return xarray_data

    @print_func_executed_decorator("computing horizontal wind components", level_begin="____", level_end="____")
    @timer_decorator("computing horizontal wind components", unit="second", level=". . . . ")
    def horizontal_wind_component(self, UV=None, UV_DIR=None,
                                  library='numexpr', xarray_data=None, wind_name="Wind",
                                  wind_dir_name="Wind_DIR", unit_direction="radian", verbose=True):
        """
        U = -np.sin(UV_DIR) * UV
        V = -np.cos(UV_DIR) * UV

        Computes U and V component from wind direction and wind speed

        Parameters
        ----------
        UV : dnarray
            Wind speed
        UV_DIR : dnarray
            Wind_direction
        working_with_xarray : boolean, optionnal
            If True, computes U and V on an xarray dataframe (provided with names of variables)
            and returns the dataframe

        Returns
        -------
        U : dnarray
            Horizontal wind speed component U
        V : string, optional
            Horizontal wind speed component V

        References
        ----------
        Liston, G. E., & Elder, K. (2006). A meteorological distribution system for high-resolution terrestrial modeling (MicroMet). Journal of Hydrometeorology, 7(2), 217-234.
        """
        if library == "xarray":
            print("____Library: xarray") if verbose else None
            return self._horizontal_wind_component_xarray(xarray_data,
                                                          wind_name=wind_name,
                                                          wind_dir_name=wind_dir_name,
                                                          unit_direction=unit_direction,
                                                          verbose=False)
        else:
            return self._horizontal_wind_component_num(UV, UV_DIR, library=library,
                                                       unit_direction=unit_direction,
                                                       verbose=verbose)

    @change_dtype_if_required_decorator(np.float32)
    def _direction_from_u_and_v_num(self, U, V, library="numexpr", verbose=True):
        if _numexpr and library == "numexpr":
            print("____Library: Numexpr") if verbose else None
            return ne.evaluate("(180 + (180/3.14159)*arctan2(U,V)) % 360")
        else:
            print("____Library: Numpy") if verbose else None
            return np.mod(180 + np.rad2deg(np.arctan2(U, V)), 360)

    def _direction_from_u_and_v_xarray(self, xarray_data=None, u_name="U", v_name="V"):
        dims = xarray_data[u_name].dims
        return xarray_data.assign(Wind_DIR=(dims,
                                            self._direction_from_u_and_v_num(xarray_data[u_name].values,
                                                                             xarray_data[v_name].values,
                                                                             verbose=False)))

    def _direction_from_u_and_v_degree(self, U=None, V=None, library="numexpr", verbose=True,
                                       xarray_data=None, u_name="U", v_name="V"):
        if library == "xarray":
            print("____Library: xarray") if verbose else None
            return self._direction_from_u_and_v_xarray(xarray_data, u_name=u_name, v_name=v_name)
        else:
            return self._direction_from_u_and_v_num(U=U, V=V, library=library, verbose=verbose)

    def direction_from_u_and_v(self, U=None, V=None, unit_output="degree",
                               library="numexpr", xarray_data=None, u_name="U", v_name="V",
                               verbose=True, name_to_print="direction from u and v", time_level=". . . . ",
                               unit="second"):
        """
        UV_DIR = 180 + (180/3.14159)*arctan2(U,V)) % 360

        Compute wind direction from U and V components

        Parameters
        ----------
        U : dnarray
            Horizontal wind speed component U
        V : string, optional
            Horizontal wind speed component V
        output_in_degree : boolean, optionnal
            If True, the wind direction is converted to degrees

        Returns
        -------
        UV_DIR : dnarray
            Wind direction
        """
        with timer_context(name_to_print, level=time_level, unit=unit, verbose=verbose):
            if unit_output == "degree":
                return self._direction_from_u_and_v_degree(U=U, V=V, library=library,
                                                           xarray_data=xarray_data, u_name=u_name, v_name=v_name,
                                                           verbose=verbose)
            else:
                raise NotImplementedError

reason: import numpy,from numba,import numexpr


text:
import orekit
from math import radians, degrees, sqrt, pi
import datetime
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import os, random

orekit.initVM()

from org.orekit.frames import FramesFactory, Frame
from org.orekit.bodies import OneAxisEllipsoid
from org.orekit.orbits import KeplerianOrbit, Orbit
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.analytical import KeplerianPropagator
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.propagation.sampling import OrekitFixedStepHandler
from org.orekit.orbits import OrbitType, PositionAngle
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.utils import IERSConventions, Constants
from org.orekit.forces.maneuvers import ConstantThrustManeuver
from org.orekit.frames import LOFType
from org.orekit.attitudes import LofOffset
from org.orekit.python import PythonEventHandler, PythonOrekitFixedStepHandler
from orekit.pyhelpers import setup_orekit_curdir
from org.orekit.forces.gravity import NewtonianAttraction
from org.orekit.utils import Constants

from org.hipparchus.geometry.euclidean.threed import Vector3D
from java.util import Arrays
from orekit import JArray_double

setup_orekit_curdir()

FUEL_MASS = "Fuel Mass"

UTC = TimeScalesFactory.getUTC()
inertial_frame = FramesFactory.getEME2000()
attitude = LofOffset(inertial_frame, LOFType.LVLH)
EARTH_RADIUS = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
MU = Constants.WGS84_EARTH_MU

class OrekitEnv:
    """
    This class uses Orekit to create an environment to propagate a satellite
    """

    def __init__(self, state, state_targ, date, duration, mass, stepT):
        """
        initializes the orekit VM and included libraries
        Params:
        _prop: The propagation object
        _initial_date: The initial start date of the propagation
        _orbit: The orbit type (Keplerian, Circular, etc...)
        _currentDate: The current date during propagation
        _currentOrbit: The current orbit paramenters
        _px: spacecraft position in the x-direction
        _py: spacecraft position in the y-direction
        _sc_fuel: The spacecraft with fuel accounted for
        _extrap_Date: changing date during propagation state
        _sc_state: The spacecraft without fuel
        """

        self._prop = None
        self._initial_date = None
        self._orbit = None
        self._currentDate = None
        self._currentOrbit = None

        self.dry_mass = mass[0]
        self.fuel_mass = mass[1]
        self.cuf_fuel_mass = self.fuel_mass
        self.initial_mass = self.dry_mass + self.fuel_mass

        self.px = []
        self.py = []
        self.pz = []
        self.a_orbit = []
        self.ex_orbit = []
        self.ey_orbit = []
        self.hx_orbit = []
        self.hy_orbit = []
        self.lv_orbit = []

        self.adot_orbit = []
        self.exdot_orbit = []
        self.eydot_orbit = []
        self.hxdot_orbit = []
        self.hydot_orbit = []

        self._sc_fuel = None
        self._extrap_Date = None
        self._targetOrbit = None
        self.target_px = []
        self.target_py = []
        self.target_pz = []

        self._orbit_tolerance = {'a': 10000, 'ex': 0.01, 'ey': 0.01, 'hx': 0.001, 'hy': 0.001, 'lv': 0.01}

        self.randomize = False
        self._orbit_randomizer = {'a': 4000.0e3, 'e': 0.2, 'i': 2.0, 'w': 10.0, 'omega': 10.0, 'lv': 5.0}
        self.seed_state = state
        self.seed_target = state_targ
        self.target_hit = False

        self.set_date(date)
        self._extrap_Date = self._initial_date
        self.create_orbit(state, self._initial_date, target=False)
        self.set_spacecraft(self.initial_mass, self.cuf_fuel_mass)
        self.create_Propagator()
        self.setForceModel()
        self.final_date = self._initial_date.shiftedBy(duration)
        self.create_orbit(state_targ, self.final_date, target=True)

        self.stepT = stepT
        self.action_space = 3  # output thrust directions
        self.observation_space = 10  # states
        self.action_bound = 0.6  # Max thrust limit
        self._isp = 3100.0

        self.r_target_state = np.array(
            [self._targetOrbit.getA(), self._targetOrbit.getEquinoctialEx(), self._targetOrbit.getEquinoctialEy(),
             self._targetOrbit.getHx(), self._targetOrbit.getHy(), self._targetOrbit.getLv()])

        self.r_initial_state = np.array([self._orbit.getA(), self._orbit.getEquinoctialEx(), self._orbit.getEquinoctialEy(),
                                  self._orbit.getHx(), self._orbit.getHy(), self._orbit.getLv()])

        self.r_target_state = self.get_state(self._targetOrbit)
        self.r_initial_state = self.get_state(self._orbit)

    def set_date(self, date=None, absolute_date=None, step=0):
        """
        Set up the date for an orekit secnario
        :param date: list [year, month, day, hour, minute, second] (optional)
        :param absolute_date: a orekit Absolute Date object (optional)
        :param step: seconds to shift the date by (int)
        :return:
        """
        if date != None:
            year, month, day, hour, minute, sec = date
            self._initial_date = AbsoluteDate(year, month, day, hour, minute, sec, UTC)
        elif absolute_date != None and step != 0:
            self._extrap_Date = AbsoluteDate(absolute_date, step, UTC)
        else:
            # no argument given, use current date and time
            now = datetime.datetime.now()
            year, month, day, hour, minute, sec = now.year, now.month, now.day, now.hour, now.minute, float(now.second)
            self._initial_date = AbsoluteDate(year, month, day, hour, minute, sec, UTC)

    def create_orbit(self, state, date, target=False):
        """
         Crate the initial orbit using Keplarian elements
        :param state: a state list [a, e, i, omega, raan, lM]
        :param date: A date given as an orekit absolute date object
        :param target: a target orbit list [a, e, i, omega, raan, lM]
        :return:
        """
        a, e, i, omega, raan, lM = state

        a += EARTH_RADIUS
        i = radians(i)
        omega = radians(omega)
        raan = radians(raan)
        lM = radians(lM)

        aDot, eDot, iDot, paDot, rannDot, anomalyDot = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0

        # Set inertial frame
        set_orbit = KeplerianOrbit(a, e, i, omega, raan, lM,
                                   aDot, eDot, iDot, paDot, rannDot, anomalyDot,
                                   PositionAngle.TRUE, inertial_frame, date, MU)

        if target:
            self._targetOrbit = set_orbit
        else:
            self._currentOrbit = set_orbit
            self._orbit = set_orbit

    def convert_to_keplerian(self, orbit):

        ko = KeplerianOrbit(orbit)
        return ko

    def set_spacecraft(self, mass, fuel_mass):
        """
        Add the fuel mass to the spacecraft
        :param mass: dry mass of spacecraft (kg, flaot)
        :param fuel_mass:
        :return:
        """
        sc_state = SpacecraftState(self._orbit, mass)
        # print(f'{sc_state.getMass()}')
        self._sc_fuel = sc_state.addAdditionalState (FUEL_MASS, fuel_mass)

    def create_Propagator(self, prop_master_mode=False):
        """
        Creates and initializes the propagator
        :param prop_master_mode: Set propagator to slave of master mode (slave default)
        :return:
        """
        # tol = NumericalPropagator.tolerances(1.0, self._orbit, self._orbit.getType())
        minStep = 0.001
        maxStep = 500.0

        position_tolerance = 60.0
        tolerances = NumericalPropagator.tolerances(position_tolerance, self._orbit, self._orbit.getType())
        abs_tolerance = JArray_double.cast_(tolerances[0])
        rel_telerance = JArray_double.cast_(tolerances[1])

        # integrator = DormandPrince853Integrator(minStep, maxStep, 1e-5, 1e-10)
        integrator = DormandPrince853Integrator(minStep, maxStep, abs_tolerance, rel_telerance)

        integrator.setInitialStepSize(10.0)

        numProp = NumericalPropagator(integrator)
        numProp.setInitialState(self._sc_fuel)
        numProp.setMu(MU)
        numProp.setOrbitType(OrbitType.KEPLERIAN)

        if prop_master_mode:
            output_step = 5.0
            handler = OutputHandler()
            numProp.setMasterMode(output_step, handler)
        else:
            numProp.setSlaveMode()

        self._prop = numProp
        self._prop.setAttitudeProvider(attitude)

    def render_plots(self, episode=1, save=True, show=True):
        """
        Renders the x-y plots of the spacecraft trajectory
        :return:
        """
        if 0 <= episode < 10:
            episode = '00'+str(episode)
        elif 10 <= episode < 100:
            episode = '0'+str(episode)
        elif episode >= 100:
            episode = str(episode)

        os.makedirs('results/' + episode, exist_ok=True)
        save_path = 'results/' + episode + '/'

        plt.clf()
        oe_fig = plt.figure(1)
        oe_params = ('sma', 'e_x', 'e_y', 'h_x', 'h_y', 'lv')
        oe = [np.asarray(self.a_orbit)/1e3, self.ex_orbit, self.ey_orbit, self.hx_orbit, self.hy_orbit, self.lv_orbit]
        oe_target = [self.r_target_state[0]/1e3, self.r_target_state[1], self.r_target_state[2], self.r_target_state[3],
                   self.r_target_state[4], self.r_target_state[5]]
        # plt.plot(np.asarray(self.px) / 1000, np.asarray(self.py) / 1000, '-b',
        #          np.asarray(self.target_px)/1000, np.asarray(self.target_py)/1000, '-r')
        # plt.xlabel("x (km)")
        # plt.ylabel("y (km)")
        # plt.tight_layout()
        # if save:
        #     plt.savefig(save_path + '2d.pdf')
        # if show: plt.show()
        fig = plt.figure(2)
        ax = fig.gca(projection='3d')
        ax.plot(np.asarray(self.px)/1000, np.asarray(self.py)/1000, np.asarray(self.pz)/1000,label='Satellite Trajectory')
        ax.plot(np.asarray(self.target_px)/1000, np.asarray(self.target_py)/1000, np.asarray(self.target_pz)/1000,
                color='red',label='Target trajectory')
        ax.legend()
        ax.set_xlabel('X (km)')
        ax.set_ylabel('Y (km)')
        ax.set_zlabel('Z (km)')
        ax.set_zlim(-7000, 7000)
        # plt.title('Inclination change maneuver')
        plt.tight_layout()
        if save:
            plt.savefig(save_path + '3d.pdf')
        plt.figure(3)
        for i in range(len(oe_params)):
            plt.subplot(3,2,i+1)
            plt.plot(oe[i])
            plt.scatter(len(oe[i]), oe_target[i], c='red')
            plt.ylabel(oe_params[i])
        plt.tight_layout()
        if save:
            plt.savefig(save_path + 'oe.pdf')

        if save:
            np.save(save_path+'oe.npy', np.array(oe))
            np.save(save_path+'xyz.npy', np.array([self.px, self.py, self.pz]))

        if show: plt.show()
        ax.clear()

    def setForceModel(self):
        """
        Set up environment force models
        """
        # force model gravity field
        newattr = NewtonianAttraction(MU)
        itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)  # International Terrestrial Reference Frame, earth fixed

        earth = OneAxisEllipsoid(EARTH_RADIUS,
                                Constants.WGS84_EARTH_FLATTENING,itrf)
        gravityProvider = GravityFieldFactory.getNormalizedProvider(8, 8)
        self._prop.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider))

    def reset(self):
        """
        Resets the orekit enviornment
        :return:
        """

        self._prop = None
        self._currentDate = None
        self._currentOrbit = None

        # Randomizes the initial orbit
        if self.randomize:
            self._orbit = None
            a_rand = self.seed_state[0]
            e_rand = self.seed_state[1]
            w_rand = self.seed_state[3]
            omega_rand = self.seed_state[4]
            lv_rand = self.seed_state[5]

            a_rand = random.uniform(self.seed_state[0]-self._orbit_randomizer['a'], self._orbit_randomizer['a']+ self.seed_state[0])
            e_rand = random.uniform(self.seed_state[1]-self._orbit_randomizer['e'], self._orbit_randomizer['e']+ self.seed_state[1])
            i_rand = random.uniform(self.seed_state[2]-self._orbit_randomizer['i'], self._orbit_randomizer['i']+ self.seed_state[2])
            w_rand = random.uniform(self.seed_state[3]-self._orbit_randomizer['w'], self._orbit_randomizer['w']+ self.seed_state[3])
            omega_rand = random.uniform(self.seed_state[4]-self._orbit_randomizer['omega'], self._orbit_randomizer['omega']+ self.seed_state[4])
            lv_rand = random.uniform(self.seed_state[5]-self._orbit_randomizer['lv'], self._orbit_randomizer['lv']+ self.seed_state[5])
            state = [a_rand, e_rand, i_rand, w_rand, omega_rand, lv_rand]
            self.create_orbit(state, self._initial_date, target=False)
        else:
            self._currentOrbit = self._orbit

        self._currentDate = self._initial_date
        self._extrap_Date = self._initial_date

        self.set_spacecraft(self.initial_mass, self.fuel_mass)
        self.cuf_fuel_mass = self.fuel_mass
        self.create_Propagator()
        self.setForceModel()

        self.px = []
        self.py = []
        self.pz = []

        self.a_orbit = []
        self.ex_orbit = []
        self.ey_orbit = []
        self.hx_orbit = []
        self.hy_orbit = []
        self.lv_orbit = []

        self.adot_orbit = []
        self.exdot_orbit = []
        self.eydot_orbit = []
        self.hxdot_orbit = []
        self.hydot_orbit = []

        state = np.array([self._orbit.getA() / self.r_target_state[0],
                          self._orbit.getEquinoctialEx(),
                          self._orbit.getEquinoctialEy(),
                          self._orbit.getHx(),
                          self._orbit.getHy(), 0, 0, 0, 0, 0])
        return state

    @property
    def getTotalMass(self):
        """
        Get the total mass of the spacecraft
        :return: dry mass + fuel mass (kg)
        """
        return self._sc_fuel.getAdditionalState(FUEL_MASS)[0] + self._sc_fuel.getMass()

    def get_state(self, orbit, with_derivatives=True):

        if with_derivatives:
            state = [orbit.getA(), orbit.getEquinoctialEx(), orbit.getEquinoctialEy(),
                     orbit.getHx(), orbit.getHy(), orbit.getLv(),
                     orbit.getADot(), orbit.getEquinoctialExDot(),
                     orbit.getEquinoctialEyDot(),
                     orbit.getHxDot(), orbit.getHyDot(), orbit.getLvDot()]
        else:
            state = [orbit.getA(), orbit.getEquinoctialEx(), orbit.getEquinoctialEy(),
                       orbit.getHx(), orbit.getHy(), orbit.getLv()]

        return state

    def step(self, thrust):
        """
        Take a propagation step
        :param thrust_mag: Thrust magnitude (Newtons, float)
        :return: spacecraft state (np.array), reward value (float), don\
        e (bbol)
        """
        thrust_mag = np.linalg.norm(thrust)
        thrust_dir = thrust / thrust_mag
        DIRECTION = Vector3D(float(thrust_dir[0]), float(thrust_dir[1]), float(thrust_dir[2]))

        if thrust_mag <= 0:
            # DIRECTION = Vector3D.MINUS_J
            thrust_mag = abs(float(thrust_mag))
        else:
            # DIRECTION = Vector3D.PLUS_J
            thrust_mag = float(thrust_mag)

        thrust_force = ConstantThrustManeuver(self._extrap_Date, self.stepT, thrust_mag, self._isp, attitude, DIRECTION)
        # thrust_force.init(SpacecraftState(self._currentOrbit, float(self.mass + self.fuel_mass)), self._extrap_Date)
        self._prop.addForceModel(thrust_force)
        currentState = self._prop.propagate(self._extrap_Date.shiftedBy(self.stepT))
        # print(f'{currentState.getMass()}')
        self.cuf_fuel_mass = currentState.getMass() - self.dry_mass
        self._currentDate = currentState.getDate()
        self._extrap_Date = self._currentDate
        self._currentOrbit = currentState.getOrbit()
        coord = currentState.getPVCoordinates().getPosition()

        self.px.append(coord.getX())
        self.py.append(coord.getY())
        self.pz.append(coord.getZ())
        self.a_orbit.append(currentState.getA())
        self.ex_orbit.append(currentState.getEquinoctialEx())
        self.ey_orbit.append(currentState.getEquinoctialEy())
        self.hx_orbit.append(currentState.getHx())
        self.hy_orbit.append(currentState.getHy())
        self.lv_orbit.append(currentState.getLv())

        reward, done = self.dist_reward(thrust)

        state_1 = [(self._currentOrbit.getA()) / self.r_target_state[0],
                   self._currentOrbit.getEquinoctialEx(), self._currentOrbit.getEquinoctialEy(),
                   self._currentOrbit.getHx(), self._currentOrbit.getHy(),
                   self._currentOrbit.getADot(), self._currentOrbit.getEquinoctialExDot(),
                   self._currentOrbit.getEquinoctialEyDot(),
                   self._currentOrbit.getHxDot(), self._currentOrbit.getHyDot()
                   ]

        self.adot_orbit.append(self._currentOrbit.getADot())
        self.exdot_orbit.append(self._currentOrbit.getEquinoctialExDot())
        self.eydot_orbit.append(self._currentOrbit.getEquinoctialEyDot())
        self.hxdot_orbit.append(self._currentOrbit.getHxDot())
        self.hydot_orbit.append(self._currentOrbit.getHyDot())

        return state_1, reward, done

    def dist_reward(self, thrust):
        """
        Computes the reward based on the state of the agent
        :param thrust: Spacecraft thrust
        :return: reward value (float)
        """
        # a, ecc, i, w, omega, E, adot, edot, idot, wdot, omegadot, Edot = state

        done = False

        # thrust_mag = np.linalg.norm(thrust)

        state = np.array([self._currentOrbit.getA(), self._currentOrbit.getEquinoctialEx(), self._currentOrbit.getEquinoctialEy(),
                          self._currentOrbit.getHx(), self._currentOrbit.getHy(), self._currentOrbit.getLv()])

        # Inclination change reward
        reward_a = np.sqrt((self.r_target_state[0] - state[0])**2) / self.r_target_state[0]
        reward_ex = np.sqrt((self.r_target_state[1] - state[1])**2)
        reward_ey = np.sqrt((self.r_target_state[2] - state[2])**2)
        reward_hx = np.sqrt((self.r_target_state[3] - state[3])**2)
        reward_hy = np.sqrt((self.r_target_state[4] - state[4])**2)

        # reward = -(reward_a + reward_hx*10 + reward_hy*10 + reward_ex + reward_ey)

        reward = -(reward_a + reward_hx*10 + reward_hy*10 + reward_ex + reward_ey*10) + (self.cuf_fuel_mass/ self.fuel_mass)
        #current
        # Inclination change
        # reward = -(reward_a*10 + reward_hx*10 + reward_hy*10 + reward_ex*10 + reward_ey)
        # Terminal staes
        if abs(self.r_target_state[0] - state[0]) <= self._orbit_tolerance['a'] and \
           abs(self.r_target_state[1] - state[1]) <= self._orbit_tolerance['ex'] and \
           abs(self.r_target_state[2] - state[2]) <= self._orbit_tolerance['ey'] and \
           abs(self.r_target_state[3] - state[3]) <= self._orbit_tolerance['hx'] and \
           abs(self.r_target_state[4] - state[4]) <= self._orbit_tolerance['hy']:
            # self.final_date.durationFrom(self._extrap_Date) <= 360:
            reward += 1
            done = True
            print('hit')
            self.target_hit = True
            return reward, done

        # if self._sc_fuel.getAdditionalState(FUEL_MASS)[0] <= 0:
        if self.cuf_fuel_mass <= 0:
            print('Ran out of fuel')
            done = True
            reward = -1
            return reward, done

        if self._currentOrbit.getA() < EARTH_RADIUS:
            reward = -1
            done = True
            print('In earth')
            return reward, done

        if self._extrap_Date.compareTo(self.final_date) >= 0:
            reward = -1
            done = True

        return reward, done

    def render_target(self):

        target_sc = SpacecraftState(self._targetOrbit)

        # Orbit time for one orbit regardless of semi-major axis
        orbit_time = sqrt(4 * pi**2 * self._targetOrbit.getA()**3 / MU) * 2.0
        minStep = 1.e-3
        maxStep = 1.e+3

        integrator = DormandPrince853Integrator(minStep, maxStep, 1e-5, 1e-10)
        integrator.setInitialStepSize(100.0)

        numProp = NumericalPropagator(integrator)
        numProp.setInitialState(target_sc)
        numProp.setMu(MU)
        numProp.setSlaveMode()

        target_prop = numProp
        earth = NewtonianAttraction(MU)
        target_prop.addForceModel(earth)

        target_date = self.final_date.shiftedBy(orbit_time)
        extrapDate = self.final_date
        stepT = 100.0

        # know this is the state for final_date + time for orbit
        while extrapDate.compareTo(target_date) <= 0:
            currentState = target_prop.propagate(extrapDate)
            coord = currentState.getPVCoordinates().getPosition()
            self.target_px.append(coord.getX())
            self.target_py.append(coord.getY())
            self.target_pz.append(coord.getZ())
            extrapDate = extrapDate.shiftedBy(stepT)

    def oedot_plots(self):
        oe_params = ('sma', 'e_x', 'e_y', 'h_x', 'h_y')

        oedot = [self.adot_orbit, self.exdot_orbit, self.eydot_orbit, self.hxdot_orbit,
                 self.hydot_orbit]

        for i in range(len(oe_params)):
            plt.subplot(3,2,i+1)
            plt.plot(oedot[i])
            plt.ylabel(oe_params[i])
        plt.title('Rate of change of OE')
        plt.show()

class OutputHandler(PythonOrekitFixedStepHandler):
    """
    Implements a custom handler for every value
    """
    def init(self, s0, t):
        """
        Initilization at every prpagation step
        :param s0: initial state (spacecraft State)
        :param t: initial time (int)
        :return:
        """
        print('Orbital Elements:')

    def handleStep(self, currentState, isLast):
        """
        Perform a step at every propagation step
        :param currentState: current spacecraft state (Spacecraft state)
        :param isLast: last step in the propagation (bool)
        :return:
        """
        o = OrbitType.KEPLERIAN.convertType(currentState.getOrbit())
        print(o.getDate())
        print('a:{:5.3f}, e:{:5.3f}, i:{:5.3f}, theta:{:5.3f}'.format(o.getA(), o.getE(),
                                                                      degrees(o.getI()), degrees(o.getLv())))
        if isLast:
            print('this was the last step ')

def main():

    from Satmind.utils import OrnsteinUhlenbeck
    year, month, day, hr, minute, sec = 2018, 8, 1, 9, 30, 00.00
    date = [year, month, day, hr, minute, sec]

    dry_mass = 500.0
    fuel_mass = 150.0
    mass = [dry_mass, fuel_mass]
    duration = 24.0 * 60.0 ** 2 * 2

    # set the sc initial state
    a = 10000.0e3  # semi major axis (m) (altitude)
    e = 0.1  # eccentricity
    i = 5.0  # inclination
    omega = 10.0  # perigee argument
    raan = 10.0  # right ascension of ascending node
    lM = 10.0  # mean anomaly
    state = [a, e, i, omega, raan, lM]

    # target state
    a_targ = 35000.0e3  # altitude
    e_targ = 0.3
    i_targ = 10.0
    omega_targ = 10.0
    raan_targ = 10.0
    lM_targ = 10.0
    state_targ = [a_targ, e_targ, i_targ, omega_targ, raan_targ, lM_targ]
    stepT = 500.0

    env = OrekitEnv(state, state_targ, date, duration, mass, stepT)
    env.render_target()

    reward = []
    s = env.reset()
    F_r, F_s, F_w = 0.0, 0.0, -3.0
    while env._extrap_Date.compareTo(env.final_date) <= 0:
        thrust_mag = np.array([F_r, F_s, F_w])
        position, r, done = env.step(thrust_mag)
        reward.append(r)

    plt.subplot(2,1,1)
    plt.plot(env.hx_orbit)
    plt.ylabel('hx')
    plt.subplot(2,1,2)
    plt.plot(env.hy_orbit)
    plt.ylabel('hy')
    plt.xlabel('Mission step 500 seconds per iteration')
    plt.tight_layout()
    plt.savefig('seesaw.pdf')
    plt.show()
    # env.render_plots(save=True, show=True, episode=4)
    # env.oedot_plots()
    # print(f'days: {(env._currentDate.durationFrom(env.final_date)/3600)/24}')
    # print(f'Done\ninc:: {env._currentOrbit.getI()} \n=====')
    # print(f'diff:   a: {(env.r_target_state[0] - env._currentOrbit.getA())/1e3} km')

if __name__ == '__main__':
    main()

ext: py

max_stars_repo_name: dkolosa/Satmind

max_issues_repo_name: dkolosa/Satmind

max_forks_repo_name: dkolosa/Satmind

reason: import numpy


text:
# -*- coding: utf-8 -*-
"""predicting_covid19.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1VjqxWNcP2N4wy_NXfLt9vUToB2YXt6Ec

# Predicting Covid19 from CT Scan Images #
The patients that visit a health care centre can be diagnosed into three categries: patients who have covid19, patients who have pneumonia and healthy patients. In this project, we want to identify these patients. We use a  GPU runtime type for efficiency.

# Packages #
First, let's run the cell below to import all the packages that we will need.
- [numpy](https://numpy.org/) is the fundamental package for scientific computing with Python.\n",
- [Keras](https://keras.io) is a deep learning API written in Python, running on top of the machine learning platform TensorFlow. It was developed with a focus on enabling fast experimentation.
- [Matplotlib](https://matplotlib.org/) is a comprehensive library for creating static, animated, and interactive visualizations in Python. It is used here for plotting images.
- The [sklearn.metrics](https://scikit-learn.org/stable/modules/classes.html#module-sklearn.metrics) module includes score functions, performance metrics and pairwise metrics and distance computations.
"""

import numpy as np
from tensorflow.keras import *
import matplotlib.pyplot as plt
from sklearn.metrics import classification_report
from keras.preprocessing.image import load_img

"""# Loading Dataset of Images #
We have a train and test dataset of the ct scan images of people's lungs in the health health center. Here we load the dataset from the drive by first copying it and the unzipping the file in the current directory.
"""

!cp "/content/drive/MyDrive/ColabNotebooks/xray.zip" .

!unzip xray.zip

"""# Reading and Preprocessing #
Here are the steps we will follow in this section:
- Reading the train and test datasets
- Convert images to grayscale
- Change the size of the images to 80x80
- Print the number of read data points and the number of the categories
- Setting the batch size to 32 (for the following sections to use)

"""

def loadDataset(train_address="/content/Data/train", test_address="/content/Data/test"
                , target_size = (80,80), color_mode="grayscale"
                , classes= ["COVID19", "NORMAL", "PNEUMONIA"], batch_size = 32, rescale=1):

  dataGenerator = preprocessing.image.ImageDataGenerator(rescale=rescale, featurewise_std_normalization=True, samplewise_std_normalization=True) #will be ready to preprocess images
  classes=['COVID19', 'NORMAL', 'PNEUMONIA']
  dataTrain = dataGenerator.flow_from_directory(train_address
                                              , target_size=target_size
                                              , color_mode=color_mode
                                              , classes=classes
                                              , batch_size=batch_size)

  dataTest = dataGenerator.flow_from_directory(test_address
                                              , target_size=target_size
                                              , color_mode=color_mode
                                              , classes=classes
                                              , batch_size=batch_size)
  return dataTrain, dataTest

train, test = loadDataset()

"""# Data Visualization #
In order to have an intuitive understanding of our datasets, we will do the following:

 1.Display an image of every class of patients

 2.Draw a histogram of the number of patients in each class

For the tasks above, we loop through the given dataset, count the number of data points in each category and display the images of every new category (three categories\classes in total as mentioned above).
"""

def plot_labels(dataset):
  images, labels = dataset.next()
  types = {0:"COVID19", 1:"NORMAL", 2:"PNEUMONIA"}
  visited_labels_indices = set() #set of where the 1 index is in the type/category
  for image, label in zip(images, labels):
    type_index = backend.argmax(label).numpy()
    if type_index not in visited_labels_indices:
      visited_labels_indices.add(type_index)
      plt.title(types[type_index])
      plt.imshow(image.reshape(80,80), cmap="gray")
      plt.show()

plot_labels(train)

def plot_distribution():
  train, test = loadDataset() #we have to reload the data since it might have been changes in the previous parts
  num_of_types = train.num_classes #the number of different label types
  labels_count = {} #making a dictionary of labels count
  for type_ in range(num_of_types):
    labels_count[type_] = 0
  for label in train.classes:
    labels_count[label] += 1
  types = {0:"COVID19", 1:"NORMAL", 2:"PNEUMONIA"}
  x = list(types.values())
  print(x)
  y = [labels_count[t] for t in range(num_of_types)]
  plt.bar(x, y, color='red')

plot_distribution()

"""# Defining Model Architecture #
Here we develop a neural network with the following properties:
- The network has 4 layers
- The loss function is categorical cross entropy
- The batch size is 32
- It has 10 epochs
- All the activation functions except the last one are relu
"""

EPOCHS, BATCH_SIZE = 10, 32
def makeNN(input_shape, output_shape, activation_f="relu", optimizer=optimizers.SGD(learning_rate=0.01), loss="categorical_crossentropy"
                      , **kwargs): #make neural network
  input_ = layers.Input(shape=input_shape)
  output_ = layers.Flatten()(input_)
  output_ = layers.Dense(2048, activation=activation_f, bias_initializer='zeros', kernel_initializer='random_normal',**kwargs)(output_)
  output_ = layers.Dense(1024, activation=activation_f, bias_initializer='zeros', kernel_initializer='random_normal',**kwargs)(output_)
  output_ = layers.Dense(2048, activation=activation_f, bias_initializer='zeros', kernel_initializer='random_normal',**kwargs)(output_)
  output_ = layers.Dense(output_shape, activation="softmax", bias_initializer='zeros', kernel_initializer='random_normal',**kwargs)(output_)
  model = models.Model(inputs=input_, outputs=output_)
  model.compile(optimizer=optimizer, loss=loss, metrics=['accuracy'])
  model.summary()
  return model

our_model = makeNN((80,80,1), train.num_classes)

"""# Traing the Model on the Test Dataset #
To train the model we need to have the labels of the test dataset. Since the labels are one-hot vectors, we will first change them to their true labels (the index of 1 in their corresponding vector).

After converting the labels to integers, we can train the model, using the architechture we defined and trained above.

## Training Using ReLu Activation Function ##
As we can seen below, after 10 epochs the relu activation function results in almost perfect accuracy in the training set. The non-linearity of relu and the few number of layers, makes it efficient in our network. However, relu can not be used in networks with many layers due to the problem of vanishing gradients.

Advantages:

- It is computationally effective as it involves simpler mathematical operations than sigmoid and tanh.
- Although it looks like a linear function, it adds non-linearity to the network, making it able to learn complex patterns.
- It doesn't suffer from the vanishing gradient problem.
- It is unbounded at the positive side. Hence removing the problem of gradient saturation.
- It provides sparsity to the network, which as a result lessens the space and time complexity. Sparse networks are always faster than dense networks as there are lesser things to compute in a sparse network. Sparsity results in concise models that often have better predictive power and less overfitting/noise.

Disadvantages:

- It suffers from the dying ReLU problem. ReLU is always going to discard the negative values i.e. the deactivations by making it 0. But because of this, the gradient of these units will also become 0 and by now we all know that 0 gradient means no weight updation during backpropagation. Simply speaking, the neurons which will go to this state will stop responding to the deviation of input or the error. This, as a result, hampers the ability of the model to fit the data properly.
-  It is non-differentiable at 0.

Uses

ReLU should be used in the hidden layers. As it is computationally less expensive than sigmoid and tanh, therefore it is usually a better choice than them. It is also to be noted that ReLU is faster than both tanh and sigmoid. Also in hidden layers, at a time only a few neurons are activated, making it efficient and easy for computation.
"""

x, y = train.next()
history = our_model.fit(train, epochs=EPOCHS)

""" Getting the actual labels: to be able to test our models we need to have the actual functions. To do so, we replace each label (one-hot vector) with the index of 1 in the vector.

"""

def get_labels(dataset):
  all_labels = []
  all_batches = int(np.ceil(len(dataset.classes)/dataset.batch_size))
  for batch_index in range(all_batches):
    images, labels = dataset.next()
    for label in labels:
      all_labels.append(backend.argmax(label))
  return np.array(all_labels) #since we need a numpy array for training our model

train, test = loadDataset()
predictions=our_model.predict(train)
actual_labels = get_labels(train)
print(backend.argmax(predictions).numpy().shape)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))

def plot_acc_loss(history):
  plt.plot(history.history['accuracy'])
  plt.title('model accuracy')
  plt.ylabel('accuracy')
  plt.xlabel('epoch')
  plt.legend(['accuracy'], loc='upper left')
  plt.show()
  plt.plot(history.history['loss'])
  plt.title('model loss')
  plt.ylabel('loss')
  plt.xlabel('epoch')
  plt.legend(['loss'], loc='upper left')
  plt.show()

plot_acc_loss(history)

"""## Testing Using ReLu Activation Fuction #

"""

predictions=our_model.predict(test)
actual_labels = get_labels(test)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))

"""## Training Using Tanh Activation Function ##
 Tanh is also like logistic sigmoid but better. The range of the tanh function is from (-1 to 1). The equation for tanh is $f(x) = 2/(1 + e^{-2x})-1$.

 Advantages

 - It has similar advantages as sigmoid but better than that because it is zero centred. The output of tanh lies between -1 and 1. Hence solving one of the issues with the sigmoid.

Disadvantages

- It also has the problem of vanishing gradient but the derivatives are steeper than that of the sigmoid. Hence making the gradients stronger for tanh than sigmoid.
- As it is almost similar to sigmoid, tanh is also computationally expensive.
    Similar to sigmoid, here also the gradients saturate.
Uses

tanh is generally used in the hidden layers. The output of tanh lies between -1 and +1. Hence the mean for the hidden layers stays near 0. This makes the learning of the hidden layer converge quickly.

Uses

 Tanh is generally used in the hidden layers. The output of tanh lies between -1 and +1. Hence the mean for the hidden layers stays near 0. This makes the learning of the hidden layer converge quickly.
"""

train, test = loadDataset()
tanh_model = makeNN((80,80,1), train.num_classes, activation_f='tanh')
x, y = train.next()
history = our_model.fit(train, epochs=EPOCHS)
train, test = loadDataset()
predictions=our_model.predict(train)
actual_labels = get_labels(train)
print(backend.argmax(predictions).numpy().shape)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))
plot_acc_loss(history)

"""## Testing Using Tanh Activation Fuction #

"""

predictions=our_model.predict(test)
actual_labels = get_labels(test)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))

"""## ReLu vs Tanh ##
As we can see above, Tanh activation function is working better than ReLu on these datasets. This is beacuse tanh converges to zero more quickly in the hidden layers. Moreover, ReLu is more prune to overfitting than tanh. We can solve this problem by variously overfitting methods. Also some learning rates and other hyperparameters are more suitable for some activation functions. Although, both of these functions have the gradient vanishing problem which is not obvious becasue of having a few layers.

# Normalizing and Optimizing the Neural Network #
To normalized the data we divide the pixels in the images by 255 (the range of numbers for each pixel is 0-255) and to improve the network, try adding more layers and chaning the learning rate. There are other things we could do i.e. changing the optimizer etc, that we will do in the following sections.
"""

train, test = loadDataset(rescale=1/255.0) #normalizing
EPOCHS, BATCH_SIZE = 10, 32
def makeNN(input_shape, output_shape, activation_f="relu", optimizer=optimizers.SGD(learning_rate=0.001), loss="categorical_crossentropy"): #make neural network
  input_ = layers.Input(shape=input_shape)
  output_ = layers.Flatten()(input_)
  output_ = layers.Dense(2048, activation=activation_f, bias_initializer='zeros', kernel_initializer='random_normal')(output_)
  output_ = layers.Dense(1024, activation=activation_f, bias_initializer='zeros', kernel_initializer='random_normal')(output_)
  output_ = layers.Dense(1024, activation=activation_f, bias_initializer='zeros', kernel_initializer='random_normal')(output_)
  output_ = layers.Dense(2048, activation=activation_f, bias_initializer='zeros', kernel_initializer='random_normal')(output_)
  output_ = layers.Dense(output_shape, activation="softmax", bias_initializer='zeros', kernel_initializer='random_normal')(output_)
  model = models.Model(inputs=input_, outputs=output_)
  model.compile(optimizer=optimizer, loss=loss, metrics=['accuracy'])
  model.summary()
  return model

our_model = makeNN((80,80,1), train.num_classes, activation_f="tanh")

"""## Fitting and Testing ##"""

x, y = train.next()
history = our_model.fit(train, epochs=EPOCHS)
predictions=our_model.predict(test)
actual_labels = get_labels(test)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))

"""# Other Optimizers #
Optimizer helps us to minimize an error function(loss function)or to maximize the efficiency of production .In this section we will apply two other optimizers: momentum and Adam

## Momentum : What is it and why is it useful? ##
Momentum is not soley an optimizer, it is an improvement on Stochastic Gradient Descent. When looking at noisy time series data such as your training/validation error graphs in tensorboard, you may often notice that the raw values are often quite noisy. Quite often, you are able to see a trend in the graph. These trends often become more obvious when you add some smoothing to the raw graph values.

Exponential Smoothing is one of the simplest ways to add smoothing to your data.

Exponential Smoothing Formula:
$$ s_0 = x_0 s_t = ax_t + (1-\alpha)s_{t-1}, t>0, 0<= \alpha <=1$$
In the above equation, momentum specifies the amount of smoothing we want. A typical value for momentum is .9. From this equation, we can see that the value at time step t takes into account previous values from time steps (t — 1) and before. The weighting of previous time steps drops off exponentially so the most recent time step has the greatest impact.The momentum optimizer uses this same equation when optimizing a loss function. In this case, the Exponential Smoothing is called Momentum and it allows the optimizer to maintain the velocity and trajectory from previous time steps to affect the present.

## Is Momentum always a good optimizer? ##
 Unfortunately, while improvements to the momentum algorithm do exist, they all run into a certain, critical, almost inescapable lower bound. That means, we can not always apply momentum and get better results.

## Implementing Momentum Optimizer on Our Model ##
 From the models below, we can see that the method of  early stopping applies for momentum too. When it is increased from 0.5 to 0.9 the model is improved but for momentum 0.99 the model has deteriorated.

 Momentum = 0.5
"""

our_model = makeNN((80,80,1), train.num_classes, activation_f="tanh", optimizer=optimizers.SGD(learning_rate=0.001, momentum = 0.5))
x, y = train.next()
history = our_model.fit(train, epochs=EPOCHS)
predictions=our_model.predict(test)
actual_labels = get_labels(test)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))

""" Momentum = 0.9

 This momentum is giving us an almost perfect model.
"""

our_model = makeNN((80,80,1), train.num_classes, activation_f="tanh", optimizer=optimizers.SGD(learning_rate=0.001, momentum = 0.9))
history = our_model.fit(train, epochs=EPOCHS)
predictions=our_model.predict(test)
actual_labels = get_labels(test)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))

""" Momentum = 0.99"""

our_model = makeNN((80,80,1), train.num_classes, activation_f="tanh", optimizer=optimizers.SGD(learning_rate=0.001, momentum = 0.99))
history = our_model.fit(train, epochs=EPOCHS)
predictions=our_model.predict(test)
actual_labels = get_labels(test)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))

"""## Adam: What is the Adam optimization algorithm? ##
Adam is an optimization algorithm that can be used instead of the classical stochastic gradient descent procedure to update network weights iterative based in training data.

Adam was presented by Diederik Kingma from OpenAI and Jimmy Ba from the University of Toronto in their 2015 ICLR paper (poster) titled “Adam: A Method for Stochastic Optimization“. I will quote liberally from their paper in this post, unless stated otherwise.

The algorithm is called Adam. It is not an acronym and is not written as “ADAM”.

    … the name Adam is derived from adaptive moment estimation.

When introducing the algorithm, the authors list the attractive benefits of using Adam on non-convex optimization problems, as follows:

-    Straightforward to implement.
-    Computationally efficient.
-    Little memory requirements.
-    Invariant to diagonal rescale of the gradients.
-    Well suited for problems that are large in terms of data and/or parameters.
-    Appropriate for non-stationary objectives.
-    Appropriate for problems with very noisy/or sparse gradients.
-    Hyper-parameters have intuitive interpretation and typically require little tuning.

### How Does Adam Work? ###

Adam is different to classical stochastic gradient descent.

Stochastic gradient descent maintains a single learning rate (termed alpha) for all weight updates and the learning rate does not change during training.

A learning rate is maintained for each network weight (parameter) and separately adapted as learning unfolds.

The method computes individual adaptive learning rates for different parameters from estimates of first and second moments of the gradients.

The authors describe Adam as combining the advantages of two other extensions of stochastic gradient descent. Specifically:

Instead of adapting the parameter learning rates based on the average first moment (the mean), Adam also makes use of the average of the second moments of the gradients (the uncentered variance).

Specifically, the algorithm calculates an exponential moving average of the gradient and the squared gradient, and the parameters beta1 and beta2 control the decay rates of these moving averages.

The initial value of the moving averages and beta1 and beta2 values close to 1.0 (recommended) result in a bias of moment estimates towards zero. This bias is overcome by first calculating the biased estimates before then calculating bias-corrected estimates.

Adam is Effective

Adam is a popular algorithm in the field of deep learning because it achieves good results fast.

Empirical results demonstrate that Adam works well in practice and compares favorably to other stochastic optimization methods.

Adam VS SGD

To better understand the paper’s implications, it is necessary to first look at the pros and cons of popular optimization algorithms Adam and SGD.

Gradient descent is the most common method used to optimize deep learning networks. First proposed in the 1950s, the technique can update each parameter of a model, observe how a change would affect the objective function, choose a direction that would lower the error rate, and continue iterating until the objective function converges to the minimum.

SGD is a variant of gradient descent. Instead of performing computations on the whole dataset — which is redundant and inefficient — SGD only computes on a small subset or random selection of data examples. SGD produces the same performance as regular gradient descent when the learning rate is low.

Essentially Adam is an algorithm for gradient-based optimization of stochastic objective functions. It combines the advantages of two SGD extensions — Root Mean Square Propagation (RMSProp) and Adaptive Gradient Algorithm (AdaGrad) — and computes individual adaptive learning rates for different parameters.

Despite the widespread popularity of Adam, recent research papers have noted that it can fail to converge to an optimal solution under specific settings. The paper Improving Generalization Performance by Switching from Adam to SGD demonstrates that adaptive optimization techniques such as Adam generalize poorly compared to SGD. This has prompted some researchers to explore new techniques that may improve on Adam
"""

our_model = makeNN((80,80,1), train.num_classes, activation_f="tanh", optimizer=optimizers.Adam(learning_rate=0.001,beta_1=0.9,
                                                                                beta_2=0.999,epsilon=1e-07,
                                                                                amsgrad=False,name="Adam"))

"""## Setting EPOCHS to 20 ##
From some epoch on, the efficiency of the model on the training set does not change which indicates the bound of the optimization scheme.
"""

EPOCHS = 20
history = our_model.fit(train, epochs=EPOCHS)
predictions=our_model.predict(test)
actual_labels = get_labels(test)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))

"""## Why are epochs necesary? ##
Generally whenever we want to optimize we use gradient descent. Gradient descent has a parameter called learning rate. In one iteration alone we can not guarantee that the gradient descent algorithm would converge to a local minima with the specified learning rate. That is the reason why we iterate again for the gradient descent to converge better.
## Do all problems need multiple epochs? ##
If the optimization is straight forward in a problem, for example if we have a linear equation to optimize, we could reach the solution with even a single pass through the dataset. That is why multiple epochs are not necessary for all problems.
## Is a large number of EPOCHS always beneficial for the problem? ##
Not necessarily. The correct answer is the number of epochs is not that significant. more important is the validation and training error. As long as these two error keeps dropping, training should continue.

For instance, if the validation error starts increasing that might be an indication of overfitting.

You should set the number of epochs as high as possible and terminate the training when validation error start increasing

# Using MSE Loss Function #
Here we want to replace categorical cross entropy with the MSE loss function. MSE is the sum of squared distances between our target variable and predicted values. As we can see below, categorical cross entropy, results in a better (more accurate model) for classification.
"""

our_model = makeNN((80,80,1), train.num_classes, activation_f="tanh", optimizer='sgd', loss=losses.MeanSquaredError())
EPOCHS = 20
history = our_model.fit(train, epochs=EPOCHS)
predictions=our_model.predict(test)
actual_labels = get_labels(test)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))

"""### Why is MSE not suitable for classification? ###
Firstly, using MSE means that we assume that the underlying data has been generated from a normal distribution (a bell-shaped curve). In Bayesian terms this means we assume a Gaussian prior. While in reality, a dataset that can be classified into categories is not from a normal distribution but a Bernoulli distribution.

If we use maximum likelihood estimation(MLE), assuming that the data is from a normal distribution(a wrong assumption, by the way), we get the MSE as a Cost function for optimizing our model.

Secondly, the MSE function is non-convex for binary classification. In simple terms, if a binary classification model is trained with MSE Cost function, it is not guaranteed to minimize the Cost function. This is because MSE function expects real-valued inputs in range(-∞, ∞), while binary classification models output probabilities in range(0,1) through the sigmoid/logistic function.

When the MSE function is passed a value that is unbounded a nice U-shaped (convex) curve is the result where there is a clear minimum point at the target value (y). On the other hand, when a bounded value from a Sigmoid function is passed to the MSE function the result is not convex; on one side the function is concave while on the other side the function convex and no clear minimum point. So, if by accident a binary classification neural network is initialized with weights which are large in magnitude such that it lands on the concave part of the MSE Cost function gradient descent will not work and consequently, weights may not update or improve very slowly.

On a final note, MSE is a good choice for a Cost function when we are doing Linear Regression (i.e fitting a line through data for extrapolation). In the absence of any knowledge of how the data is distributed assuming normal/gaussian distribution is perfectly reasonable.

# Regularization #
One of the major aspects of training your machine learning model is avoiding overfitting. The model will have a low accuracy if it is overfitting. This happens because your model is trying too hard to capture the noise in your training dataset. By noise we mean the data points that don’t really represent the true properties of your data, but random chance. Learning such data points, makes your model more flexible, at the risk of overfitting.

A modern approach to reducing generalization error is to use a larger model that may be required to use regularization during training that keeps the weights of the model small. These techniques not only reduce overfitting, but they can also lead to faster optimization of the model and better overall performance.

There are two ways to approach an overfit model:

-  Reduce overfitting by training the network on more examples.
-   Reduce overfitting by changing the complexity of the network.

Below is a list of five of the most common additional regularization methods.

- 1.Activity Regularization: Penalize the model during training base on the magnitude of the activations.
- 2.Weight Constraint: Constrain the magnitude of weights to be within a range or below a limit.
- 3.Dropout: Probabilistically remove inputs during training.
- 4.Noise: Add statistical noise to inputs during training.
- 5.Early Stopping: Monitor model performance on a validation set and stop training when performance degrades.

## L2 Regularization ##
We reimplemnt our model with the L2 Regularization in Keras with a value of 0.0001. L1 regularization adds the penalty term in cost function by adding the absolute value of weight(Wj) parameters, while L2 regularization adds the squared value of weights in the cost function. This improves the accuracy, as seen below.

Loss function with L1 regularization:
$$ Loss = Error(y,\hat y) + \lambda \sum_{i=1}^{N} |w_i|$$

Loss function with L2 regularization:
$$ Loss = Error(y,\hat y) + \lambda \sum_{i=1}^{N} w_i^{2}$$
"""

our_model = makeNN((80,80,1), train.num_classes, activation_f="tanh", optimizer=optimizers.SGD(learning_rate=0.001, momentum = 0.99)
                            ,kernel_regularizer=regularizers.l2(1e-4), bias_regularizer=regularizers.l2(1e-4))
history = our_model.fit(train, epochs=EPOCHS)
predictions=our_model.predict(test)
actual_labels = get_labels(test)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))

"""## Dropout ##
Ensembles of neural networks with different model configurations are known to reduce overfitting, but require the additional computational expense of training and maintaining multiple models.

Dropout is a regularization method that approximates training a large number of neural networks with different architectures in parallel.

During training, some number of layer outputs are randomly ignored or “dropped out.” This has the effect of making the layer look-like and be treated-like a layer with a different number of nodes and connectivity to the prior layer. In effect, each update to a layer during training is performed with a different “view” of the configured layer.

### How to Dropout? ###

Dropout is implemented per-layer in a neural network.

It can be used with most types of layers, such as dense fully connected layers, convolutional layers, and recurrent layers such as the long short-term memory network layer.

Dropout may be implemented on any or all hidden layers in the network as well as the visible or input layer. It is not used on the output layer.

A new hyperparameter is introduced that specifies the probability at which outputs of the layer are dropped out, or inversely, the probability at which outputs of the layer are retained.

A common value is a probability of 0.5 for retaining the output of each node in a hidden layer and a value close to 1.0, such as 0.8, for retaining inputs from the visible layer.

#Applying Dropout #
Applying dropout as below, results in a better accuracy as we expected based on the explanations above.
"""

our_model.add(Dropout(0.1))
history = our_model.fit(train, epochs=EPOCHS)
predictions=our_model.predict(test)
actual_labels = get_labels(test)
print(classification_report(actual_labels, backend.argmax(predictions).numpy()))

"""# References #
- [Keras](https://keras.io/)
- [Tensorflow Quick Start for Beginners](https://www.tensorflow.org/tutorials/quickstart/beginner)
- [Gentle Introduction to the Adam Optimization Algorithm for Deep Learning](https://machinelearningmastery.com/adam-optimization-algorithm-for-deep-learning/#:~:text=Adam%20is%20an%20optimization%20algorithm,iterative%20based%20in%20training%20data.&text=The%20algorithm%20is%20called%20Adam,derived%20from%20adaptive%20moment%20estimation.)
- [Why Momentum Really Works](https://distill.pub/2017/momentum/)
- [‘Fast as Adam & Good as SGD’ ](https://medium.com/syncedreview/iclr-2019-fast-as-adam-good-as-sgd-new-optimizer-has-both-78e37e8f9a34)
- [why-use-the-momentum-optimizer-with-minimal-code](https://medium.com/analytics-vidhya/why-use-the-momentum-optimizer-with-minimal-code-example-8f5d93c33a53)
-[What is epoch and How to choose the correct number of epoch](https://medium.com/@upendravijay2what-is-epoch-and-how-to-choose-the-correct-number-of-epoch-d170656adaaf)
- [Why Using Mean Squared Error(MSE) Cost Function for Binary Classification is a Bad Idea?](https://towardsdatascience.com/why-using-mean-squared-error-mse-cost-function-for-binary-classification-is-a-bad-idea-933089e90df7)
- [Regularization in Machine Learning](https://towardsdatascience.com/regularization-in-machine-learning-76441ddcf99a)
- [How to Avoid Overfitting in Deep Learning Neural Networks](https://machinelearningmastery.com/introduction-to-regularization-to-reduce-overfitting-and-improve-generalization-error/)
- [L1 vs L2 Regularization: The intuitive difference](https://medium.com/analytics-vidhya/l1-vs-l2-regularization-which-is-better-d01068e6658c)
- [A Gentle Introduction to Dropout for Regularizing Deep Neural Networks](https://machinelearningmastery.com/dropout-for-regularizing-deep-neural-networks/)
"""

max_stars_repo_name: pegimdp/AI-projects

max_issues_repo_name: pegimdp/AI-projects

max_forks_repo_name: pegimdp/AI-projects


text:
import numpy as np
from skimage import transform
from skimage.feature import match_template
from skimage.color import rgb2gray

def find_max_correlation(image, template):
    """Find the location of the maximum correlation between the template and
    the image.

    Parameters
    ----------

    image: M, N  ndarray
        The image to search within for the template. Colour images must first
        be converted to a channel.
    template: m, n ndarray
        The template image to search for. Must be single channel.

    Returns
    -------
    (x,y): tuple of x and y coordinates in the image.

    """
    score = match_template(image, template)
    ij = np.unravel_index(np.argmax(score), score.shape)
    x, y = ij[::-1]
    return x, y

def extract_patches(image, windows, pad=0):
    """Extract grayscale patches from the image at the given locations.

    Parameters
    ----------

    image: M,N,3 ndarray
        The three channel image to extract the windows from.
    windows: n_windows*2,2 ndarray
        X,Y coordinates of starting and finishing points of each rectangular
        window.
    pad: integer (default=0)
        The amount to pad the window size from the given locations.

    Returns
    -------
    patches: list of ndarrays
        The patches extracted from the image.

    """
    patches = []
    coords = []
    max_rows, max_cols = image.shape[:2]
    n_windows = int(windows.shape[0] / 2)

    for i in range(n_windows):
        rows = windows[i*2:(i*2 + 2), 0] + [-pad, pad]
        cols = windows[i*2:(i*2 + 2), 1] + [-pad, pad]
        rows[rows < 0] = 0
        rows[rows > max_rows] = max_rows
        cols[cols < 0] = 0
        cols[cols > max_cols] = max_cols
        this_patch = image[rows[0]:rows[1], cols[0]:cols[1]]
        patches.append(rgb2gray(this_patch))
        coords.append(np.vstack((rows, cols)).T)

    if pad > 0:
        return patches, np.vstack(coords)
    else:
        return patches

class Registrator(object):
    """Transform a reference image to match a baseline image.

    Parameters
    ----------

    windows: (n_windows*2),2 ndarray
        x,y coordinates for focal points.
    base_image: M,N,3 ndarray
        Baseline image to match other images to.
    pad: integer (default=400)
        Size of search space outside the defined window in images to be
        registered.

    """
    def __init__(self, windows, base_image, pad=200):
        self.windows = windows
        self.pad = pad
        self.templates = extract_patches(base_image, windows)

        # Determine suitable transform from the number of points
        n_points = len(self.templates)

        n_points = min(n_points, 4) # uses projective transformation for >= 4 points

        tform_types = ['translate', 'similarity', 'affine', 'projective']
        self.tform_type = tform_types[n_points - 1]

    def match(self, image):
        """Matches an image to the baseline image.

        Parameters
        ----------

        image: M,N,3 ndarray
            Image to be matched.

        Returns
        -------

        transformed_image: M,N,3 ndarray
            The input image transformed to match the baseline image at the
            selected points.

        """
        search_windows, search_coords = extract_patches(image,
                                                        self.windows,
                                                        pad=self.pad)

        shifts = []

        for window, template in zip(search_windows, self.templates):
            shifts.append(find_max_correlation(window, template))

        shifts = np.vstack(shifts)

        # Convert coordinates from padded windows to absolute position
        delta = search_coords[::2, [1, 0]] - self.windows[::2, [1, 0]]
        points1 = self.windows[::2, [1, 0]]
        points2 = points1 - shifts - delta

        if self.tform_type == 'translate':
            x, y = (-shifts - delta).ravel()
            match_tform = transform.SimilarityTransform(translation=(-x, -y))
        else:
            match_tform = transform.estimate_transform(self.tform_type,
                                                       points2,
                                                       points1)
        return (transform.warp(image, match_tform), match_tform)

    def __call__(self, image):
        return self.match(image)

max_stars_repo_name: SamHames/aperturesynth

max_issues_repo_name: SamHames/aperturesynth

max_forks_repo_name: SamHames/aperturesynth


text:
import keras
from keras.datasets import cifar10
from keras.preprocessing.image import ImageDataGenerator
from keras.models import Sequential
from keras.layers import Dense, Dropout, Activation, Flatten
from keras.layers import Conv2D, MaxPooling2D
from keras.models import model_from_json
import cv2
from scipy.misc import imresize
imgDimension=100

json_file = open('model (1).json', 'r')
loaded_model_json = json_file.read()
json_file.close()
loaded_model = model_from_json(loaded_model_json)
loaded_model.load_weights("model.h5")
print("Loaded model from disk")
loaded_model.compile(loss='categorical_crossentropy', optimizer='rmsprop', metrics=['accuracy'])
img1=r"C:\Users\PIYUSH\Desktop\Images\TCImages\TC001_39.png"
img2=r"C:\Users\PIYUSH\Desktop\Images\TCImages\TS001_11.png"
img1 = cv2.imread(img1)

img1=imresize(img1,[imgDimension,imgDimension,3])
img1 = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY)
img1=img1.reshape(1,100,100,1)
score1 = loaded_model.predict_classes(img1, verbose=1)
#score2 = loaded_model.predict_classes(img2, verbose=1)
if score1==[0]:
    print("NOT AUTISTIC")

reason: from scipy